[1]UC Berkeley [2]Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LIG, 38000
Grenoble, France [3]Singapore University of Technology and Design

# A Quadratic Speedup in Finding Nash Equilibria of Quantum Zero-Sum Games

Francisca Vasconcelos<EMAIL_ADDRESS>Emmanouil-Vasileios Vlatakis-
Gkaragkounis<EMAIL_ADDRESS>Panayotis Mertikopoulos  Georgios
Piliouras  Michael I. Jordan

###### Abstract

Recent developments in domains such as non-local games, quantum interactive
proofs, and quantum generative adversarial networks have renewed interest in
quantum game theory and, specifically, quantum zero-sum games. Central to
classical game theory is the efficient algorithmic computation of Nash
equilibria, which represent optimal strategies for both players. In 2008, Jain
and Watrous proposed the first classical algorithm for computing equilibria in
quantum zero-sum games using the Matrix Multiplicative Weight Updates (MMWU)
method to achieve a convergence rate of $\mathcal{O}(d/\epsilon^{2})$
iterations to $\epsilon$-Nash equilibria in the $4^{d}$-dimensional
spectraplex. In this work, we propose a hierarchy of quantum optimization
algorithms that generalize MMWU via an extra-gradient mechanism. Notably,
within this proposed hierarchy, we introduce the Optimistic Matrix
Multiplicative Weights Update (OMMWU) algorithm and establish its average-
iterate convergence complexity as $\mathcal{O}(d/\epsilon)$ iterations to
$\epsilon$-Nash equilibria. This quadratic speed-up relative to Jain and
Watrous’ original algorithm sets a new benchmark for computing $\epsilon$-Nash
equilibria in quantum zero-sum games.

###### Contents

  1. 1 Introduction
    1. 1.1 Motivation
    2. 1.2 Prior Work
    3. 1.3 Our Contributions and Methodology
    4. 1.4 Paper Organization
  2. 2 Preliminaries
    1. 2.1 Quantum Zero-Sum Games
    2. 2.2 Quantum Nash Equilibria
    3. 2.3 Convexity, Smoothness, & Duality
    4. 2.4 Bregman Divergence & Mirror Map
    5. 2.5 Mirror & Proximal Steps
      1. 2.5.1 Mirror Step
      2. 2.5.2 Proximal Step
      3. 2.5.3 Regularizers
  3. 3 Algorithms for Quantum Zero-Sum Games
    1. 3.1 Matrix Dual Averaging (MDA) Method
    2. 3.2 Matrix Mirror Prox (MMP)
    3. 3.3 Optimistic Matrix Mirror-Prox (OMMP) Methods
  4. 4 Convergence Analysis
    1. 4.1 Monotonicity of $\mathcal{F}$
    2. 4.2 Error Decomposition
    3. 4.3 Optimality Analysis of $(A_{t})$
    4. 4.4 Optimality Analysis of $(C_{t})$
    5. 4.5 An Upper Bound for $(B_{t})$
    6. 4.6 Average-Iterate Error Upper Bound
    7. 4.7 Regularizer-Induced Dimensionality Dependence
  5. 5 Conclusions
  6. 6 Acknowledgements
  7. A Complex Matrix Differentiation and Gradients
  8. B Mirror & Proximal Maps
  9. C Properties of Finite-Valued Zero-Sum Quantum Games
    1. C.1 Properties of Spectraplexes
    2. C.2 Properties of $\mathcal{F}$
  10. D Convergence Analysis Proofs

##  1 Introduction

###  1.1 Motivation

Ever since von Neumann’s groundbreaking work in the 1920’s [1], game theory
has become a foundational pillar of modern mathematics, economics, and
computer science. The central solution concept in non-cooperative game theory
is that of a _Nash equilibrium_ , i.e., a strategy of joint play in which no
single player can benefit from a unilateral deviation [2]. Despite its
mathematical significance, the computational perspective on general Nash
equilibria is murkier, since even approximating such an equilibrium is a PPAD-
complete problem for general classical games. On the other hand, Nash
equilibria in two-player zero-sum games (i.e., where one player’s gain is the
other’s loss) can be computed in polynomial time, and they have found
practical applications in a wide variety of domains, from resource allocation
problems [3], to political strategy [4], and the training of machine learning
models [5, 6, 7].

When the players’ actions are the states of a quantum system, classical game
theory no longer applies, and one must consider _quantum games_. These games
constitute a natural framework for the study of the exchange and processing of
quantum information, resulting in various useful applications to theoretical,
computational, and cryptographic settings. The original interest in quantum
games dates back to seminal work by Bell [8] and CSHS [9], explaining the
Einstein-Podolsky-Rosen (EPR) Paradox. Since then, substantial work on quantum
“non-local games” [10, 11, 12, 13, 14] has lead to novel methods for testing
whether two or more spatially separated, non-communicating quantum devices can
generate correlations not reproducible by any pair of classical
devices—resulting in entanglement-based tests of quantum advantage. Much of
this work culminated in establishing MIP∗=RE [15], which, among many notable
results, demonstrated that there exists an efficient reduction from the
Halting Problem to deciding whether a two-player non-local game has entangled
value at most 1 or $\frac{1}{2}$. Furthermore, it was recently proven that
non-local games can be compiled into single-prover interactive games [16]. In
general, quantum interactive proofs with competing provers can be modeled as
competitive quantum refereed games [17, 18, 19] and multi-prover quantum
interactive proofs [20, 21, 22, 23] can be modeled as cooperative quantum
games. Finally, there has been substantial work on quantum coin-flipping [24,
25, 26, 27], a game model where two players directly exchange quantum
information, studying how two parties with competing interests can carry out a
fair coin flip across a quantum communication channel.

In this work, we consider the specific class of quantum _zero-sum_ games, in
which two players are in direct competition and do not share entanglement. One
of the earliest works in quantum game theory [28] studied a “matching
pennies”-type zero-sum game to prove that quantum strategies are at least as
good as and, in some cases, can outperform classical strategies. Further work
on quantum zero-sum games led to a proof that the quantum complexity class
QRG(1), of problems having one-turned quantum refereed games, is contained in
PSPACE [29]. Finally, in the context of quantum machine learning, there has
been substantial recent interest in quantum generative adversarial networks
(QGANs) [30, 31], for which the training of the competing generator and
discriminator networks can be modeled as a quantum zero-sum game.

As in the classical setting, finding Nash equilibria of _general_ quantum
games is computationally prohibitive—PPAD-complete to be exact [32]. However,
the situation is more favorable for quantum _zero-sum_ games. As demonstrated
by Jain and Watrous [29], Nash equilibria of $d$-qubit two-player quantum
zero-sum games can be calculated to $\epsilon$-accuracy in polynomial time and
linear dependence on the number of qubits (logarithmic in the spectraplex
dimensionality), i.e., $\mathcal{O}(d/\epsilon^{2})$, even with payoff-based
information on the players’ side [33]. In view of this, equilibrium strategies
in quantum zero-sum games can in principle be computed algorithmically, thus
providing the required guarantees of implementability for proof schemes or
cryptographic ciphers that rely on the computation of quantum Nash equilibria.
Our goal is to pursue this objective further, bringing tools from online
learning and computational learning theory to improve the design of algorithms
for the computation of approximate Nash equilibria in quantum zero-sum games.

###  1.2 Prior Work

Our work builds upon the work of Jain and Watrous [29] for non-interactive
quantum zero-sum games. To provide some context, [29] leveraged feedback in
the form of quantum channels (governed by superoperators) to introduce the
Matrix Multiplicative Weights Update (MMWU) algorithm, a two-player variant of
the matrix exponentiated gradient algorithm of Tsuda et al. [34], which is
itself a special case of Mirror Descent [35, 36]. Beyond the setting of
quantum zero-sum games, variants of the MMWU algorithm have also been used in
many important applications such as proving QIP=PSPACE [37], solving SDPs
[38], finding balanced separators [39], enhancing spectral sparsification
[40], covariance optimization tasks [41, 42, 43], and matrix learning [44].
For the quantum zero-sum games of interest in this work, despite recent work
on no-regret learning dynamics for zero-sum and general quantum games [33,
45], to the best of our knowledge, the $\mathcal{O}(d/\epsilon^{2})$
convergence rate of the MMWU algorithm is the state of the art for finding
$\epsilon$-approximate Nash equilibria. In this work, we ask:

Can we improve on the performance of MMWU in quantum zero-sum games?

More precisely, our paper’s goal and main contribution is to provide an
accelerated version of the MMWU algorithm which achieves $\epsilon$-accuracy
in linear, $O(d/\epsilon)$ time, instead of $O(d/\epsilon^{2})$.

To achieve this quadratic speedup, our point of departure is the literature on
classical _finite_ games. The classical version of the MMWU algorithm—known,
among other names, as the Multiplicative Weights Update (MWU) method
[46]—similarly achieves an $\mathcal{O}_{d}(1/\epsilon^{2})$ iteration
complexity.111Hereon, we will use the $\mathcal{O}_{d}(\cdot)$ notation to
hide the dimensionality dependence of the iteration complexity. There will be
explicit discussion of dimensionality dependence in the quantum setting in
Section 4.7. However, Nemirovski [47] and Auslender & Teboulle [48] showed
that, by a suitable modification of the MWU method—now commonly known as the
_mirror-prox_ family of algorithms—it is possible to accelerate this rate to
$\mathcal{O}_{d}(1/\epsilon)$, a rate which was recently shown by Ouyang and
Xu [49] to be order-optimal in the setting of min-max convex-concave problems.

The main algorithmic insight of the mirror-prox template is the combination of
an extra-gradient step in the spirit of Korpelevich [50] with an iterate-
averaging mechanism à la Polyak & Juditsky [51]. Coupled with the smoothness
of the underlying objective, this combination yields a
$\mathcal{O}(1/\epsilon)$ convergence rate, albeit with a possibly suboptimal
dependence on the dimension of the problem. In the case of the simplex—which
is the state space of classical, finite games—the dependence on the dimension
can be improved dramatically by instantiating the mirror-prox with an entropic
regularizer, which essentially boils down to a variant of MWU “with advice”
[47]. The only notable drawback of the resulting algorithm is that it requires
not one, but two oracle queries per iteration, which doubles the computation
cost, and requires additional coordination / communication from the players’
side. Building on an original idea by Popov [52], this “cost doubling” issue
was mitigated by the so-called _optimistic mirror descent_ (OMD) proposal of
Raklhin and Sridharan [53], which optimistically reused past gradient
information as a surrogate for the extra-gradient step of the mirror-prox
algorithm. In this way, when applied to mixed extensions of classical, finite
games, OMD ultimately achieves an order-optimal convergence rate in terms of
both $\epsilon$ and $d$, all the while using a single oracle call per
iteration. We summarize this hierarchy of optimization methods and the
relevant trade-offs in Figure 1.

###  1.3 Our Contributions and Methodology

Our work focuses on developing an accelerated algorithm with similar, linear-
time performance guarantees in the _quantum_ setting. The main contributions
of this work are thus twofold: a) we propose the _Optimistic Matrix
Multiplicative Weights Update (OMMWU) algorithm_ (Algorithm 1) for finding
$\epsilon$-Nash equilibria in quantum zero-sum games; and b) we show that
OMMWU achieves the following iteration complexity guarantees.

######  Corollary 1 (OMMWU Iteration Complexity)

In a $4^{d}$-dimensional spectraplex, OMMWU computes average-iterate
$\epsilon$-Nash equilibria in $\mathcal{O}(d/{\epsilon})$ iterations.

In order to achieve this result, we forge a link between quantum zero-sum
games and established works in classical game theory and optimization. In
addition to OMMWU—which surpasses the $\mathcal{O}(d/\epsilon^{2})$ rate of
MMWU with a single gradient call per iteration—we also propose a hierarchical
family of algorithms and techniques that can be applied to various
semidefinite programming problems (quantum or otherwise) and thus may be of
independent interest.

Methodology. In terms of methodology, we move away from the _channel-based_
view of quantum zero-sum games in [29] and instead adopt a _gradient-based_
perspective, similar to prior work on learning in general quantum games [54]
(see also [45] connecting the same online learning, regret minimization
perspective to new classes of quantum correlated equilibria). In this
framework, Alice and Bob’s feedback is characterized not by traditional
quantum superoperators, but rather by mathematical gradient operators, which
comprise the game’s gradient feedback operator $\mathcal{F}$. We prove several
important properties of $\mathcal{F}$, namely that Lipschitz continuity and
monotonicity. In doing so, we conclude that solving a quantum zero-sum game to
a desired accuracy essentially boils down to solving a smooth, semidefinite
convex-concave problem with first-order oracle information.

Algorithm 1 Optimistic Matrix Multiplicative Weights Update (OMMWU)

Accuracy Parameter: $\epsilon$

Regularization Function: $h:\\{\mathcal{A},\mathcal{B}\\}\to\mathbb{R}$

Strong Convexity Parameter (of $h$): $\mu_{h}$

Diameter (of $\mathcal{BD}_{h}$):
$\mathcal{D}_{h}=\sup_{X,Y\in\mathcal{C}}\mathcal{BD}_{h}(X\|Y)$

Lipschitz Parameter (of $\mathcal{F}$): $\gamma_{\mathcal{F}}$

Proximal Map:
$\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{A},\mathcal{B}}(X,Y):={\arg\min}_{C\in\mathcal{C}}\\{\langle
Y,C-X\rangle-\frac{1}{\eta}\mathcal{BD}_{h}(C||X)\\}.$

$\eta\leftarrow\ {\mu_{h}}/(2\gamma_{\mathcal{F}})$ $\triangleright$ Step Size

$N\leftarrow\lceil\mathcal{D}_{h}/(\eta\cdot\epsilon)\rceil$ $\triangleright$
Number of Rounds

$(\alpha_{0},\beta_{0})\leftarrow\left(\frac{1}{2^{n}}\mathbb{1}_{\mathcal{A}},\frac{1}{2^{m}}\mathbb{1}_{\mathcal{B}}\right)$
$\triangleright$ State Initialization

$(\hat{\alpha}_{0},\hat{\beta}_{0})\leftarrow(\alpha_{0},\beta_{0})$
$\triangleright$ Intermediate State Initialization

for $t\in[1,N-1]$ do

$\alpha_{t+1}\leftarrow\Lambda(\log\hat{\alpha}_{t}+\eta\mathcal{F}_{\alpha}(\beta_{t}))$
$\triangleright$ State Updates

$\beta_{t+1}\leftarrow\Lambda(\log\hat{\beta}_{t}+\eta\mathcal{F}_{\beta}(\alpha_{t}))$

$\hat{\alpha}_{t+1}\leftarrow\Lambda(\log\hat{\alpha}_{t}+\eta\mathcal{F}_{\alpha}(\beta_{t+1}))$
$\triangleright$ Momentum Updates

$\hat{\beta}_{t+1}\leftarrow\Lambda(\log\hat{\beta}_{t}+\eta\mathcal{F}_{\beta}(\alpha_{t+1}))$

end for

return
$\bigg{(}\bar{\alpha}=\frac{1}{N}\sum_{t=0}^{N-1}\alpha_{t},\bar{\beta}=\frac{1}{N}\sum_{t=0}^{N-1}\beta_{t}\bigg{)}$

Our approach is based on endowing the problem’s feasible region— a product of
spectraplexes —with an OMD template in the spirit of Rakhlin and Sridharan
[53]. A difficulty that arises when trying to combine these elements is that
the product structure of standard exponential weight algorithms does not carry
over automatically to the non-commutative matrix variables that arise in the
quantum setting. The key observation that allows us to retain the speed-up of
a mirror prox method and the state-of-art dimensionality dependence is that
the von Neumann entropy (which is the underlying regularizer of the MMWU
algorithm) can be encoded in a matrix exponentiation step, as per the original
MMWU algorithm (see also [55] for generalizations o this connection for
optimization over symmetric cones). This allows us to implement an optimistic
update structure in the problem’s dual space, where the primary arithmetic
operation is ordinary addition—which is an abelian operation on the space of
Hermitian matrices. Thus, drawing inspiration from a classical hierarchy of
optimization methods, we propose a _design of quantum zero-sum game
algorithms_ (depicted in Figure 2) achieving an optimized convergence rate
with a single gradient call per iteration.

Within this design, we demonstrate that the MMWU algorithm of [29]
instantiates (via a von Neumann entropy regularizer) the _Matrix Dual
Averaging (MDA)_ method (Algorithm 4). Following a classical proof technique
in the spirit of Ene and Nguyên [56], we also prove that another of the
proposed algorithms in our hierarchy—the Optimistic Matrix Mirror Prox (OMMP)
method (Algorithm 5)—obtains a quadratic speedup relative to the convergence
rate of MDA, while still requiring only one gradient call per iteration:

######  Theorem 1 (Main Result)

The OMMP method computes $\epsilon$-Nash equilibria in
$\mathcal{O}_{d}({1}/{\epsilon})$ steps.

Thus, OMMP establishes a new benchmark for computing $\epsilon$-Nash finite-
valued quantum zero-sum games. It is by instantiating OMMP with a von Neumann
entropy regularizer, that the OMMWU algorithm and $\mathcal{O}(d/{\epsilon})$
iteration complexity is achieved.

Proof Insights. The core idea of an _optimistic_ algorithm is to leverage the
gradient from the previous iteration as a forecast for the subsequent
iteration’s gradient. This is predicated on the presumption that if the the
game’s adversary employs a consistent (or stabilized) strategy, then the
discrepancy between the gradients of successive iterations will be negligible.
Therefore, the formulation of our algorithm is crafted by devising a parallel
to such methods, circumventing the limitations imposed by the non-
commutativity of the matrix operator, which is a common hurdle in quantum
problems. Upon establishing the algorithm in this manner, we demonstrate that
the (quantum) superoperators, corresponding to the method’s gradients, exhibit
two significant characteristics: _monotonicity_ and _Lipschitz
continuity/smoothness_.

Monotonicity is a valuable property that allows us to: $(i)$ correlate any
Nash equilibrium with a more manageable variational inequality problem, and
$(ii)$, akin to convexity in single-agent optimization, it guarantees a
consistent progression towards the equilibrium.

Lipschitz continuity, in contrast, is the essential quality that accelerates
the convergence rate from $O(1/\epsilon^{2})$ to $O(1/\epsilon)$. To
comprehend the significance of this attribute, let’s view the situation
through the lens of the minimization participant—named Bob, with his
maximization counterpart being Alice. Bob finds himself involved in an online
convex minimization game. In the absence of any predictive information about
the forthcoming functions, which are influenced by Alice’s choices, folklore
results suggests that the optimal rate is $\Theta(1/\sqrt{T})$, or
equivalently $\Theta(1/\epsilon^{2})$ iterations. However, should Alice’s
strategy remain relatively static, Bob can anticipate that the upcoming
minimization challenges will be similar, affording him a critical advantage
for achieving a rate of $O(1/T)$, or $O(1/\epsilon)$.

This insight reveals that the algorithm’s performance is deeply tied to the
Lipschitz parameter of the gradient operator. Here, the concept of mirror
proxies plays a vital role in the proof. Utilizing mirror steps, distinct
regularizers enable either more assertive or more refined movements within the
strategy space of the spectraplex for both parties. Contrary to the
exponential dependency implied by a Frobenius norm, we show that the Von
Neumann entropy regularizer maintains a dimension-independent Lipschitz
constant under a norm that is congruent with the specific needs of our
problem.

Future Work. In addition to our novel methodology, we conclude with the
following conjecture that aims to tighten our method’s performance guarantees
and is likely of interest to both the quantum and optimization communities:

###### Conjecture 1

The rates of Theorem 1 are tight: specifically, there exists a quantum two-
player, zero-sum game for which OMMP methods require $\Omega_{d}(1/\epsilon)$
iterations to compute $\epsilon$-Nash equilibria.

###  1.4 Paper Organization

The paper is divided as follows. Section 2 reviews preliminaries, such as the
quantum zero-sum games setup and mathematical concepts necessary for our
proofs. Section 3 presents the series of quantum optimization methods for
computing approximate Nash equilibria of quantum zero-sum games, as depicted
in Figure 2. It begins with the MMWU method of [29] and shows how the
algorithm is an instantiation of the more general Matrix Dual Averaging (MDA)
method. It then introduces the Matrix Mirror Prox (MMP) method, which achieves
the optimal rate, but requires two gradients per iteration. The section
concludes with our proposed OMMP algorithm, which requires only a single
gradient call per iteration. Finally, Section 4 elaborates the convergence
analysis of OMMP algorithm, proving that it achieves an
$\mathcal{O}_{d}(1/\epsilon)$ rate. The paper concludes with a discussion of
different implementations of the OMMP method, via different regularizers. We
show that a von Neumann entropy regularizer, instantiating the OMMWU
algorithm, achieves a better dimensionality dependence than the Frobenius
regularizer. Thus, the OMMWU algorithm achieves the desired
$\mathcal{O}(d/\epsilon)$ and a quadratic speedup relative to MMWU. For
concision, we defer some details and proofs to the Appendix.

Figure 1: Design of Classical Zero-Sum Game Algorithms. This diagram provides
the update rules for and relationships between different learning algorithms
for classical zero-sum games. The left-hand side presents the Bregman
generalized methods, parameterized by distance-generating/regularization
function $h$. The right-hand side presents instantiations based on the
$\ell_{2}$ norm and entropy function, inducing an orthogonal projection
$\textrm{Orth}\Pi$ and logit map $\Lambda$, respectively. Note that moving
from $h=\ell_{2}$ to $h=$entropy results in a logarithmic improvement in the
total number of rounds $T$, by reducing the dependence on the simplex
dimension $d$. Furthermore, moving from the mirror map $\textnormal{Mir}\Pi$
of the Dual Averaging method (a “lazy” variant of the classic Mirror Descent
Ascent algorithm) to the proximal maps $\textnormal{Prox}\Pi$ of the Mirror
Prox method results in a quadratic improvement in convergence, achieving the
desired $O(1/\epsilon)$ rate. Furthermore, by reusing the “past gradient,” the
Single-Call Mirror Prox method retains the rate of Mirror Prox, but reduces
the total number of gradient calls per iteration from two to one.  Figure 2:
Design of Quantum Zero-Sum Game Algorithms. This diagram provides the update
rules for and relationships between the learning algorithms for quantum zero-
sum games as proposed in this work. The left-hand side presents the Bregman
generalized methods, parameterized by distance-generating/regularization
function $h$. The right-hand side presents instantiations based on the von
Neumann entropy function and Frobenius ($\ell_{2}$) norm, inducing an
orthogonal projection $\textrm{Orth}\Pi$ and logit map $\Lambda$,
respectively. Note that moving from $\ell_{2}$ to an entropy function as the
regularizer results in a logarithmic improvement in the total number of rounds
$T$, by reducing the dependence on the spectraplex dimension $D=4^{d}$.
Furthermore, moving from the mirror map $\textnormal{Mir}\Pi$ of the Matrix
Dual Averaging method (Jain and Watrous’ MMWU proposal [29]) to the proximal
maps $\textnormal{Prox}\Pi$ of the Matrix Mirror Prox method results in a
quadratic improvement in convergence, achieving the desired $O(1/\epsilon)$
rate. Furthermore, by reusing the “past gradient,” the Single-Call Matrix
Mirror Prox method retains the rate of Matrix Mirror Prox, but reduces the
total number of gradient calls per iteration from two to one.

##  2 Preliminaries

###  2.1 Quantum Zero-Sum Games

In this work, we consider the restricted class of non-cooperative zero-sum
quantum games. Within this setting, two competing players, Alice and Bob, each
transmit a mixed quantum state to a referee (henceforth referred to as Roger),
who performs a joint measurement on both states. Roger also possesses a
utility function that attributes payoffs to the players based on the
measurement outcomes. Since this is a zero-sum game, Bob’s payoff will be the
negative value of Alice’s—e.g., if Alice wins $5, Bob will lose $5.

Concretely, Alice and Bob, play the game by independently preparing multi-
qubit quantum states $\alpha$ and $\beta$, respectively, to be sent to Roger.
It is important to note that these are mixed states, corresponding to “mixed
strategies,” and that Alice and Bob do not share entanglement. Thus, Alice and
Bob’s set of possible moves is constrained to the set of density matrices
lying in their respective _spectraplexes_. We will assume Alice prepares
$n$-qubit states in her spectraplex,

$\displaystyle\mathcal{A}=\\{\alpha\in\mathcal{H}^{n}_{+}:\textnormal{Tr}(\alpha)=1\\},$
(1)

while Bob prepares $m$-qubit states in his spectraplex,

$\displaystyle\mathcal{B}=\\{\beta\in\mathcal{H}^{m}_{+}:\textnormal{Tr}(\beta)=1\\},$
(2)

where $\mathcal{H}^{d}_{+}$ denotes the set of $2^{d}\times
2^{d}$-dimensional, positive Hermitian matrices. In Section C.1, we show that
all spectraplexes are compact and convex sets. Denote the joint set of their
states as the direct sum

$\displaystyle\Psi=(\alpha,\beta)\in\mathcal{C}=\mathcal{A}\oplus\mathcal{B},$
(3)

with tensor product denoted by $\Psi^{\otimes}=\alpha\otimes\beta$. As a
general notational note, we will use lower-case variables (e.g. $\alpha$ or
$\beta$) to refer to states in the individual spectraplexes of either Alice
($\mathcal{A}$) or Bob ($\mathcal{B}$) and upper-case variables (e.g. $\Psi$
or $\Phi$) to denote joint states, which lie in the joint space
($\mathcal{C}$).

Alice and Bob send their states to Roger, who performs a joint measurement.
This joint measurement has a finite number of possible measurement outcomes,
$\Omega$, defining an $(n+m)$-qubit _positive operator-valued measurement
(POVM)_ ,

$\displaystyle\\{P_{\omega}\\}_{\omega\in\Omega},\text{ where
}\sum_{\omega\in\Omega}P_{\omega}=\mathbb{1}.$ (4)

Thus, Roger’s probability of observing outcome $\omega\in\Omega$ is

$p_{\omega}(\Psi)=\textnormal{Tr}\left(P_{\omega}^{\dagger}\Psi^{\otimes}\right)$
(5)

Dependent on the measurement outcome, Roger gives Alice and Bob a
reward/payoff. Let

$\mathcal{U}:\Omega\rightarrow\mathbb{R}$ (6)

denote Alice’s _utility_ function, which we will assume to be finite-valued
and have image $[-1,1]$ under normalization. We define Alice’s _payoff
observable_ as

$\displaystyle U=\sum_{\omega\in\Omega}\mathcal{U}(\omega)P_{\omega},$ (7)

such that Alice’s _expected payoff_ is

$u_{\text{Alice}}:=u(\Psi_{t})=\textnormal{Tr}(U^{\dagger}\Psi^{\otimes})=\sum_{\omega\in\Omega}\mathcal{U}(\omega)\textnormal{Tr}(P_{\omega}^{\dagger}\Psi^{\otimes})=\sum_{\omega\in\Omega}\mathcal{U}(\omega)p_{\omega}(\Psi).$
(8)

Since the game is zero-sum, this implies that Bob’s utility is $-\mathcal{U}$,
with payoff observable $-U$ and expected payoff $u_{\text{Bob}}:=-u(\Psi)$. In
playing the game, Alice and Bob aim to maximize their expected payoff, written
as $\max_{\alpha}u(\alpha,\beta)$ and $\max_{\beta}u(\alpha,\beta)$,
respectively. Since the game is zero-sum, in maximizing his expected payoff,
Bob can be seen as trying to minimize Alice’s payoff:
$\max_{\beta}-u(\alpha,\beta)=\min_{\beta}u(\alpha,\beta).$

Finally, since $u(\alpha,\beta)$ is a bilinear function and the spectraplexes
are compact and convex sets, von Neumann’s Min-Max Theorem [1] states

$\max_{\alpha}\min_{\beta}u(\alpha,\beta)=\min_{\beta}\max_{\alpha}u(\alpha,\beta).$
(9)

###  2.2 Quantum Nash Equilibria

A _Nash equilibrium_ of the quantum game is a pair of states
$(\alpha^{*},\beta^{*})$, such that each player has no incentive to change to
a different state unilaterally:

$\displaystyle u(\alpha^{*},\beta^{*})$ $\displaystyle\geq
u(\alpha,\beta^{*}),\>\>\forall\alpha\in\mathcal{A}$ (10) $\displaystyle
u(\alpha^{*},\beta^{*})$ $\displaystyle\leq
u(\alpha^{*},\beta),\>\>\forall\beta\in\mathcal{B}.$ (11)

Since $\mathcal{C}$ is convex and $u$ is linear in $\alpha$ and $\beta$, the
existence of Nash equilibria follows from Debreu’s equilibrium existence
theorem [57].

In playing the game, Alice and Bob receive feedback from Roger. In the
original work on quantum zero-sum games, Jain and Watrous [29] define Alice’s
feedback as a quantum channel, given by the superoperator
$\Xi:\mathcal{B}\rightarrow\mathcal{A}$, applied solely to Bob’s state. The
superoperator is parameterized by the payoff observable as

$\Xi(\beta)=\textnormal{Tr}_{\mathcal{B}}[U(\mathbb{1}_{\mathcal{A}}\otimes\beta^{\top})],$
(12)

such that the expected utility can be expressed as

$u(\alpha,\beta)=\textnormal{Tr}[U(\alpha\otimes\beta)]=\textnormal{Tr}\left[\alpha\>\textnormal{Tr}_{\mathcal{B}}[U(\mathbb{1}_{\mathcal{A}}\otimes\beta)]\right]=\textnormal{Tr}[\alpha\>\Xi(\beta^{\top})].$
(13)

Similarly, Bob’s feedback is defined by the quantum channel, corresponding to
adjoint superoperator $\Xi^{*}:\mathcal{A}\rightarrow\mathcal{B}$, applied
solely to Alice’s state. This adjoint superoperator is defined as

$\Xi^{*}(\alpha)=\textnormal{Tr}_{\mathcal{A}}[U(\alpha^{\top}\otimes\mathbb{1}_{\mathcal{B}})].$
(14)

and is uniquely determined by the condition

$u(\alpha,\beta)=\textnormal{Tr}[\alpha\>\Xi(\beta^{\top})]=\textnormal{Tr}[\Xi^{*}(\alpha^{\top})\>\beta].$
(15)

In this work, we will move away from thinking about Roger’s feedback in terms
of quantum channels towards a gradient-based view, as leveraged in recent work
on learning general quantum games [54]. It is mathematically equivalent to
think of the channel-based feedback in terms of gradients of the expected
utility with respect to Alice and Bob’s states. For a review of complex
differentiation and gradients, refer to Appendix A. Namely, we define Alice’s
feedback, or _payoff gradient_ , as

$\displaystyle\mathcal{F}_{\alpha}(\beta)$
$\displaystyle=\nabla_{\alpha^{\top}}u(\alpha,\beta)=\nabla_{\alpha^{\top}}\textnormal{Tr}\left[U^{\dagger}(\alpha\otimes\beta)\right]=\textnormal{Tr}_{\mathcal{B}}\left[U^{\dagger}(\mathbb{1}_{n}\otimes\beta)\right],$
(16)

which is equivalent to $\Xi(\beta^{\top})$ of Equation 12. Similarly Bob’s
feedback is defined as the payoff gradient

$\displaystyle\mathcal{F}_{\beta}(\alpha)$
$\displaystyle=-\nabla_{\beta^{\top}}u(\alpha,\beta)=-\nabla_{\beta^{\top}}\textnormal{Tr}\left[U^{\dagger}(\alpha\otimes\beta)\right]=-\textnormal{Tr}_{\mathcal{A}}\left[U^{\dagger}(\alpha\otimes\mathbb{1}_{m})\right],$
(17)

which is equivalent to $-\Xi^{*}(\alpha^{\top})$ of Equation 14. Similarly to
Equation 15, we have that

$\displaystyle
u(\alpha,\beta)=\textnormal{Tr}\left[{\alpha}\>{\mathcal{F}_{\alpha}(\beta)}\right]=-\textnormal{Tr}\left[{\beta}\>{\mathcal{F}_{\beta}(\alpha)}\right].$
(18)

Under the gradient-based view,
$\textnormal{Tr}[\alpha\>\mathcal{F}_{\alpha}(\beta)]=\textnormal{Tr}[\alpha\>\nabla_{\alpha^{\top}}u(\alpha,\beta)]$
can be interpreted as the directional derivative of $u$, as a function of
Alice’s state, in the direction of Alice’s state. Similarly,
$-\textnormal{Tr}[\mathcal{F}_{\beta}(\alpha)\>\beta]=\textnormal{Tr}[\nabla_{\beta^{\top}}u(\alpha,\beta)\>\beta]$
can be interpreted as the directional derivative of $u$, as a function of
Bob’s state, in the direction of Bob’s state. As will be demonstrated
throughout the work, moving away from a channel-based to this gradient-based
view of the quantum zero-sum game feedback proves invaluable for mapping
gradient-based insights from the classical games and optimization literature
directly to the quantum setting.

Henceforth, we will denote the joint state as $\Psi=(\alpha,\beta)$, the joint
feedback of both players, as follows:222In a slight abuse of notation, let
$\mathcal{F}(\Psi)=\mathcal{F}(\alpha,\beta)$.

$\displaystyle\mathcal{F}(\Psi)=\Big{(}\mathcal{F}_{\alpha}(\beta),\mathcal{F}_{\beta}(\alpha)\Big{)},$
(19)

and the game’s Nash equilibrium as $\Psi^{*}=(\alpha^{*},\beta^{*})$. Thus, as
depicted in Figure 3, in each round $t$ of the game, after playing the joint
state $\Psi_{t}$, Alice will receive the individual payoff gradient
$\mathcal{F}_{\alpha_{t}}(\beta_{t})$ and Bob will receive
$\mathcal{F}_{\beta_{t}}(\alpha_{t})$, comprising joint feedback
$\mathcal{F}(\Psi_{t})$.

Figure 3: An illustration of a single round of the quantum zero-sum game
(gradient-based view).

By using standard arguments [58], it can be shown that the solutions of the
variational inequality

$\displaystyle\textnormal{Tr}\left[{(\Psi-\Psi^{*})}\>{\mathcal{F}(\Psi^{*})}\right]\leq
0,\>\>\>\>\forall\>\Psi\in\mathcal{C},$ (20)

are the game’s Nash equilibria, $\Psi^{*}$.333Note that in the direct sum
notation, for matrices $A_{1}$, $A_{2}$, $B_{1}$, and $B_{2}$, where
$\text{dim}(A_{1})=\text{dim}(A_{2})=n\times n$ and
$\text{dim}(B_{1})=\text{dim}(B_{2})=m\times m$, we have that
$\textnormal{Tr}[(A_{1},B_{1})(A_{2},B_{2})]=\textnormal{Tr}[A_{1}A_{2}]+\textnormal{Tr}[B_{1}B_{2}]$.
Specifically, solutions satisfying Equation 20 are known as strong solutions
to the variational inequality. Alternatively, weak solutions to the
variational inequality are $\Psi^{*}\in\mathcal{C}$ satisfying

$\textnormal{Tr}\left[{(\Psi^{*}-\Psi)}\>{\mathcal{F}(\Psi)}\right]\geq
0,\>\>\>\>\forall\>\Psi\in\mathcal{C}.$ (21)

We offer intuitive visualizations of both the weak and strong solutions of
variational inequalities as vector inner products in Figure 4.

Figure 4: Intuitive illustrations of the strong and weak solutions $\Psi^{*}$
of the variational inequalities given in Equation 20 and Equation 21,
respectively. For strong solutions, the inner product between $\Psi-\Psi^{*}$
and $\mathcal{F}(\Psi^{*})$ must be non-negative for all $\Psi$ in the joint
spectraplex $\mathcal{C}$. Meanwhile, for weak solutions, the inner product
between $\Psi^{*}-\Psi$ and $\mathcal{F}(\Psi)$ must be non-positive for all
$\Psi\in\mathcal{C}$.

As we will discuss in Section 4 and prove in Section C.2, $\mathcal{F}$ is
both monotone (Lemma 8) and Lipschitz-continuous (Lemma 9). Additionally, in
Appendix C, we prove the following two important lemmas:

###### Lemma 1

If $\mathcal{F}$ is monotone, then a strong solution is a weak solution.

###### Lemma 2

If $\mathcal{F}$ is Lipschitz continuous, then a weak solution is a strong
solution.

Therefore, by the monotonicity and Lipschitzness of $\mathcal{F}$, the strong
and weak solutions (described, respectively, by Equation 20 and Equation 21)
are equally valid criterion for Nash equilibria of the finite-valued zero-sum
quantum games under consideration. We will generally leverage the weak-
solution criterion in our analysis.

Finally, in this work, we focus on the iteration complexity of our methods for
computing an $\epsilon$-approximate Nash equilibrium,
$\tilde{\Psi}^{*}=(\tilde{\alpha}^{*},\tilde{\beta}^{*})$, which—as defined by
[29]—must satisfy

$\max_{\alpha}u(\alpha,\tilde{\beta}^{*})-\epsilon||U||_{\infty}\leq
u(\tilde{\alpha}^{*},\tilde{\beta}^{*})\leq\min_{\beta}u(\tilde{\alpha}^{*},\beta)+\epsilon||U||_{\infty},$

where $||U||_{\infty}$ denotes the infinity norm of the payoff observable $U$.
From this definition and the maximal utility of Nash equilibria it follows
that

$\displaystyle\epsilon||U||_{\infty}\geq
u(\tilde{\alpha}^{*},\tilde{\beta}^{*})-\min_{\beta}u(\tilde{\alpha}^{*},\beta)$
$\displaystyle\geq 0$ (22)
$\displaystyle\epsilon||U||_{\infty}\geq\max_{\alpha}u(\alpha,\tilde{\beta}^{*})-u(\tilde{\alpha}^{*},\tilde{\beta}^{*})$
$\displaystyle\geq 0.$ (23)

Summing these equations gives the duality gap,

$\mathcal{DG}(\alpha^{\prime},\beta^{\prime})=\max_{\alpha}u(\alpha,\beta^{\prime})-\min_{\beta}u(\alpha^{\prime},\beta),$
(24)

which effectively serves as a measure of proximity between any proposed
states, $(\alpha^{\prime},\beta^{\prime})$, and an exact Nash equilibrium,
$(\alpha^{*},\beta^{*})$. From Equation 22 it follows that for an
$\epsilon$-Nash equilibrium, the duality gap is bounded as

$\epsilon^{\prime}=2\epsilon||U||_{\infty}\geq\mathcal{DG}(\alpha^{\prime},\beta^{\prime})\geq
0.$ (25)

Therefore, to evaluate the error of a proposed solution $\tilde{\Psi}$,
relative to any true Nash equilbrium, we leverage the notion of an
_error/merit function_ , which corresponds with the duality gap in the setting
of min-max optimization:

$\displaystyle\textnormal{Error}(\tilde{\Psi})=\sup_{\Psi\in\mathcal{C}}\textnormal{Tr}\left[{(\Psi-\tilde{\Psi})}\>{\mathcal{F}(\Psi)}\right].$
(26)

It is easy to check that $\textnormal{Error}(\tilde{\Psi})\geq 0$ for all
$\tilde{\Psi}$ and $\textnormal{Error}(\tilde{\Psi})=0$ if and only if
$\tilde{\Psi}=\Psi^{*}$.

###  2.3 Convexity, Smoothness, & Duality

Crucial to our analysis will be functional notions of convexity, smoothness,
and duality [59, 60]. We begin by defining convexity and strong convexity of a
differentiable function $f:\mathcal{C}\rightarrow\mathbb{R}$.

######  Definition 1 (Convex Function)

A differentiable function $f:\mathcal{C}\rightarrow\mathbb{R}$ is _convex_ if

$\displaystyle f(X)\geq f(Y)+\langle\nabla f(Y),X-Y\rangle,\hskip
7.22743pt\forall X,Y\in\mathcal{C}.$ (27)

######  Definition 2 ($\mu$-Strongly Convex Function)

A convex differentiable function $f:\mathcal{C}\rightarrow\mathbb{R}$ is
_$\mu$ -strongly convex_ with respect to norm $\|\cdot\|$ if

$\displaystyle f(X)\geq f(Y)+\langle\nabla
f(Y),X-Y\rangle+\frac{\mu}{2}\|X-Y\|^{2},\hskip 7.22743pt\forall
X,Y\in\mathcal{C},$ (28)

or equivalently, following from the definition of the Bregman divergence in
Equation 40,

$\displaystyle\mathcal{BD}_{f}(X\|Y)\geq\frac{\mu}{2}\|X-Y\|^{2},\hskip
7.22743pt\forall X,Y\in\mathcal{C}.$ (29)

Since strong convexity implies standard convexity (with $\mu=0$), the optimal
points (minima/maxima) of any strongly convex function must satisfy first-
order optimality.

######  Fact 1 (First-Order Optimality Condition)

For an optimization problem,

$\displaystyle\min f(X)\>\>\text{ subject to }\>\>X\in\mathcal{X},$

with convex differentiable $f$, a feasible point $X$ is optimal if and only if

$\displaystyle\langle\nabla f(X),Y-X\rangle\geq 0,\>\>\forall
Y\in\mathcal{X}.$

Another property of interest for the convex differentiable function $f$ is
smoothness, which implies Lipschitz continuity of $\nabla f$.

######  Definition 3 ($\beta$-Smooth Function)

A convex differentiable function $f:\mathcal{C}\rightarrow\mathbb{R}$ is
$\beta$-smooth with respect to norm $\|\cdot\|$ if

$\displaystyle f(X)\leq f(Y)+\langle\nabla
f(Y),X-Y\rangle+\frac{\beta}{2}\|X-Y\|^{2},\hskip 7.22743pt\forall
X,Y\in\mathcal{C}.$ (30)

or equivalently, following from the definition of the Bregman divergence in
Equation 40,

$\displaystyle\mathcal{BD}_{f}(X\|Y)\leq\frac{\beta}{2}\|X-Y\|^{2},\hskip
7.22743pt\forall X,Y\in\mathcal{C}.$ (31)

Before establishing the Lipschitz condition of $\nabla f$, we introduce the
_dual norm_.

######  Definition 4 (Dual Norm)

If the generic norm of a matrix $X$ is denoted $\|X\|$, the dual norm is
defined as

$\displaystyle\|X\|_{*}=\sup\\{\langle Y,X\rangle:\|Y\|\leq 1\\}.$ (32)

As shown in [61, 62], the Frobenius norm $\|\cdot\|_{F}$ is self-dual, meaning
$\|X\|_{*}=\|X\|_{F}$.

######  Proposition 1 ($\beta$-Lipschitz)

A convex differentiable function $f:\mathcal{C}\rightarrow\mathbb{R}$ is
$\beta$-smooth with respect to norm $\|\cdot\|$ if its gradient $\nabla f$ is
$\beta$-Lipschitz continuous,

$\displaystyle\|\nabla f(X)-\nabla f(Y)\|_{*}\leq\beta\|X-Y\|\hskip
7.22743pt\forall X,Y\in\mathcal{C}.$ (33)

For the equivalence of Proposition 1 and Definition 3, see [63].
Geometrically, $\beta$-smoothness can be interpreted as a convex quadratic
_upper_ bound,

$\displaystyle g_{X}^{+}(Y)=f(Y)+\langle\nabla
f(Y),X-Y\rangle+\frac{\beta}{2}\|X-Y\|^{2},$ (34)

of $f$ such that $g_{X}^{+}(X)=f(X)$ and $g_{X}^{+}(Y)\geq f(Y),\>\forall
Y\in\mathcal{C}$. Meanwhile, $\mu$-strong convexity can be interpreted as a
convex quadratic _lower_ bound,

$\displaystyle g_{X}^{-}(Y)=f(Y)+\langle\nabla
f(Y),X-Y\rangle+\frac{\mu}{2}\|X-Y\|^{2},$ (35)

of $f$ such that $g_{X}^{-}(X)=f(X)$ and $g_{X}^{-}(Y)\leq f(Y),\>\forall
Y\in\mathcal{C}$. In this way, strong convexity is dual to smoothness.

Fenchel conjugacy relates two equivalent representations of a convex function.
In the standard representation, a convex function $f$ is represented by pairs
$(X,f(X))$ consisting of coordinates $X$ and corresponding function
evaluations $f(X)$. Alternatively, as depicted in Figure 5, a convex function
can be represented by its tangent at each point. These tangents are encoded as
the pairs $(\theta,f^{*}(\theta))$, where $\theta=\nabla f(X)$ for some $X$ is
the gradient and $-f^{*}(\theta)$ is the intercept. The function $f^{*}$,
which calculates the tangent’s intercept $f^{*}(\theta)$ from a given gradient
$\theta$, is known as the Fenchel conjugate function, formally defined as
follows.

######  Definition 5 (Fenchel Conjugate)

The _Fenchel conjugate_ of a function $f:\mathcal{C}\rightarrow\mathbb{R}$ is
the function $f^{*}:\mathcal{C}\rightarrow\mathbb{R}$ such that

$\displaystyle f^{*}(\theta)=\max_{X\in\mathcal{C}}\left\\{\langle
X,\theta\rangle-f(X)\right\\}.$ (36)

Relating back to Definition 4 of the dual norm, if it were the case that for
all $X$, $f(X)=0$, then $f^{*}(\theta)=\max_{X\in\mathcal{C}}\langle
X,\theta\rangle=||\theta||_{*}$. Therefore, the Fenchel conjugate can be seen
as maximizing the dual norm while also minimizing the penalty induced by
function $f$. In this way, $f$ can be seen as a “regularizer” function.

The definition of Fenchel conjugacy immediately implies the _Fenchel-Young
Inequality_ ,

$\displaystyle\forall X\in\mathcal{C},\hskip 7.22743ptf^{*}(\theta)\geq\langle
X,\theta\rangle-f(X).$ (37)

Furthermore, taking the gradient of the Fenchel conjugate should result in the
value of $X$ which maximizes $\langle X,\theta\rangle-f(X)$, resulting in
_Danskin’s Theorem_ , as follows. For more details, refer to Proposition 4.5.1
of [64].

######  Proposition 2 (Danskin’s Theorem)

Let $f:\mathcal{C}\rightarrow\mathbb{R}$ be a strongly convex function. For
all $X\in\mathcal{C}$, we have

$\displaystyle\nabla
f^{*}(\theta)=\arg\min_{X\in\mathcal{C}}\left\\{f(X)-\langle
X,\theta\rangle\right\\}.$

Finally, as stated below and explained further in Corollary 3.5.11 of [65],
there is a duality between the the convexity of $f$ and the smoothness of
$f^{*}$.

######  Lemma 3 (Strong-Smooth Duality)

A closed convex function $f:\mathcal{C}\rightarrow\mathbb{R}$ is
$\alpha$-strongly convex with respect to norm $\|\cdot\|$ if and only if
$f^{*}$ is $\frac{1}{\alpha}$-smooth with respect to the corresponding dual
norm $\|\cdot\|_{*}$.

Figure 5: Illustration of a function $f$ and its Fenchel conjugate $f^{*}$.
The function $f$ is represented by a collection of pairs $(X,f(X))$, while its
Fenchel conjugate $f^{*}$ is visualized as the set of tangents to $f$,
characterized by pairs $(\theta,f^{*}(\theta))$. [Adapted from Figure 2.3 in
[60].]

###  2.4 Bregman Divergence & Mirror Map

Critical to understanding our proposed algorithm and its improved convergence
over Jain and Watrous’ MMWU proposal are the notions of “mirror maps” and
Bregman divergence.

Historically, the first proposals for constrained Gradient Descent Ascent
leveraged Euclidean ($\ell_{2}$ norm) projections in order to ensure that the
algorithm output would lie within the feasible set. However, it was eventually
realized that for certain problems, as we will discuss in the following
subsection, regularization based on this $\ell_{2}$ norm results in suboptimal
constants in the convergence rates and, by leveraging different distance
metrics, faster rates can be achieved. In order to develop algorithms which
could encompass any such desired distance metric, the Bregman divergence was
introduced.

Figure 6: Bregman Divergence. The Bregman divergence $\mathcal{BD}_{h}(U||V)$
between two states in the primal space $U,V\in\mathcal{X}^{d}$, relative to
regularizer $h$, is given by the difference in the dual space
($\mathbb{R}^{2^{d}\times 2^{d}}$) between $h(U)$ and its linear approximation
via the first-order Taylor expansion of $h$ around $V$, $\tilde{h}_{V}(U)$.

The Bregman divergence is parameterized by a “distance-generating” or
“regularizer” function. Formally, a _distance-generating function_ (DGF) on
the $d$-dimensional spectraplex $\mathcal{H}_{+}^{d}$ is defined as a convex,
lower semi-continuous function
$h:\mathbb{C}^{d}\rightarrow\mathbb{R}\cup\\{\infty\\}$, where the _effective_
domain of $h$ is $\mathcal{X}^{d}$ and $\mathbb{C}^{d}$ is the space of
$2^{d}\times 2^{d}$ complex matrices. Denoting the space of $2^{d}\times
2^{d}$ real matrices as $\mathbb{R}^{d}$, let $\partial
h:\mathcal{H}_{+}^{d}\rightarrow\mathbb{R}^{d}$ denote the _subdifferential_
of $h$ defined as444The subdifferential is the set-valued generalization of
the derivative for convex functions which are not necessarily differentiable.

$\partial h(X):=\left\\{V\in\mathbb{R}^{d}:h(Y)\geq
h(X)+\textnormal{Tr}[(Y-X)^{\dagger}V],\>\>\forall\>Y\in\mathcal{H}_{+}^{d}\right\\}.$
(38)

The _domain of subdifferentiability_ of $h$ is defined as

$\textrm{dom}\partial h:=\\{X\in\mathcal{X}^{d}:\partial
h(X)\neq\emptyset\\}.$ (39)

Having defined the DGF $h$ and its subdifferential $\partial h$, we are ready
to define the Bregman divergence. Given two states
$U,V\in\mathcal{H}_{+}^{d}$, the Bregman divergence induced by $h$ is the non-
negative function $\mathcal{BD}_{h}:\textrm{dom}h\times\textrm{dom}\partial
h\rightarrow\mathbb{R}$ defined as

$\mathcal{BD}_{h}(U||V):=h(U)-h(V)-\textnormal{Tr}[(U-V)^{\dagger}\partial
h(V)],$ (40)

where the subdifferential maps from states in the _feasible region_
$\mathcal{H}_{+}^{d}$ of the _primal space_ $\mathbb{C}^{d}$, to the _dual
space_ of subgradients $\mathbb{R}^{d}$. The function $\partial h$, which
performs a primal-to-dual mapping, is often referred to as a _mirror map_.

Intuitively, as illustrated in Figure 6, the Bregman divergence
$\mathcal{BD}_{h}(U||V)$ between two states $U$ and $V$ in the primal space,
relative to regularizer $h$, is given by the difference in the dual space
between $h(U)$ and its linear approximation via the first-order Taylor
expansion of $h$ around $V$,

$\widetilde{h}_{V}(U)=h(V)+\textnormal{Tr}[(U-V)^{\dagger}\partial h(V)].$
(41)

Since $h$ is convex, this difference is always non-negative. Note that
$\mathcal{BD}_{h}(U||V)$ is not necessarily equivalent to
$\mathcal{BD}_{h}(V||U)$.

Finally, our analysis will use the well-known _Three-Point Identity_ ,

$\textnormal{Tr}[(Z-X)^{\dagger}\nabla_{X^{\top}}\mathcal{BD}_{h}(X||Y)]=\mathcal{BD}_{h}(Z||Y)-\mathcal{BD}_{h}(X||Y)-\mathcal{BD}_{h}(Z||X),$
(42)

which follows from the definition of the Bregman divergence [66].

###  2.5 Mirror & Proximal Steps

We will now introduce “mirror steps” and “proximal steps,” which incorporate
the mirror map generalization of distance in constrained optimization
updates.555For further background, we refer the interested reader to [66] and
Chapter 4 of [67]. These steps implement the algorithmic state updates of the
classical games algorithms discussed in the Prior Work (Section 1.2). Similar
to a standard constrained optimization gradient step, the goal of these
procedures is to produce an update that lies in the feasible set, closer to
the optimum. However, by generalizing the Euclidean distance of traditional
Gradient Descent Ascent steps, mirror and proximal steps enable incorporation
of domain-specific knowledge about the geometry of the optimization task into
the updates. This can lead to substantial improvements in convergence rate.

For intuition, let us first consider a classical single-agent case as an
illustrative example. When optimizing a function $f$ over a
$(d-1)$-dimensional simplex, restricted by the $\ell_{\infty}$-norm, standard
Projected Subgradient Descent (PSD) reaches a convergence rate of
$\mathcal{O}(\sqrt{d/T})$. In contrast, classical Mirror Descent Ascent using
a negative Shannon entropy regularizer significantly improves the rate to
$\mathcal{O}(\sqrt{\log d/T})$.

Our objective is to adapt this methodology to the more challenging,
exponentially larger $(4^{d}-1)$-dimensional spectraplex
$\mathcal{H}^{d}_{+}$, in a min-max framework. Specifically, we aim to
incorporate domain-specific regularization techniques to achieve an
$\mathcal{O}(d/T)$ convergence rate. This would offer an exponential
dimensionality improvement from standard Projected Subgradient Ascent Descent,
which has an $\mathcal{O}(2^{d}/\sqrt{T})$ convergence rate [68, 69], and a
quadratic speed-up over the Jain-Watrous algorithm, which has an
$\mathcal{O}(d/\sqrt{T})$ convergence rate [29].

Note that intuitions for the following definitions and proofs are provided in
Appendix B.

Figure 7: Mirror & Proximal Step.Visualizations of the mirror and proximal
steps described in Algorithm 2 and Algorithm 3, respectively. [Adapted from
Figures 4.1 and 4.2 in [67].]

####  2.5.1 Mirror Step

Given a function $f:\mathbb{C}^{d}\rightarrow\mathbb{R}$ that we aim to
optimize, convex feasible set $\mathcal{H}_{+}^{d}\subset\mathbb{C}^{d}$, DGF
$h$, subgradients $\partial h$, and initial primal state
$X_{t}\in\mathcal{H}_{+}^{d}$, each call to a mirror step oracle generates an
updated state $X_{t+1}$ via Algorithm 2 (as depicted in Figure 7). In Step
(3), the mirror map, $\partial h^{-1}$, can be implemented as follows.

######  Definition 6 (Mirror Map)

For a state $D$ in the dual space $\mathcal{D}$, the _mirror map_ , or
_regularized best response_ ,

$\displaystyle\textnormal{Mir}\Pi_{\mathcal{P}}^{h}(D):=\underset{P\in\mathcal{P}}{\arg\max}\\{\textnormal{Tr}[D^{\dagger}P]-h(P)\\}$
(43)

implements the inverse mirror map, $\partial
h^{-1}:\mathcal{D}\rightarrow\mathcal{P}$, mapping $D\in\mathcal{D}$ to the
corresponding state $P$ in the feasible region $\mathcal{P}$ of the primal
space.

Algorithm 2 Mirror Step

  1. 1.

Initial primal state $X_{t}\in\mathcal{H}_{+}^{d}$ is mapped to a dual state
$\nabla h(X_{t})\in\partial h(X_{t})\subset\mathbb{R}^{d}$.

  2. 2.

Gradient $\nabla f:\mathbb{C}^{d}\rightarrow\mathbb{R}^{d}$ is calculated with
respect to $X_{t}$.

  3. 3.

In the dual space, the state is updated according to the gradient with step
size $\eta$, resulting in a new dual state

$\nabla h(X_{t+1})=\nabla h(X_{t})+\eta\nabla f(X_{t})\>\>\in\mathbb{R}^{d}.$

  4. 4.

Mapping back to the primal space, this corresponds to the updated feasible
state

$X_{t+1}\in\partial h^{-1}(\nabla h(X_{t+1})+\eta\nabla
f(X_{t}))\>\>\in\mathcal{H}_{+}^{d}.$

Algorithm 3 Proximal Step

  1. 1.

From initial feasible primal state $X_{t}\in\mathcal{H}_{+}^{d}$, updated
feasible primal state $Y_{t+1}\in\mathcal{H}_{+}^{d}$ is generated via the
Mirror Step procedure of Algorithm 2.

  2. 2.

Initial feasible primal state $X_{t}$ is mapped to a dual state $\nabla
h(X_{t})\in\partial h(X_{t})\subset\mathbb{R}^{d}$.

  3. 3.

Gradient $\nabla f:\mathbb{C}^{d}\rightarrow\mathbb{R}^{d}$ is calculated with
respect to $Y_{t+1}$.

  4. 4.

In the dual space, the state is updated according to the gradient with step
size $\eta$, resulting in a new dual state

$\nabla h(X_{t+1})=\nabla h(X_{t})+\eta\nabla
f(Y_{t+1})\>\>\in\mathbb{R}^{d}.$

  5. 5.

Mapping back to the primal space, this corresponds to the updated feasible
state

$X_{t+1}\in\partial h^{-1}(\nabla h(X_{t})+\eta\nabla
f(Y_{t+1}))\>\>\in\mathcal{H}_{+}^{d}.$

Therefore, the full mirror step procedure, mapping from $X_{t}$ to $X_{t+1}$,
is achieved via the mirror map of $Y_{t}=\nabla h(X_{t})+\eta\nabla f(X_{t})$
on the feasible set $\mathcal{H}_{+}^{d}$, as

$X_{t+1}=\textnormal{Mir}\Pi_{\mathcal{H}_{+}^{d}}^{h}(\nabla
h(X_{t})+\eta\nabla f(X_{t})).$ (44)

Alternatively, the mirror step procedure can be expressed as a _mirror
proximal map_.

######  Lemma 4 (Mirror Step)

For state $X_{t}$ in the feasible set $\mathcal{H}_{+}^{d}$ of primal space
$\mathbb{C}^{d}$,

$\displaystyle\textnormal{MirProx}\Pi_{\mathcal{X}^{d}}^{h,\eta,f}(X_{t}):=\underset{X_{t+1}\in\mathcal{H}_{+}^{d}}{\arg\min}\>\left\\{\textnormal{Tr}[\nabla
f(X_{t})^{\dagger}X_{t+1}]+\frac{1}{\eta}\mathcal{BD}_{h}(X_{t+1}||X_{t})\right\\},$
(45)

implements the mirror step mapping to an updated feasbile set primal state
$X_{t+1}\in\mathcal{H}_{+}^{d}$.

####  2.5.2 Proximal Step

We refer to the mirror prox projection as a _proximal map_. Its generalization
comprises the _proximal step_ oracle, described in Algorithm 3 and depicted in
Figure 7.

######  Definition 7 (Proximal Map)

For a state $P$ in the feasible set of the primal space $\mathcal{P}$ and a
state $D$ in the dual space, the _proximal map_ ,

$\displaystyle\textnormal{Prox}\Pi^{h}_{\mathcal{P}}(P,D):$
$\displaystyle=\underset{P^{*}\in\mathcal{P}}{\arg\min}\left\\{\textnormal{Tr}[D^{\dagger}(P^{*}-P)]+\mathcal{BD}_{h}(P^{*}||P)\right\\}$
$\displaystyle=\underset{P^{*}\in\mathcal{P}}{\arg\min}\left\\{\textnormal{Tr}[D^{\dagger}P^{*}]+\mathcal{BD}_{h}(P^{*}||P)\right\\},$
(46)

implements the proximal mapping to an updated feasible set primal state
$P^{*}\in\mathcal{P}$.

Applying this definition to our problem of interest yields the following map.

######  Lemma 5 (Proximal Step)

For state $X_{t}\in\mathcal{H}_{+}^{d}$ in the feasible set of the primal
space,

$\displaystyle\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{H}_{+}^{d}}(X_{t},\nabla
f(Y_{t+1})):$
$\displaystyle=\underset{X_{t+1}\in\mathcal{H}_{+}^{d}}{\arg\min}\left\\{\textnormal{Tr}[\nabla
f(Y_{t+1})^{\dagger}(X_{t+1}-X_{t})]+\frac{1}{\eta}\mathcal{BD}_{h}(X_{t+1}||X_{t})\right\\}$
$\displaystyle=\underset{X_{t+1}\in\mathcal{H}_{+}^{d}}{\arg\min}\left\\{\textnormal{Tr}[\nabla
f(Y_{t+1})^{\dagger}X_{t+1}]+\frac{1}{\eta}\mathcal{BD}_{h}(X_{t+1}||X_{t})\right\\}$
(47)

implements the proximal step mapping to an updated feasbile set primal state
$X_{t+1}\in\mathcal{H}_{+}^{d}$.

Note that $\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{H}_{+}^{d}}(X_{t},\nabla
f(X_{t}))=\textnormal{MirProx}\Pi_{\mathcal{H}_{+}^{d}}^{h,\eta,f}(X_{t})$,
which, by Lemma 4, implies that the proximal step is a generalization of the
mirror step and, thus, can be used to implement a mirror step. Thus, the
proximal step procedure of Algorithm 3 can be achieved algorithmically via two
proximal maps as:

$\displaystyle Y_{t+1}$
$\displaystyle=\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{H}_{+}^{d}}(X_{t},\nabla
f(X_{t}))$ (48) $\displaystyle X_{t+1}$
$\displaystyle=\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{H}_{+}^{d}}(X_{t},\nabla
f(Y_{t+1})).$ (49)

Note that the proximal step can be alternatively written as a single update:

$\displaystyle X_{t+1}$
$\displaystyle=\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{H}_{+}^{d}}\bigg{(}X_{t},\nabla
f\Big{(}\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{H}_{+}^{d}}(X_{t},\nabla
f(X_{t}))\Big{)}\bigg{)}.$ (50)

Standard gradient-based optimization methods are generally based on a single
step in the direction of the gradient. In contrast, proximal step methods are
based on a two-step approach:

  1. 1.

From initial point $X_{t}$, apply a “look-ahead” step, resulting in a
“prediction” $Y_{t+1}$.

  2. 2.

Use prediction $Y_{t+1}$ in a “correction” step that achieves the update
$X_{t+1}$.

As we will discuss in Section 3.3, replacing mirror with proximal steps
provides the desired acceleration from $\mathcal{O}_{d}(1/\epsilon^{2})$ to
$\mathcal{O}_{d}(1/\epsilon)$. Furthermore, the usage of the Bregman
divergence enables a geometry-aware projection and, specifically, can be used
to achieve a logarithmic improvement on dependence of the dimensionality of
the spectraplex, as will be discussed in Section 4.7.

####  2.5.3 Regularizers

In order to guarantee the well-posedness of the mirror and proximal
projections, we will restrict $h$ to the class of functions which are strongly
convex. In the classical games setting, as discussed in the Prior Work
(Section 1.2), the most widely used regularizers are the Euclidean
($\ell_{2}$) norm ($h(x)=\frac{1}{2}||x||_{2}^{2}$) and the negative entropy
($h(x)=\sum_{i}x_{i}\log(x_{i})$). In this work, we will leverage the quantum
analogs of these classical regularizers: the Frobenius norm
($h(X)=\frac{1}{2}||X||_{F}^{2}$) and the von Neumann entropy
($h(X)=\textnormal{Tr}[X\log X]$). As highlighted in the following lemmas
(proved in Appendix B), the mirror and proximal maps take well-known forms
when these regularizers are used.

######  Lemma 6 (Bregman Divergence of the Frobenius Norm)

If $h$ is the squared Frobenius norm,

$h({X})=\frac{1}{2}||X||_{F}^{2}=\frac{1}{2}\textnormal{Tr}[X^{\dagger}X],$
(51)

then the Bregman divergence is the squared Frobenius distance,

$\mathcal{BD}_{h}(X||Y)=\tfrac{1}{2}||X-Y||_{F}^{2},$ (52)

the corresponding mirror map (regularized best-response) is the orthogonal
projection,

$\displaystyle\textnormal{Mir}\Pi^{h}_{\mathcal{C}}(Y)={\arg\min}_{C\in\mathcal{C}}||Y-C||_{F}^{2}=\textrm{Orth}\Pi_{\mathcal{C}}(Y),$
(53)

and the corresponding proximal map is

$\displaystyle\textnormal{Prox}\Pi^{h}_{\mathcal{C}}(X,Y)={\arg\min}_{C\in\mathcal{C}}||X+Y-C||_{F}^{2}=\textrm{Orth}\Pi_{\mathcal{C}}(X+Y).$
(54)

######  Lemma 7 (Bregman Divergence of the von Neumann Entropy)

If $h$ is the von Neumann entropy,

$h({X})=\textnormal{Tr}[X\log X],$ (55)

then the Bregman divergence is the quantum relative entropy

$\mathcal{BD}_{h}(X||Y)=\textnormal{Tr}[X(\log X-\log Y)],$ (56)

the corresponding mirror map (regularized best-response) is the logit map,

$\displaystyle\textnormal{Mir}\Pi^{h}_{\mathcal{C}}(Y)=\frac{\exp(Y)}{\textnormal{Tr}[\exp(Y)]}=\Lambda(Y),$
(57)

and the corresponding proximal map is

$\displaystyle\textnormal{Prox}\Pi^{h}_{\mathcal{C}}(X,Y)=\frac{\exp(\log
X+Y)}{\textnormal{Tr}[\exp(\log X+Y)]}=\Lambda(\log X+Y).$ (58)

In Section 3, we will leverage these mirror and proximal steps (viewed as
oracles) to craft efficient iterative algorithms for finding Nash equilbria of
zero-sum quantum games.

##  3 Algorithms for Quantum Zero-Sum Games

As previously discussed, our goal is to efficiently find Nash equilibria,
$\Psi^{*}=(\alpha^{*},\beta^{*})$, satisfying the variational inequality of
Equation 20. This section describes existing algorithms, as well as our
proposed algorithm for finding these Nash equilibria. Section 4 will
demonstrate that our proposed algorithm achieves a quadratic speedup relative
to existing algorithms.

###  3.1 Matrix Dual Averaging (MDA) Method

In 2008, Jain and Watrous introduced a parallel algorithm to find approximate
Nash-equilibria of quantum zero-sum games [29]. Their algorithm was based on
the Matrix Multiplicative Weights Update (MMWU) method [46, 34, 70]. Here, we
will instead present and assess a more general method, which we will refer to
as Matrix Dual Averaging (MDA),666The name is a reference to the analogous
algorithm for classical games, i.e., Nesterov’s Dual Averaging [68]. Dual
Averaging is sometimes also known as “Lazy” Mirror Descent, since the key
difference from standard Mirror Descent is that averaging happens in the dual
space (rather than mapping from primal to dual space and back in each
iteration). which, with different DGFs, parameterizes a family of algorithms
that include Jain and Watrous’ MMWU algorithm.

MDA algorithms are played over several rounds, $t\in[0,N]$, of the quantum
zero-sum game. Suppose that, in round $t$, Alice and Bob play states
$\Psi_{t}=(\alpha_{t},\beta_{t})$.777Note that, in the first round ($t=0$),
Alice and Bob play the maximally mixed state
$\Psi_{0}=(\frac{\mathbb{1}_{n}}{2^{n}},\frac{\mathbb{1}_{m}}{2^{m}})$. After
performing his joint measurement, Roger returns feedback to Alice and Bob in
the form of their respective payoff gradients,
$\mathcal{F}_{\alpha}(\Psi_{t})$ and $\mathcal{F}_{\beta}(\Psi_{t})$, as
defined in Equation 16 and Equation 17. Alice and Bob use this feedback to
individually update their states so as to iteratively maximize their
individual payoff. The general MDA family of algorithms (Algorithm 4) can be
expressed concisely as the updates:

$\displaystyle\Phi_{t+1}$ $\displaystyle=\Phi_{t}+\mathcal{F}(\Psi_{t})$ (59)
$\displaystyle\Psi_{t+1}$
$\displaystyle=\textnormal{Mir}\Pi^{h}_{\mathcal{C}}(\eta\>\Phi_{t+1}),$ (60)

where $\eta$ is a chosen step-size and $h$ is a regularizing, convex DGF. As
we will see, it is the choice of DGF $h$ which instantiates different MDA
algorithms, including MMWU. Furthermore, relating back to the mirror step
picture of Figure 7, the intermediate states $\Phi$ are averaged in the dual
space, but are mapped onto feasible states $\Psi$ of the primal space by the
mirror projection $\textnormal{Mir}\Pi_{\mathcal{C}}^{h}$. Overall, the
algorithm only incorporates first-order derivatives in the update steps.

Algorithm 4 Matrix Dual Averaging

Accuracy Parameter: $\epsilon$

Regularization Function: $h:\\{\mathcal{A},\mathcal{B}\\}\to\mathbb{R}$

Diameter (of $\mathcal{BD}_{h}$):
$\mathcal{D}_{h}=\sup_{X,Y\in\mathcal{C}}\mathcal{BD}_{h}(X\|Y)$

Regularized Best Response:
$\textnormal{Mir}\Pi_{\mathcal{A},\mathcal{B}}^{h}(\cdot)=\arg\max_{x\in\\{\mathcal{A},\mathcal{B}\\}}\\{\langle\cdot,x\rangle-h(x)\\}$

$\eta\leftarrow\ {\mu_{h}}/(2\gamma_{\mathcal{F}})$ $\triangleright$ Step Size

$N\leftarrow\lceil\mathcal{D}_{h}/(\eta\cdot\epsilon^{2})\rceil$
$\triangleright$ Number of Rounds

$(\alpha_{0},\beta_{0})\leftarrow\left(\frac{1}{2^{n}}\mathbb{1}_{\mathcal{A}},\frac{1}{2^{m}}\mathbb{1}_{\mathcal{B}}\right)$
$\triangleright$ State Initialization

for $t\in[1,N-1]$ do

$\mathcal{CF}_{\alpha}^{(t)}\leftarrow\sum_{i=0}^{t-1}\mathcal{F}_{\alpha}(\Psi_{i})$
$\triangleright$ Calculate Cumulative Feedback

$\mathcal{CF}_{\beta}^{(t)}\leftarrow\sum_{i=0}^{t-1}\mathcal{F}_{\beta}(\Psi_{i})$

$\alpha_{t}\leftarrow\textnormal{Mir}\Pi_{\mathcal{A}}^{h}\left(\eta\>\mathcal{CF}_{\alpha}^{(t)}\right)$
$\triangleright$ State Updates

$\beta_{t}\leftarrow\textnormal{Mir}\Pi_{\mathcal{B}}^{h}\left(\eta\>\mathcal{CF}_{\beta}^{(t)}\right)$

end for

return
$\bigg{(}\bar{\alpha}=\frac{1}{N}\sum_{t=0}^{N-1}\alpha_{t},\bar{\beta}=\frac{1}{N}\sum_{t=0}^{N-1}\beta_{t}\bigg{)}$

We will now specify how each player chooses their next state via the
previously described MDA method. Note that both Alice and Bob perform the same
procedure, so we will only specify the procedure from Alice’s perspective.
Upon receiving feedback $\mathcal{F}_{\alpha}(\Psi_{t})$ from Roger, Alice
calculates her cumulative feedback,

$\mathcal{CF}_{\alpha}^{(t)}=\sum_{i=0}^{t-1}\mathcal{F}_{\alpha}(\Psi_{i}),$
(61)

which is simply the sum of the payoff gradients across all prior rounds. To
calculate the next state she will play, $\alpha_{t+1}$, Alice re-weights the
cummulative feedback by step-size $\eta$ and performs a mirror projection onto
her spectraplex,

$\alpha_{t+1}=\textnormal{Mir}\Pi^{h}_{\mathcal{A}}(\eta\>\mathcal{CF}_{\alpha}^{(t)})=\arg\max_{\alpha\in\mathcal{A}}\\{\langle\eta\>\mathcal{CF}_{\alpha}^{(t)},\alpha\rangle-h(\alpha)\\}.$
(62)

We will now specify how different choices of $h$ lead to different
instantiations of the MDA algorithm. For example, following from Lemma 6, if
the regularizer is the squared Frobenius norm,

$h(\alpha)=\frac{1}{2}\|\alpha\|_{F}^{2},$ (63)

then the mirror projection maps Alice’s next state onto the orthogonal
projection

$\alpha_{t+1}=\textrm{Orth}\Pi_{\mathcal{A}}(\eta\>\mathcal{CF}_{\alpha}^{(t)}).$
(64)

Meanwhile, following from Lemma 7, if the regularizer is the negative von
Neumann entropy,

$h(\alpha)=\textnormal{Tr}[\alpha\log\alpha],$ (65)

then the mirror projection maps Alice’s next state onto the logit map

$\alpha_{t+1}=\Lambda_{\mathcal{A}}\left(\eta\>\mathcal{CF}_{\alpha}^{(t)}\right)=\frac{\text{exp}\left(\eta\>\mathcal{CF}_{\alpha}^{(t)}\right)}{\textnormal{Tr}\left(\text{exp}\left(\eta\>\mathcal{CF}_{\alpha}^{(t)}\right)\right)}$
(66)

and, thus, the MMWU algorithm of [29]. However, as proven in [29],

######  Theorem 2 (MDA Rate)

MDA methods compute an $\epsilon$-Nash equilibrium in
$\mathcal{O}_{d}({1}/{\epsilon^{2}})$ steps.

Thus, MMWU for quantum zero-sum games obtains a worse convergence rate than
the best known algorithms for classical zero-sum games, which obtain an
$\mathcal{O}_{d}(1/\epsilon)$ rate, as described in Section 1.2 and depicted
in Figure 1. In the remainder of this work, we will show how the classically
efficient Mirror Prox algorithm and its variants can be lifted to the quantum
zero-sum game setting in order to obtain the desired
$\mathcal{O}_{d}(1/\epsilon)$ rate.

###  3.2 Matrix Mirror Prox (MMP)

The classical Mirror Prox (MP) algorithm is described in Section 1.2, with
updates given in Figure 1. In each step of the iterative MP optimization
procedure, the gradient of the players’ utilities is calculated followed by
the execution of a mirror step (Algorithm 2), as illustrated in Figure 7. The
mirror step can be implemented via a mirror map of the dual gradient, as in
Equation 44, or mirror map of the primal state, as in Equation 45. The
addition of proximal steps—or extra-gradient steps—is crucial as they
introduce a refinement to the basic gradient step. After the initial gradient
is computed and the first mirror step is taken, the algorithm does not
immediately proceed with this information. Instead, it calculates a second
gradient at the point reached after the first mirror step. This second
gradient offers a more accurate direction for the update because it accounts
for the immediate effects of the initial step, thus incorporating more
information about the curvature and constraints of the utility space. This
two-pronged gradient approach is what distinguishes the MP algorithm from
simpler gradient descent methods, offering a sophisticated tool for navigating
the intricate landscapes of variational inequalities.

Drawing inspiration from classical MP, we propose an analog for the quantum
zero-sum games of interest. This quantum analog, which we will refer to as the
Matrix Mirror Prox (MMP) method, has updates:

$\displaystyle\Phi_{t+1}$
$\displaystyle=\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{C}}(\Psi_{t},\mathcal{F}(\Psi_{t}))$
(67) $\displaystyle\Psi_{t+1}$
$\displaystyle=\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{C}}(\Psi_{t},\mathcal{F}(\Phi_{t+1})).$
(68)

Notice that these update steps are of the same form as the proximal step
updates given in Equation 48 and Equation 49. Alternatively, the updates can
be expressed purely in terms of the joint states $\Psi_{t}$:

$\displaystyle\Psi_{t+1}$
$\displaystyle=\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{C}}\bigg{(}\Psi_{t},\mathcal{F}\Big{(}\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{C}}(\Psi_{t},\mathcal{F}(\Psi_{t}))\Big{)}\bigg{)}.$
(69)

As discussed in Section 2.5, the incorporation of proximal steps in the Matrix
Mirror Prox algorithm offers two significant advantages over the conventional
Projected Subgradient Descent Ascent method:

  1. 1.

Geometry-Awareness: MMP replaces the Frobenius norm of PGDA with a generalized
Bregman divergence. Similar to MDA, utilizing a mirror map with this Bregman
divergence enables MMP to adapt to the inherent geometry of the problem space.
This makes the algorithm highly versatile, capable of effectively handling a
wide array of problems, including those with intricate constraints.

  2. 2.

Intermediate Proximal Step: Proximal steps in MMP introduce an intermediate
“prediction” state (i.e. state $Y_{t}$ in Section 2.5.2), which plays a
crucial role. At a high level, this prediction state extrapolates the future
trajectory of the descent/ascent, thereby guiding the subsequent state update.
This leverages the intuition that the gradient at the next state offers a more
informative update gradient than that at the current state.

Importantly, this combination of geometry-aware mirror maps and forward-
looking intermediate proximal steps accelerate the convergence of MMP with von
Neumann regularization to min-max points, thereby making MMP a potent tool for
solving min-max problems efficiently.

As is the case for classical MP (discussed in Section 1.2), a key downside of
MMP is that it requires two gradient calls per iteration. Specifically, in
each round $t$, to calculate the updated joint state $\Psi_{t+1}$, Alice and
Bob must query Roger with two distinct sets of states, $\Psi_{t}$ and
$\Phi_{t+1}$, so as to obtain the respective feedback gradients
$\mathcal{F}(\Psi_{t})$ and $\mathcal{F}(\Phi_{t+1})$. Classically, there
exist variants of MP with the same convergence rate that require only a single
gradient call per iteration. The remainder of this work will focus on a
proposed classically-inspired variant of MMP which requires only a single
gradient call per iteration, and a proof that it obtains the desired
$\mathcal{O}_{d}(1/\epsilon)$ convergence rate.

###  3.3 Optimistic Matrix Mirror-Prox (OMMP) Methods

Algorithm 5 Optimistic Matrix Mirror Prox (OMMP)

Accuracy Parameter: $\epsilon$

Regularization Function: $h:\\{\mathcal{A},\mathcal{B}\\}\to\mathbb{R}$

Strong Convexity Parameter (of $h$): $\mu_{h}$

Diameter (of $\mathcal{BD}_{h}$):
$\mathcal{D}_{h}=\sup_{X,Y\in\mathcal{C}}\mathcal{BD}_{h}(X\|Y)$

Lipschitz Parameter (of $\mathcal{F}$): $\gamma_{\mathcal{F}}$

Proximal Map:
$\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{A},\mathcal{B}}(X,Y):={\arg\min}_{C\in\mathcal{C}}\\{\langle
Y,C-X\rangle-\frac{1}{\eta}\mathcal{BD}_{h}(C||X)\\}.$

$\eta\leftarrow\ {\mu_{h}}/(2\gamma_{\mathcal{F}})$ $\triangleright$ Step Size

$N\leftarrow\lceil\mathcal{D}_{h}/(\eta\cdot\epsilon)\rceil$ $\triangleright$
Number of Rounds

$(\alpha_{0},\beta_{0})\leftarrow\left(\frac{1}{2^{n}}\mathbb{1}_{\mathcal{A}},\frac{1}{2^{m}}\mathbb{1}_{\mathcal{B}}\right)$
$\triangleright$ State Initialization

$(\hat{\alpha}_{0},\hat{\beta}_{0})\leftarrow(\alpha_{0},\beta_{0})$
$\triangleright$ Optimistic-Momentum Initialization

for $t\in[1,N-1]$ do

$\alpha_{t+1}\leftarrow\textnormal{Prox}\Pi_{\mathcal{A}}^{h,\eta}(\hat{\alpha}_{t},\mathcal{F}_{\alpha}(\Psi_{t}))$
$\triangleright$ State Updates

$\beta_{t+1}\leftarrow\textnormal{Prox}\Pi_{\mathcal{B}}^{h,\eta}(\hat{\beta}_{t},\mathcal{F}_{\beta}(\Psi_{t}))$

$\hat{\alpha}_{t+1}\leftarrow\textnormal{Prox}\Pi_{\mathcal{A}}^{h,\eta}(\hat{\alpha}_{t},\mathcal{F}_{\alpha}(\Psi_{t+1}))$
$\triangleright$ Momentum Updates

$\hat{\beta}_{t+1}\leftarrow\textnormal{Prox}\Pi_{\mathcal{B}}^{h,\eta}(\hat{\beta}_{t},\mathcal{F}_{\beta}(\Psi_{t+1}))$

end for

return
$\bigg{(}\bar{\alpha}=\frac{1}{N}\sum_{t=0}^{N-1}\alpha_{t},\bar{\beta}=\frac{1}{N}\sum_{t=0}^{N-1}\beta_{t}\bigg{)}$

We will now discuss the main contribution of this work, a novel approach to
finding Nash equilbria of quantum zero-sum games, which we refer to as the
Optimistic Matrix Mirror Prox (OMMP) method. This algorithm draws inspiration
from seminal classical game theory work on Optimistic Mirror Prox
(OMP)888There are a variety of classical algorithms of this general form,
known as Popov’s modified Arrow-Hurwicz algorithm, Past Extra-Gradient,
Optimistic Extra-Gradient, or the Optimistic Gradient Variant. methods, as
described in the Prior Work (Section 1.2). The OMMP algorithm, detailed in
Algorithm 5, can be succinctly represented through the following update
equations:

$\displaystyle\Psi_{t+1}$
$\displaystyle=\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{C}}(\Phi_{t},\mathcal{F}(\Psi_{t}))$
(70) $\displaystyle\Phi_{t+1}$
$\displaystyle=\textnormal{Prox}\Pi^{h,\eta}_{\mathcal{C}}(\Phi_{t},\mathcal{F}(\Psi_{t+1}))$
(71)

Here, $\eta$ is a predetermined step size, and $\textnormal{Prox}\Pi$ refers
to the proximal map, as elaborated in Algorithm 3 and Definition 7.

Although the update equations for OMMP closely resemble those of the standard
MMP algorithm, a critical difference exists. As alluded to in the previous
section, the OMMP variant necessitates only a single gradient call per
iteration. While it may appear that both Alice and Bob require two feedback
gradients, $\mathcal{F}(\Psi_{t})$ and $\mathcal{F}(\Psi_{t+1})$, for each
round $t$, they will have already acquired $\mathcal{F}(\Psi_{t})$ during
round $t-1$. This allows them to reuse $\mathcal{F}(\Psi_{t})$ and only make a
query for a new gradient with one state, $\Psi_{t+1}$.

To gain some intuition for this algorithm, we will revisit the original
proposal of “optimism” in this context, due to Popov [52]. In conventional
Grdient Descent Ascent, both players update their strategies concurrently,
moving in directions opposite to their individual gradients—akin to a tug-of-
war. Optimistic Gradient Descent Ascent (OGDA) enhances this by incorporating
a “look-ahead” feature. Players consider not only the current gradient but
also anticipate the opponent’s forthcoming strategic adjustments. This
anticipatory or “optimistic” approach is formalized by adding an extra term to
the update rule, essentially an estimate of the future gradient.

This anticipatory element is also how OMMP diverges from standard MMP. Rather
than employing a new prediction, optimistic OMMP utilizes a past “inertia”
point, $\Phi_{t}$. Unlike MMP, which focuses solely on future states, OMMP
integrates past gradient information to stabilize its learning trajectory.

##  4 Convergence Analysis

We will now prove one of the main results of the paper: See 1 This establishes
that our proposed OMMP method yields a quadratic speedup compared to previous
MDA-based approaches. Concretely, we will leverage the proof structure of Ene
and Nguyên [56] to demonstrate that, with an appropriately chosen step-size
$\eta$ and for a given accuracy parameter $\epsilon$, the _average iterate_

$\displaystyle\overline{\Psi}_{T}=\frac{1}{T}\sum_{t=1}^{T}\Psi_{t}$ (72)

converges to an $\epsilon$-Nash equilibrium in $O_{d}(1/\epsilon)$
iterations.999 It’s important to note the difference in the notation used to
describe algorithmic convergence rates between the optimization and computer
science communities. In the realm of optimization, _convergence rates_ such as
$\mathcal{O}(e^{-\rho T})$, $\mathcal{O}(1/T)$, and $\mathcal{O}(1/\sqrt{T})$
are commonly employed to indicate how swiftly an algorithm approximates the
optimal solution as the iteration count $T$ increases. In contrast, the
computer science literature tends to recast these convergence rates as number
of iterations required to achieve an $\epsilon$-accurate solution, thereby
giving a more direct account of computational efficiency. Expressed in terms
of accuracy parameter $\epsilon$, this leads to expressions like
$\mathcal{O}(\log(1/\epsilon))$, $\mathcal{O}(1/\epsilon)$, and
$\mathcal{O}(1/\epsilon^{2})$.

We will conclude with a discussion of the dimensionality-dependence hidden in
the $\mathcal{O}_{d}(\cdot)$ notation. In particular in Section 4.7, we will
show that the von Neumann entropy regularizer reduces OMMP’s dimensionality-
dependence logarithmically relative to the Frobenius norm regularizer. Since
OMMP with the von Neumann entropy regularizer instatiates the OMMWU algorithm,
we thus prove that OMMWU obtains the desired $\mathcal{O}(d/\epsilon)$
iteration complexity and, thus, achieves a quadratic speedup to MMWU’s
$\mathcal{O}(d/\epsilon^{2})$ iteration complexity.

###  4.1 Monotonicity of $\mathcal{F}$

Critical to our analysis is the monotonicity of the feedback operator
$\mathcal{F}$, which is a key property of the finite-valued quantum zero sum
games of interest.

######  Lemma 8 (Mononicity of $\mathcal{F}$)

$\mathcal{F}$ is monotone or, equivalently,

$\displaystyle\textnormal{Tr}\left[{(\mathcal{F}(X)-\mathcal{F}(Y))}\>{(X-Y)}\right]\geq
0,\>\>\forall X,Y\in\mathcal{C}.$ (73)

The proof of this result employs techniques involving the Pauli decomposition
of the feedback operator, thereby deviating from the conventional proof for
monotonicity of classical normal game operators. We defer the proof of
$\mathcal{F}$’s monotonicity (in addition to many other properties) to Section
C.2.

###  4.2 Error Decomposition

Following from our discussion of Nash equilibria in Section 2.2 and the
natural error function, given in Equation 26, the error of the average iterate
is defined as

$\displaystyle\textnormal{Error}(\overline{\Psi}_{T})$
$\displaystyle=\sup_{Z\in\mathcal{C}}\textnormal{Tr}\left[{\mathcal{F}(Z)}\>{(\overline{\Psi}_{T}-Z)}\right]=\frac{1}{T}\sup_{Z\in\mathcal{C}}\>\sum_{t=1}^{T}\textnormal{Tr}\left[{\mathcal{F}(Z)}\>{(\Psi_{t}-Z)}\right].$
(74)

From the monotonicity of $\mathcal{F}$, given in Equation 73, for $\Psi_{t}$
and some state $Z\in\mathcal{C}$,

$\displaystyle\textnormal{Tr}\left[{\mathcal{F}(Z)}\>{(\Psi_{t}-Z)}\right]\leq\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Psi_{t}-Z)}\right],$
(75)

which implies that the average-iterate error is upper bounded as

$\displaystyle\textnormal{Error}(\overline{\Psi}_{T})\leq\frac{1}{T}\sup_{Z\in\mathcal{C}}\>\sum_{t=1}^{T}\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Psi_{t}-Z)}\right].$
(76)

In order to produce an explicit bound, we begin by decomposing the terms of
the summation, as follows, resulting in a sum of three distinct terms:
$(A_{t})$, $(B_{t})$, and $(C_{t})$.

$\displaystyle\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Psi_{t}-Z)}\right]=\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Psi_{t}-\Phi_{t}+\Phi_{t}-Z)}\right]$
$\displaystyle=\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Psi_{t}-\Phi_{t})+\mathcal{F}(\Psi_{t})(\Phi_{t}-Z)}\right]=\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})(\Psi_{t}-\Phi_{t})}\>{+}\right]\textnormal{Tr}\left[{(\Psi_{t})}\>{(\Phi_{t}-Z)}\right]$
$\displaystyle=\textnormal{Tr}\left[{(\mathcal{F}(\Psi_{t})-\mathcal{F}(\Psi_{t-1})+\mathcal{F}(\Psi_{t-1}))}\>{(\Psi_{t}-\Phi_{t})}\right]+\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Phi_{t}-Z)}\right]$
$\displaystyle=\textnormal{Tr}\left[{(\mathcal{F}(\Psi_{t})-\mathcal{F}(\Psi_{t-1}))}\>{(\Psi_{t}-\Phi_{t})+\mathcal{F}(\Psi_{t-1})(\Psi_{t}-\Phi_{t})}\right]+\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Phi_{t}-Z)}\right]$
$\displaystyle=\underbrace{\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Phi_{t}-Z)}\right]}_{(A_{t})}+\underbrace{\textnormal{Tr}\left[{(\mathcal{F}(\Psi_{t})-\mathcal{F}(\Psi_{t-1}))}\>{(\Psi_{t}-\Phi_{t})}\right]}_{(B_{t})}+\underbrace{\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t-1})}\>{(\Psi_{t}-\Phi_{t})}\right]}_{(C_{t})}.$
(77)

In the following analysis, we will show that terms $(A_{t})$ and $(C_{t})$ can
be upper bounded via the optimality of the proximal updates in the OMMP
algorithm. Furthermore, the term $(B_{t})$ can be upper bounded via the strong
convexity, smoothness, and duality of the OMMP proximal updates.

###  4.3 Optimality Analysis of $(A_{t})$

We begin by upper bounding the term $(A_{t})$. Recall from the definition of
the OMMP algorithm in Section 3.3 that $\Phi_{t}$ is the outcome of the
proximal projection of Equation 71, for a regularization function $h$.
Therefore, $\Phi_{t}$ can be expressed as

$\displaystyle\Phi_{t}=\arg\min_{Z\in\mathcal{C}}\phi_{t}(Z),$ (78)

where $\phi_{t}(Z)$ is the proximal optimization term

$\displaystyle\phi_{t}(Z)=\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(Z-\Phi_{t-1})}\right]+\frac{1}{\eta}\mathcal{BD}_{h}(Z\|\Phi_{t-1}).$
(79)

Note that $\phi_{t}$ is differentiable, with gradient

$\displaystyle\nabla_{Z^{\top}}\phi_{t}(Z)$
$\displaystyle=\mathcal{F}(\Psi_{t})+\frac{1}{\eta}\nabla_{Z^{\top}}\mathcal{BD}_{h}(Z\|\Phi_{t-1}).$
(80)

Using the definitions of Section 2.3, in the supplement (Appendix D) we
further show that:

###### Proposition 3

For $\mu$-strongly convex regularization function $h$ with respect to norm
$\|\cdot\|$, $\phi_{t}$ is $\frac{\mu}{\eta}$-strongly convex with respect to
norm $\|\cdot\|$.

Since $\phi_{t}$ is strongly convex and differentiable it satisfies the First-
Order Optimality Condition (1) for feasible solution momentum $\Phi_{t}$,

$\displaystyle\textnormal{Tr}\left[{\nabla_{\Phi_{t}^{\top}}\phi_{t}(\Phi_{t})}\>{(Z-\Phi_{t})}\right]\geq
0,\>\>\forall Z\in\mathcal{C}.$ (81)

Plugging the gradient of Equation 80 into this expression gives

$\displaystyle\textnormal{Tr}\left[{\Big{(}\mathcal{F}(\Psi_{t})+\frac{1}{\eta}\nabla_{\Phi_{t}^{\top}}\mathcal{BD}_{h}(\Phi_{t}||\Phi_{t-1})\Big{)}}\>{(\Phi_{t}-Z)}\right]\leq
0$
$\displaystyle\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Phi_{t}-Z)}\right]+\frac{1}{\eta}\textnormal{Tr}\left[{\nabla_{\Phi_{t}^{\top}}\mathcal{BD}_{h}(\Phi_{t}||\Phi_{t-1})}\>{(\Phi_{t}-Z)}\right]\leq
0$
$\displaystyle\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Phi_{t}-Z)}\right]\leq\frac{1}{\eta}\textnormal{Tr}\left[{-\nabla_{\Phi_{t}^{\top}}\mathcal{BD}_{h}(\Phi_{t}||\Phi_{t-1})}\>{(\Phi_{t}-Z)}\right].$
(82)

From our definition of $(A_{t})$ in Equation 77 and the three-point identity
for the Bregman divergence given in Equation 42, it follows that, $\forall
Z\in\mathcal{C}$,

$\displaystyle(A_{t})$
$\displaystyle\leq\frac{1}{\eta}\Big{(}\mathcal{BD}_{h}(Z||\Phi_{t-1})-\mathcal{BD}_{h}(\Phi_{t}||\Phi_{t-1})-\mathcal{BD}_{h}(Z||\Phi_{t})\Big{)}.$
(83)

###  4.4 Optimality Analysis of $(C_{t})$

We will now upper bound $(C_{t})$, similarly to how we upper bounded
$(A_{t})$. Recall from the definition of the OMMP algorithm in Section 3.3
that $\Psi_{t}$ is the outcome of the proximal projection of Equation 70, for
a regularization function $h$. Therefore, $\Psi_{t}$ can be expressed as

$\displaystyle\Psi_{t}=\arg\min_{Z\in\mathcal{C}}\psi_{t}(Z),$ (84)

where $\psi_{t}(Z)$ is the proximal optimization term

$\displaystyle\psi_{t}(Z)=\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t-1})}\>{(Z-\Phi_{t-1})}\right]+\frac{1}{\eta}\mathcal{BD}_{h}(Z\|\Phi_{t-1}).$
(85)

Note that $\phi_{t}$ is differentiable, with gradient

$\displaystyle\nabla_{Z^{\top}}\psi_{t}(Z)$
$\displaystyle=\mathcal{F}(\Psi_{t-1})+\frac{1}{\eta}\nabla_{Z^{\top}}\mathcal{BD}_{h}(Z\|\Phi_{t-1}).$
(86)

Similarly with the previous section, in the supplement (Appendix D) we further
show that:

###### Proposition 4

For $\mu$-strongly convex regularization function $h$ with respect to norm
$\|\cdot\|$, $\psi_{t}$ is $\frac{\mu}{\eta}$-strongly convex with respect to
norm $\|\cdot\|$.

Since $\psi_{t}$ is strongly convex and differentiable, it satifies the First-
Order Optimality Condition (1) for a feasible solution state $\Psi_{t}$,

$\displaystyle\textnormal{Tr}\left[{\nabla_{\Psi_{t}^{\top}}\psi_{t}(\Psi_{t})}\>{(Z-\Psi_{t})}\right]\geq
0,\>\>\forall Z\in\mathcal{C}.$ (87)

Plugging the gradient of Equation 86 into this expression gives, $\forall
Z\in\mathcal{C}$,

$\displaystyle\textnormal{Tr}\left[{\Big{(}\mathcal{F}(\Psi_{t-1})+\frac{1}{\eta}\nabla_{\Psi_{t}^{\top}}\mathcal{BD}_{h}(\Psi_{t}||\Phi_{t-1})\Big{)}}\>{(\Psi_{t}-Z)}\right]\leq
0$ (88)
$\displaystyle\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t-1})}\>{(\Psi_{t}-Z)}\right]\leq\frac{1}{\eta}\textnormal{Tr}\left[{\nabla_{\Psi_{t}^{\top}}\mathcal{BD}_{h}(\Psi_{t}||\Phi_{t-1})}\>{(Z-\Psi_{t})}\right].$
(89)

Since the expression holds for all $Z\in\mathcal{C}$, we can choose to set
$Z=\Phi_{t}$, which implies

$\displaystyle\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t-1})}\>{(\Psi_{t}-\Phi_{t})}\right]\leq\frac{1}{\eta}\textnormal{Tr}\left[{\nabla_{\Psi_{t}^{\top}}\mathcal{BD}_{h}(\Psi_{t}||\Phi_{t-1})}\>{(\Phi_{t}-\Psi_{t})}\right].$
(90)

From our definition of $(C_{t})$ in Equation 77 and the Three-Point Tdentity
for the Bregman divergence, given in Equation 42, it follows that

$\displaystyle(C_{t})\leq\frac{1}{\eta}\Big{(}\mathcal{BD}_{h}(\Phi_{t}||\Phi_{t-1})-\mathcal{BD}_{h}(\Psi_{t}||\Phi_{t-1})-\mathcal{BD}_{h}(\Phi_{t}||\Psi_{t})\Big{)}.$
(91)

###  4.5 An Upper Bound for $(B_{t})$

By the Cauchy-Schwarz inequality for Frobenius norm, we can upper bound
$(B_{t})$ as

$\displaystyle(B_{t})$
$\displaystyle=\textnormal{Tr}\left[{(F(\Psi_{t})-F(\Psi_{t-1}))}\>{(\Psi_{t}-\Phi_{t})}\right]\leq\|F(\Psi_{t})-F(\Psi_{t-1})\|_{F}\|\Psi_{t}-\Phi_{t}\|_{F}.$
(92)

or more generally, with respect to norms $(\|\cdot\|,\|\cdot\|_{*})$

$\displaystyle(B_{t})$
$\displaystyle=\textnormal{Tr}\left[{(F(\Psi_{t})-F(\Psi_{t-1}))}\>{(\Psi_{t}-\Phi_{t})}\right]\leq\|F(\Psi_{t})-F(\Psi_{t-1})\|_{*}\|\Psi_{t}-\Phi_{t}\|.$
(93)

To further manipulate this upper bound for $(B_{t})$, we leverage the ideas
from Fenchel duality and convexity introduced in Section 2.3.

Following from Equation 84 and Equation 79, respectively, $\Psi_{t}$ is the
minimizer of $\psi_{t}$ and $\Phi_{t}$ is the minimizer of $\phi_{t}$. Re-
expressing $\phi_{t}$ in terms of $\psi_{t}$,

$\displaystyle\phi_{t}(U)=\psi_{t}(U)+\textnormal{Tr}\left[{\left(\mathcal{F}(\Psi_{t})-\mathcal{F}(\Psi_{t-1})\right)}\>{(U-\Phi_{t-1})}\right]$
(94)

implies that $\Phi_{t}$ can be expressed in terms of a minimization of
$\psi_{t}$,

$\displaystyle\Phi_{t}=\arg\min_{U\in\mathcal{C}}\left\\{\psi_{t}(U)+\textnormal{Tr}\left[{\left(\mathcal{F}(\Psi_{t})-\mathcal{F}(\Psi_{t-1})\right)}\>{U}\right]\right\\}.$
(95)

Since we assumed that $h$ is $\mu$-strongly convex, by Proposition 4
$\psi_{t}$ must be $\tfrac{\mu}{\eta}$-strongly convex. Denoting the Fenchel
conjugate of $\psi_{t}$ as $\psi_{t}^{*}$, since $\psi_{t}$ is strongly
convex, Danskin’s Theorem (Proposition 2) states that $\nabla\psi_{t}^{*}$ is
well-defined, for all $V\in\mathcal{C}$, as

$\displaystyle\nabla\psi^{*}_{t}(V)=\arg\min_{U\in\mathcal{C}}\left\\{\psi_{t}(U)-\textnormal{Tr}\left[{V}\>{U}\right]\right\\}.$
(96)

In the case that $V=0$, Equation 78 implies that

$\displaystyle\nabla\psi^{*}_{t}(0)=\arg\min_{U\in\mathcal{C}}\psi_{t}(U)=\Psi_{t}.$
(97)

Meanwhile, if $V=\mathcal{F}(\Psi_{t-1})-\mathcal{F}(\Psi_{t})$, Equation 95
implies that

$\displaystyle\nabla\psi^{*}_{t}\left(\mathcal{F}(\Psi_{t-1})-\mathcal{F}(\Psi_{t})\right)=\arg\min_{U\in\mathcal{C}}\left\\{\psi_{t}(U)+\textnormal{Tr}\left[{\left(\mathcal{F}(\Psi_{t})-\mathcal{F}(\Psi_{t-1})\right)}\>{U}\right]\right\\}=\Phi_{t}.$
(98)

Therefore,

$\displaystyle\|\Psi_{t}-\Phi_{t}\|_{F}$
$\displaystyle=\|\nabla\psi^{*}_{t}(0)-\nabla\psi^{*}_{t}\left(\mathcal{F}(\Psi_{t-1})-\mathcal{F}(\Psi_{t})\right)\|_{F}.$
(99)

Furthermore, since $\psi_{t}$ is $\frac{\mu}{\eta}$-strongly convex with
respect the norm $\|\cdot\|$, Strong-Smooth Duality (Lemma 3) implies that
$\psi^{*}_{t}$ is $\frac{\eta}{\mu}$-smooth with respect the norm
$\|\cdot\|_{*}$. By Proposition 1, this implies that $\nabla\psi^{*}_{t}$ is
$\frac{\eta}{\mu}$-Lipschitz continuous with respect $\|\cdot\|_{*}$,

$\displaystyle\|\nabla\psi^{*}_{t}(X)-\nabla\psi^{*}_{t}\left(Y\right)\|=\|\nabla\psi^{*}_{t}(X)-\nabla\psi^{*}_{t}\left(Y\right)\|_{**}$
$\displaystyle\leq\frac{\eta}{\mu}\|X-Y\|_{*}.$ (100)

Equivalently for a pair of norms, we get Plugging this and Equation 99 into
Equation 92 produces the final $(B_{t})$ upper bound,

$(B_{t})\leq\frac{\eta}{\mu}\|\mathcal{F}(\Psi_{t})-\mathcal{F}(\Psi_{t-1})\|_{*}^{2}.$
(101)

###  4.6 Average-Iterate Error Upper Bound

Combining the average-iterate error upper bound of Equation 76 and its
decomposition, as expressed in Equation 77, we have that

$\displaystyle\textnormal{Error}(\overline{\Psi}_{T})$
$\displaystyle\leq\frac{1}{T}\sup_{Z\in\mathcal{C}}\sum_{t=1}^{T}\textnormal{Tr}\left[{\mathcal{F}(\Psi_{t})}\>{(\Psi_{t}-Z)}\right]=\frac{1}{T}\sup_{Z\in\mathcal{C}}\sum_{t=1}^{T}(A_{t})+(B_{t})+(C_{t}).$

Plugging in the upper bounds for $(A_{t})$ from Equation 83, $(B_{t})$ from
Equation 101, and $(C_{t})$ from Equation 91 gives

$\displaystyle\textnormal{Error}(\overline{\Psi}_{T})$
$\displaystyle\leq\frac{1}{T}\sup_{Z\in\mathcal{C}}\sum_{t=1}^{T}\bigg{(}\frac{1}{\eta}\Big{(}\mathcal{BD}_{h}(Z\|\Phi_{t-1})-\left(\mathcal{BD}_{h}(Z\|\Phi_{t})+\mathcal{BD}_{h}(\Psi_{t}\|\Phi_{t-1})+\mathcal{BD}_{h}(\Phi_{t}\|\Psi_{t})\right)\Big{)}$
$\displaystyle\hskip
101.17755pt+\frac{\eta}{\mu}\|\mathcal{F}(\Psi_{t})^{\top}-\mathcal{F}(\Psi_{t-1})^{\top}\|_{*}^{2}\bigg{)}$
$\displaystyle\leq\frac{1}{T}\sup_{Z\in\mathcal{C}}\underbrace{\frac{1}{\eta}\sum_{t=1}^{T}\big{(}\mathcal{BD}_{h}(Z\|\Phi_{t-1})-\mathcal{BD}_{h}(Z\|\Phi_{t})\big{)}}_{(I)}-\frac{1}{\eta}\sum_{t=1}^{T}\big{(}\mathcal{BD}_{h}(\Psi_{t}\|\Phi_{t-1})+\mathcal{BD}_{h}(\Phi_{t}\|\Psi_{t})\big{)}$
$\displaystyle\hskip
101.17755pt+\frac{\eta}{\mu}\sum_{t=1}^{T}\|\mathcal{F}(\Psi_{t})^{\top}-\mathcal{F}(\Psi_{t-1})^{\top}\|_{*}^{2}.$

Via telescoping sums, the first term $(I)$ can be reduced to

$\displaystyle(I)=\frac{1}{\eta}\sum_{t=1}^{T}\big{(}\mathcal{BD}_{h}(Z\|\Phi_{t-1})-\mathcal{BD}_{h}(Z\|\Phi_{t})\big{)}=\frac{1}{\eta}\Big{(}\mathcal{BD}_{h}(Z\|\Phi_{0})-\mathcal{BD}_{h}(Z\|\Phi_{T})\Big{)}.$
(102)

Additionally, since $h$ is $\mu$-strongly convex with respect $\|\cdot\|$,

$\displaystyle\mathcal{BD}_{h}(X\|Y)\geq\frac{\mu}{2}\|X-Y\|^{2}.$ (103)

Therefore, we can re-express the error as

$\displaystyle\textnormal{Error}(\overline{\Psi}_{T})$
$\displaystyle\leq\frac{1}{T}\sup_{Z\in\mathcal{C}}\frac{1}{\eta}\mathcal{BD}_{h}(Z\|\Phi_{0})\underbrace{-\frac{\mu}{2\eta}\|Z-\Phi_{T}\|^{2}+-\frac{\mu}{2\eta}\sum_{t=1}^{T}\big{(}\|\Phi_{t-1}-\Psi_{t}\|^{2}+\|\Phi_{t}-\Psi_{t}\|^{2}\big{)}}_{(II)}$
$\displaystyle\hskip
101.17755pt+\underbrace{\frac{\eta}{\mu}\sum_{t=1}^{T}\|\mathcal{F}(\Psi_{t})^{\top}-\mathcal{F}(\Psi_{t-1})^{\top}\|_{*}^{2}}_{(III)},$

and proceed to upper bound terms $(II)$ and $(III)$. We begin by expanding
$(II)$ as

$\displaystyle(II)$
$\displaystyle=-\frac{\mu}{2\eta}\|\Phi_{T}-Z\|^{2}-\frac{\mu}{2\eta}\sum_{t=1}^{T}\big{(}\|\Phi_{t-1}-\Psi_{t}\|^{2}+\|\Phi_{t}-\Psi_{t}\|^{2}\big{)}$
$\displaystyle=\frac{\mu}{2\eta}\bigg{(}\|\Phi_{0}-\Psi_{0}\|^{2}-\|\Phi_{T}-\Psi_{T}\|^{2}-\|\Phi_{T}-Z\|^{2}-\sum_{t=1}^{T}\big{(}\|\Phi_{t-1}-\Psi_{t}\|^{2}+\|\Phi_{t-1}-\Psi_{t-1}\|^{2}\big{)}\bigg{)}$

Since we initialized $\Phi_{0}=\Psi_{0}$,

$\displaystyle\|\Phi_{0}-\Psi_{0}\|^{2}=0.$ (104)

Furthermore, since $\|\Phi_{T}-Z\|^{2}\geq 0$ for all $Z\in\mathcal{C}$,
including $Z=\Psi_{T}$,

$\displaystyle(II)$
$\displaystyle\leq-\frac{\mu}{2\eta}\sum_{t=1}^{T}\big{(}\|\Phi_{t-1}-\Psi_{t}\|^{2}+\|\Phi_{t-1}-\Psi_{t-1}\|^{2}\big{)}.$
(105)

Defining

$\displaystyle\Gamma:=\sum_{t=1}^{T}\left(\|\Psi_{t}-\Phi_{t-1}\|^{2}+\|\Psi_{t-1}-\Phi_{t-1}\|^{2}\right)$
(106)

implies the following upper bound on $(II)$

$\displaystyle(II)$ $\displaystyle\leq-\frac{\mu}{2\eta}\Gamma.$ (107)

In order to bound term $(III)$, we will leverage an important property of
finite-valued quantum zero-sum games. Namely, in Appendix D, we prove several
properties of the $\mathcal{F}$ operator, including the following lemma
regarding Lipschitz continuity.

######  Lemma 9 (Lipschitz Continuity of $\mathcal{F}$)

For finite-valued quantum zero-sum games, $\mathcal{F}$ is a Lipschitz
continuous operator, meaning there exists a Lipschitz constant
$\gamma\in\mathbb{R}$ such that

$\displaystyle\|\mathcal{F}(Z)-\mathcal{F}(Z^{\prime})\|_{*}\leq\gamma\|Z-Z^{\prime}\|,\>\>\forall\>Z,Z^{\prime}\in\mathcal{C}.$
(108)

For the case of $(\|\cdot\|_{F},\|\cdot\|_{F})$ and
$(\|\cdot\|_{\infty},\|\cdot\|_{1})$, $\gamma_{F,F}=\mathcal{O}(2^{d})$ while
$\gamma_{\infty,1}=\mathcal{O}(1)$.

Leveraging the Lipschitz continuity of $\mathcal{F}$, as given in Equation
108, and the fact that for an arbitrary norm $\|\cdot\|$

$\displaystyle\|X-Y\|^{2}\leq(\|X\|+\|Y\|)^{2}\leq\|X\|^{2}+\|Y\|^{2}+2\|X\|\|Y\|\leq
2\|X\|^{2}+2\|Y\|^{2},$ (109)

$(III)$ can be bounded as

$\displaystyle(III)$
$\displaystyle=\frac{\eta}{\mu}\sum_{t=1}^{T}\|\mathcal{F}(\Psi_{t})-\mathcal{F}(\Psi_{t-1})\|_{*}^{2}\leq\frac{\eta}{\mu}\sum_{t=1}^{T}\gamma^{2}\|\Psi_{t}-\Psi_{t-1}\|^{2}$
$\displaystyle=\frac{\gamma^{2}\eta}{\mu}\sum_{t=1}^{T}\|\Psi_{t}-\Phi_{t-1}+\Phi_{t-1}-\Psi_{t-1}\|^{2}$
$\displaystyle=\frac{\gamma^{2}\eta}{\mu}\sum_{t=1}^{T}\|(\Psi_{t}-\Phi_{t-1})-(\Psi_{t-1}-\Phi_{t-1})\|^{2}$
$\displaystyle\leq\frac{2\gamma^{2}\eta}{\mu}\Big{(}\sum_{t=1}^{T}\|\Psi_{t}-\Phi_{t-1}\|^{2}+\|\Psi_{t-1}-\Phi_{t-1}\|^{2}\Big{)}$
$\displaystyle=\frac{2\gamma^{2}\eta}{\mu}\>\Gamma.$ (110)

Therefore, combining the bound on $(II)$, from Equation 107, with that of
$(III)$, from Section 4.6, results in

$\displaystyle(II)+(III)$
$\displaystyle\leq\left(\frac{2\gamma^{2}\eta}{\mu}-\frac{\mu}{2\eta}\right)\Gamma.$
(111)

Note that there exist step-sizes $\eta$, such that $(II)+(III)\leq 0$. Namely,

$\displaystyle\eta\leq\frac{\mu}{2\gamma}\>\>\implies\>\>(II)+(III)\leq 0.$
(112)

Recall that our objective was to upper bound
$\textnormal{Error}(\overline{\Psi}_{T})$, where

$\displaystyle\textnormal{Error}(\overline{\Psi}_{T})$
$\displaystyle\leq\frac{1}{T}\sup_{Z\in\mathcal{C}}\frac{1}{\eta}\mathcal{BD}_{h}(\Phi_{0}\|Z)+(II)+(III),$
(113)

If we choose the step-size such that $\eta\leq\frac{\mu}{2\gamma}$, the upper
bound simplifies to

$\textnormal{Error}(\overline{\Psi}_{T})\leq\frac{1}{T}\sup_{Z\in\mathcal{C}}\frac{1}{\eta}\mathcal{BD}_{h}(\Phi_{0}\|Z).$
(114)

Since $\mathcal{BD}_{h}(\Phi_{0}\|Z)$ is bounded above by the Bregman diameter
of $\mathcal{C}$, i.e.

$\displaystyle\mathcal{BD}_{h}(\Phi_{0}\|Z)\leq\sup_{X,Y\in\mathcal{C}}\mathcal{BD}_{h}(X\|Y),$
(115)

it has no time dependence, meaning

$\displaystyle\textnormal{Error}(\overline{\Psi}_{T})=\mathcal{O}_{d}\bigg{(}\frac{1}{T}\bigg{)},$

where $\mathcal{O}_{d}$ is used to refer to the convergence dependence on the
spectraplex dimension $d$, governed by the Bregman diameter of $\mathcal{C}$,
as given in Equation 115. Therefore, setting
$\textnormal{Error}(\overline{\Psi}_{T})$ to a desired accuracy parameter,
i.e. $\textnormal{Error}(\overline{\Psi}_{T})=\epsilon$, we obtain that the
number of steps necessary for average-iterate convergence to an
$\epsilon$-Nash equilibrium is

$\displaystyle T=\mathcal{O}_{d}\left(\frac{1}{\epsilon}\right).$ (116)

###  4.7 Regularizer-Induced Dimensionality Dependence

Following from this analysis, for a fixed dimension $d$, the family of
algorithms encompassed by OMMP can compute an $\epsilon$-Nash Equilibrium in
$\mathcal{O}_{d}(1/\epsilon)$ steps. Note, however, that the
$\mathcal{O}_{d}(\cdot)$ notation hides the dimension-dependence of the
result, incurred by the $\frac{1}{\eta}\mathcal{BD}_{h}(\Phi_{0}\|Z)$ term. In
order to achieve an explicit bound, we will demonstrate how different OMMP
regularizers affect the Bregman divergence $\mathcal{BD}_{h}(\Phi_{0}\|Z)$ and
step-size $\eta$. In particular, we consider two regularizers—the Frobenius
norm and the von Neumann entropy—which are quantum analogs of two of the most
widely used classical regularizers (entropy and Euclidean norm). Note that
OMMP instantiated with the Frobenius norm implements the Optimistic Matrix
Extra-Gradient (OMEG) algorithm (updates given in Figure 2), while the von
Neumann entropy regularizer implements the Optimistic Matrix Multiplicative
Weights Update (OMMWU) algorithm.

We will first consider the Bregman divergence term, noting that in OMMP we
initialize
$\Phi_{0}=\left(\frac{1}{2^{n}}\mathbb{1}_{\mathcal{A}},\frac{1}{2^{m}}\mathbb{1}_{\mathcal{B}}\right)$.
We offer an explicit proof for the Frobenius norm regularizer:

###### Proposition 5

The Frobenius norm regularizer induces a maximum Bregman divergence that is
constant, i.e.
$\max_{X\in\mathcal{C}}\mathcal{BD}_{h}(\Phi_{0}||X)=2=\mathcal{O}(1)$.

#####  Proof.

Using the derivation of the Bregman divergence for the Frobenius norm in
Equation 52 and the fact that $\forall X\in\mathcal{C}$,
$\textnormal{Tr}(X)=\textnormal{Tr}(\alpha)+\textnormal{Tr}(\beta)\leq 1+1=2$,

$\displaystyle\max_{X\in\mathcal{C}}\mathcal{BD}_{h}(\Phi_{0}||X)$
$\displaystyle=\max_{X\in\mathcal{C}}\tfrac{1}{2}||\Phi_{0}-X||_{F}^{2}=\tfrac{1}{2}\max_{X\in\mathcal{C}}\left(||\Phi_{0}||_{F}^{2}-2\textnormal{Tr}(\Phi_{0}^{\dagger}X)+||X||_{F}^{2}\right)$
$\displaystyle\leq\tfrac{1}{2}\max_{X\in\mathcal{C}}\left(2+||X||_{F}^{2}\right)\leq\tfrac{1}{2}\left(2+2\right)=2.$

$\square$
In the case of the von Neumann entropy regularizer, as discussed in [71], for
a $(4^{d}-1)$-dimensional spectraplex, the corresponding bound on the Bregman
divergence is equal to the maximum quantum entropy for a $2d$-qubit system,
and hence linear in $d$—more precisely, equal to
$\log(4^{d}-1)=\mathcal{O}(d)$.

We now consider the step-size. Following from Equation 112, the optimal step-
size is dependent on $h$’s strong convexity parameter ($\mu_{h}$) and
$\mathcal{F}$’s Lipschitz constant ($\gamma_{\mathcal{F}}$), as

$\displaystyle\eta=\Theta\left({\mu_{h}/\gamma_{\mathcal{F}}}\right).$ (117)

In scenarios utilizing the Frobenius regularizer, the modulus of strong
convexity is 1, while $\gamma_{\text{Lipsch}}$, the Lipschitz constant with
respect to $\|\cdot\|_{F}$ of the gradient for games with bounded utilities in
$[-1,1]$, is of the order of $\mathcal{O}(\sqrt{4^{d}})=\mathcal{O}(2^{d})$.
In contrast, for the von Neumann entropy, the modulus of strong convexity is 1
with respect to the Shatten 1-norm (the matrix analog of the $\ell_{1}$-norm)
[72]. Moreover, for games with bounded utilities, the relative Lipschitz
constant with respect to Shatten 1-norm $\|\cdot\|_{1}$ is of the order of
$\mathcal{O}(1)$.

In light of the preceding discussions, we conclude with the dimension-
dependent iteration complexities of OMMWU and OMEG for achieving
$\epsilon$-Nash equilibria of quantum zero-sum games:

See 1

######  Corollary 2 (OMEG Iteration Complexity)

In the $4^{d}$-dimensional spectraplex, OMEG computes average-iterate
$\epsilon$-Nash equilibria in $\mathcal{O}(2^{d}/{\epsilon})$ iterations.

Not only does OMMWU match the dimension-dependece of MMWU, but it also
logarithmically outperform the dimension dependence of OMEG. Thus, between the
Frobenius norm (OMEG) and von Neumann entropy (OMMWU), the latter is the clear
choice for efficient implementation of the OMMP method for quantum zero-sum
games played in the spectraplex.

##  5 Conclusions

In this work, we proposed a “design” of classical algorithms for finding
$\epsilon$-Nash equibliria of finite-valued quantum zero-sum games, served to
unify and clarify the landscape of quantum zero-sum game algorithms. We
demonstrated that the best prior method—Jain and Watrous’ Matrix
Multiplicative Weight Updates (MMWU) algorithm [29]–fits within the framework
as an instantiation of our proposed Matrix Dual Averaging (MDA) method. To
improve upon MMWU, we leveraged a gradient-based view, which enabled us to
draw from the large classical game theory and optimization literature for
solving variational inequalities of monotone, Lipschitz operators. Via this
gradient-based perspective, proximal steps, and “optimism” we proposed and
proved convergence guarantees of the Optimistic Matrix Mirror Prox (OMMP)
method. When instantiated with the von Neumann entropy regularizer, OMMP
implements the Optimistic Matrix Multiplicative Weights Update (OMMWU)
algorithm, which we proved to have dimensionality dependence equivalent to
that of MMWU. However, OMMWU establishes a new benchmark in convergence rate,
achieving a quadratic speed-up relative MMWU, while still requiring only one
gradient call per iteration.

With these contributions several challenges remain and interesting open
problems are introduced:

  1. 1.

Tightness of Convergence Rates: Our first conjecture posits that the
$O(1/\epsilon)$ convergence rate of 1-MMP is tight, but empirical or
theoretical confirmation remains elusive. Establishing the tightness would
clarify the limitations of our method. See 1

  2. 2.

Dimensionality of Quantum Games: Our work focuses on finite quantum zero-sum
games. Extending the applicability of our algorithms to games in infinite
dimensional Hilbert spaces remains an open challenge that could have
theoretical applications to operator analysis and practical implications to
derive new optimal transport techniques.

  3. 3.

Robustness to Noise: Quantum systems are inherently prone to noise. Important
for practical implementations is how our method performs with different kind
of feedback under noisy measurement conditions.

  4. 4.

Computational Scaling Challenges: Our proposal requires knowledge of the
entire system density matrix in each iteration, which is computationally
prohibitive for large system sizes. Therefore, it would be valuable to explore
efficient classical representations of the system state or propose entirely
quantum protocols/algorithms which can eliminate the need to classically
learn/store the system state.

  5. 5.

Going Beyond 2-Player Zero-Sum: Early forays into quantum potential games [73]
and quantum zero-sum polymatrix games [45] are suggestive of exciting
possibilities for further research with novel classes of games and solution
concepts.

##  6 Acknowledgements

Part of this work was done while the authors were visiting the Simons
Institute for the Theory of Computing and the Institute for Mathematical
Sciences at the National University of Singapore. We thank Rahul Jain for
useful discussions. F.V. is supported by NSF grant DGE-2146752 and the
Vannevar Bush Faculty Fellowship Program grant N00014-21-1-2941. E.V. is
grateful for financial support by the Post-Doctoral FODSI-Simons Fellowship.
P.M. is a member of the Archimedes Unit—Athena RC—Department of Mathematics,
National & Kapodistrian University of Athens, and gratefully acknowledges
financial support by the French National Research Agency (ANR) in the
framework of the “Investissements d’avenir” program (ANR-15- IDEX-02), the
LabEx PERSYVAL (ANR-11-LABX-0025-01), MIAI@Grenoble Alpes (ANR-19- P3IA-0003),
and project MIS 5154714 of the National Recovery and Resilience Plan Greece
2.0 funded by the European Union under the NextGenerationEU Program. G.P.
would like to acknowledge National Research Foundation, Singapore and DSO
National Laboratories under its AI Singapore Program (AISG Award No:
AISG2-RP-2020-016), NRF 2018 Fellowship NRF-NRFF2018- 07, NRF2019-NRF-ANR095
ALIAS grant, grant PIESGP-AI-2020-01, AME Programmatic Fund (Grant
No.A20H6b0151) from the Agency for Science, Technology and Research (A*STAR)
and Provost’s Chair Professorship grant RGEPPV2101. M.J. is partially
supported by NSF Grant IIS-1901252 and Vannevar Bush Faculty Fellowship
program under grant number N00014-21-1-2941.

## References

  * [1] John Von Neumann.  “Zur Theorie der Gesellschaftsspiele”.  Mathematische Annalen 100, 295–320 (1928).
  * [2] John Nash.  “Non-Cooperative Games”.  Annals of Mathematics 54, 286–295 (1951).  url: http://www.jstor.org/stable/1969529.
  * [3] Parvathy S. Pillai and Shrisha Rao.  “Resource Allocation in Cloud Computing Using the Uncertainty Principle of Game Theory”.  IEEE Systems Journal 10, 637–648 (2016).
  * [4] Steven J. Brams.  “Game theory and politics”.  Free Press. (1975).  url: https://cir.nii.ac.jp/crid/1130000793730346368.
  * [5] David Silver, Aja Huang, Chris J. Maddison, Arthur Guez, Laurent Sifre, George van den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, Sander Dieleman, Dominik Grewe, John Nham, Nal Kalchbrenner, Ilya Sutskever, Timothy Lillicrap, Madeleine Leach, Koray Kavukcuoglu, Thore Graepel, and Demis Hassabis.  “Mastering the game of Go with deep neural networks and tree search”.  Nature 529, 484–489 (2016).
  * [6] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  “Generative Adversarial Nets”.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems.  Volume 27.  Curran Associates, Inc. (2014).  url: https://proceedings.neurips.cc/paper_files/paper/2014/file/5ca3e9b122f61f8f06494c97b1afccf3-Paper.pdf.
  * [7] Lerrel Pinto, James Davidson, Rahul Sukthankar, and Abhinav Gupta.  “Robust Adversarial Reinforcement Learning”.  In International Conference on Machine Learning.  Pages 2817–2826.  PMLR (2017).
  * [8] J. S. Bell.  “On the Einstein Podolsky Rosen paradox”.  Physics Physique Fizika 1, 195–200 (1964).
  * [9] John F. Clauser, Michael A. Horne, Abner Shimony, and Richard A. Holt.  “Proposed Experiment to Test Local Hidden-Variable Theories”.  Phys. Rev. Lett. 23, 880–884 (1969).
  * [10] N. David Mermin.  “Simple unified form for the major no-hidden-variables theorems”.  Phys. Rev. Lett. 65, 3373–3376 (1990).
  * [11] Lucien Hardy.  “Nonlocality for two particles without inequalities for almost all entangled states”.  Phys. Rev. Lett. 71, 1665–1668 (1993).
  * [12] P.K. Aravind.  “Bell’s theorem without inequalities and only two distant observers”.  In Optical Fiber Communication Conference and International Conference on Quantum Information.  Page EAPA7.  Optica Publishing Group (2001).
  * [13] Daniel M. Greenberger, Michael A. Horne, and Anton Zeilinger.  “Going beyond bell’s theorem” (2007).  arXiv:0712.0921.
  * [14] Ben W. Reichardt, Falk Unger, and Umesh Vazirani.  “Classical command of quantum systems”.  Nature 496, 456–460 (2013).
  * [15] Zhengfeng Ji, Anand Natarajan, Thomas Vidick, John Wright, and Henry Yuen.  “MIP* = RE”.  Commun. ACM 64, 131–138 (2021).
  * [16] Yael Kalai, Alex Lombardi, Vinod Vaikuntanathan, and Lisa Yang.  “Quantum Advantage from Any Non-local Game”.  In Barna Saha and Rocco A. Servedio, editors, Proceedings of the 55th Annual ACM Symposium on Theory of Computing, STOC 2023, Orlando, FL, USA, June 20-23, 2023.  Pages 1617–1628.  ACM (2023).
  * [17] G. Gutoski.  “Upper bounds for quantum interactive proofs with competing provers”.  In 20th Annual IEEE Conference on Computational Complexity (CCC’05).  Pages 334–343.  (2005).
  * [18] Gus Gutoski and John Watrous.  “Quantum interactive proofs with competing provers”.  In Volker Diekert and Bruno Durand, editors, STACS 2005.  Pages 605–616.  Berlin, Heidelberg (2005). Springer Berlin Heidelberg.
  * [19] Gus Gutoski and John Watrous.  “Toward a general theory of quantum games”.  In Proceedings of the Thirty-Ninth Annual ACM Symposium on Theory of Computing.  Page 565–574.  STOC ’07New York, NY, USA (2007). Association for Computing Machinery.
  * [20] Hirotada Kobayashi and Keiji Matsumoto.  “Quantum multi-prover interactive proof systems with limited prior entanglement”.  Journal of Computer and System Sciences 66, 429–450 (2003).
  * [21] R. Cleve, P. Hoyer, B. Toner, and J. Watrous.  “Consequences and limits of nonlocal strategies”.  In Proceedings. 19th IEEE Annual Conference on Computational Complexity, 2004.  Pages 236–249.  (2004).
  * [22] Julia Kempe, Hirotada Kobayashi, Keiji Matsumoto, and Thomas Vidick.  “Using Entanglement in Quantum Multi-Prover Interactive Proofs”.  computational complexity 18, 273–307 (2009).
  * [23] Julia Kempe, Hirotada Kobayashi, Keiji Matsumoto, Ben Toner, and Thomas Vidick.  “Entangled Games Are Hard to Approximate”.  SIAM Journal on Computing 40, 848–877 (2011).
  * [24] Andris Ambainis.  “A New Protocol and Lower Bounds for Quantum Coin Flipping”.  In Proceedings of the Thirty-Third Annual ACM Symposium on Theory of Computing.  Page 134–142.  STOC ’01New York, NY, USA (2001). Association for Computing Machinery.
  * [25] R. W. Spekkens and Terry Rudolph.  “Quantum Protocol for Cheat-Sensitive Weak Coin Flipping”.  Phys. Rev. Lett. 89, 227901 (2002).
  * [26] Carlos Mochon.  “Quantum weak coin flipping with arbitrarily small bias” (2007).  arXiv:0711.4114.
  * [27] Carl A. Miller.  “The Impossibility of Efficient Quantum Weak Coin Flipping”.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing.  Page 916–929.  STOC 2020New York, NY, USA (2020). Association for Computing Machinery.
  * [28] David A. Meyer.  “Quantum Strategies”.  Phys. Rev. Lett. 82, 1052–1055 (1999).
  * [29] Rahul Jain and John Watrous.  “Parallel Approximation of Non-interactive Zero-sum Quantum Games”.  In 2009 24th Annual IEEE Conference on Computational Complexity.  Pages 243–253.  (2009).
  * [30] Pierre-Luc Dallaire-Demers and Nathan Killoran.  “Quantum generative adversarial networks”.  Phys. Rev. A 98, 012324 (2018).
  * [31] Shouvanik Chakrabarti, Yiming Huang, Tongyang Li, Soheil Feizi, and Xiaodi Wu.  “Quantum Wasserstein Generative Adversarial Networks”.  In Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada.  Pages 6778–6789.  (2019).  url: https://proceedings.neurips.cc/paper/2019/hash/f35fd567065af297ae65b621e0a21ae9-Abstract.html.
  * [32] John Bostanci and John Watrous.  “Quantum game theory and the complexity of approximating quantum Nash equilibria”.  Quantum 6, 882 (2022).
  * [33] Kyriakos Lotidis, Panayotis Mertikopoulos, Nicholas Bambos, and Jose Blanchet.  “Payoff-based learning with matrix multiplicative weights in quantum games”.  In NeurIPS ’23: Proceedings of the 37th International Conference on Neural Information Processing Systems.  (2023).
  * [34] Koji Tsuda, Gunnar Rätsch, and Manfred K. Warmuth.  “Matrix Exponentiated Gradient Updates for On-line Learning and Bregman Projection”.  Journal of Machine Learning Research 6, 995–1018 (2005).  url: http://jmlr.org/papers/v6/tsuda05a.html.
  * [35] Arkadi Semen Nemirovski and David Berkovich Yudin.  “Problem complexity and method efficiency in optimization”.  Wiley. New York, NY (1983).
  * [36] Amir Beck and Marc Teboulle.  “Mirror descent and nonlinear projected subgradient methods for convex optimization”.  Operations Research Letters 31, 167–175 (2003).
  * [37] Rahul Jain, Zhengfeng Ji, Sarvagya Upadhyay, and John Watrous.  “QIP = PSPACE”.  J. ACM 58, 30:1–30:27 (2011).
  * [38] Sanjeev Arora and Satyen Kale.  “A Combinatorial, Primal-Dual Approach to Semidefinite Programs”.  In Proceedings of the Thirty-Ninth Annual ACM Symposium on Theory of Computing.  Page 227–236.  STOC ’07New York, NY, USA (2007). Association for Computing Machinery.
  * [39] Li Chen, Rasmus Kyng, Maximilian Probst Gutenberg, and Sushant Sachdeva.  “A Simple Framework for Finding Balanced Sparse Cuts via APSP”.  In Symposium on Simplicity in Algorithms (SOSA).  Pages 42–55.  SIAM (2023).
  * [40] Zeyuan Allen-Zhu, Zhenyu Liao, and Lorenzo Orecchia.  “Spectral Sparsification and Regret Minimization Beyond Matrix Multiplicative Updates”.  In Proceedings of the Forty-Seventh Annual ACM Symposium on Theory of Computing.  Page 237–245.  STOC ’15New York, NY, USA (2015). Association for Computing Machinery.
  * [41] Panayotis Mertikopoulos, E. Veronica Belmega, and Aris L. Moustakas.  “Matrix exponential learning: Distributed optimization in MIMO systems”.  In 2012 IEEE International Symposium on Information Theory Proceedings.  Pages 3028–3032.  (2012).
  * [42] Panayotis Mertikopoulos and Aris L. Moustakas.  “Learning in an Uncertain World: MIMO Covariance Matrix Optimization With Imperfect Feedback”.  IEEE Transactions on Signal Processing 64, 5–18 (2016).
  * [43] Olivier Bilenne, Panayotis Mertikopoulos, and Elena Veronica Belmega.  “Fast Optimization With Zeroth-Order Feedback in Distributed, Multi-User MIMO Systems”.  IEEE Transactions on Signal Processing 68, 6085–6100 (2020).
  * [44] Sham M Kakade, Shai Shalev-Shwartz, and Ambuj Tewari.  “Regularization techniques for learning with matrices”.  The Journal of Machine Learning Research 13, 1865–1890 (2012).
  * [45] Wayne Lin, Georgios Piliouras, Ryann Sim, and Antonios Varvitsiotis.  “No-regret learning and equilibrium computation in quantum games” (2023).  arXiv:2310.08473.
  * [46] Sanjeev Arora, Elad Hazan, and Satyen Kale.  “The Multiplicative Weights Update Method: a Meta-Algorithm and Applications”.  Theory of Computing 8, 121–164 (2012).
  * [47] Arkadi Nemirovski.  “Prox-Method with Rate of Convergence O(1/t) for Variational Inequalities with Lipschitz Continuous Monotone Operators and Smooth Convex-Concave Saddle Point Problems”.  SIAM Journal on Optimization 15, 229–251 (2004).
  * [48] Alfred Auslender and Marc Teboulle.  “Interior projection-like methods for monotone variational inequalities”.  Mathematical Programming 104, 39–68 (2005).
  * [49] Yuyuan Ouyang and Yangyang Xu.  “Lower complexity bounds of first-order methods for convex-concave bilinear saddle-point problems”.  Mathematical Programming 185, 1–35 (2021).
  * [50] Galina M Korpelevich.  “The extragradient method for finding saddle points and other problems”.  Matecon 12, 747–756 (1976).  url: https://cir.nii.ac.jp/crid/1571698600143951616.
  * [51] Boris Teodorovich Polyak and Anatoli Juditsky.  “Acceleration of stochastic approximation by averaging”.  SIAM Journal on Control and Optimization 30, 838–855 (1992).
  * [52] Leonid Denisovich Popov.  “A modification of the Arrow-Hurwicz method for search of saddle points”.  Mathematical notes of the Academy of Sciences of the USSR 28, 845–848 (1980).
  * [53] Sasha Rakhlin and Karthik Sridharan.  “Optimization, learning, and games with predictable sequences”.  Advances in Neural Information Processing Systems26 (2013).
  * [54] Kyriakos Lotidis, Panayotis Mertikopoulos, and Nicholas Bambos.  “Learning in quantum games” (2023).  arXiv:2302.02333.
  * [55] Ilayda Canyakmaz, Wayne Lin, Georgios Piliouras, and Antonios Varvitsiotis.  “Multiplicative updates for online convex optimization over symmetric cones” (2023).
  * [56] Alina Ene and Huy Lê Nguyên.  “Adaptive and Universal Algorithms for Variational Inequalities with Optimal Convergence”.  Proceedings of the AAAI Conference on Artificial Intelligence 36, 6559–6567 (2022).
  * [57] Gerard Debreu.  “A Social Equilibrium Existence Theorem”.  Proceedings of the National Academy of Sciences 38, 886–893 (1952).
  * [58] Gesualdo Scutari, Daniel P. Palomar, Francisco Facchinei, and Jong-shi Pang.  “Convex Optimization, Game Theory, and Variational Inequality Theory”.  IEEE Signal Processing Magazine 27, 35–49 (2010).
  * [59] Stephen Boyd and Lieven Vandenberghe.  “Convex Optimization”.  Cambridge University Press. (2004).
  * [60] Shai Shalev-Shwartz.  “Online Learning and Online Convex Optimization”.  Foundations and Trends in Machine Learning 4, 107–194 (2012).
  * [61] John Watrous.  “Theory of quantum information”.  2.3 Norms of operators, lecture notes, University of Waterloo (2011).
  * [62] Carl D. Meyer and Ian Stewart.  “Matrix Analysis and Applied Linear Algebra, Second Edition”.  Society for Industrial and Applied Mathematics. Philadelphia, PA (2023).
  * [63] Xingyu Zhou.  “On the Fenchel Duality between Strong Convexity and Lipschitz Continuous Gradient” (2018).  arXiv:1803.06573.
  * [64] D. Bertsekas, A. Nedic, and A. Ozdaglar.  “Convex Analysis and Optimization”.  Athena Scientific optimization and computation series. Athena Scientific. (2003).
  * [65] Constantin Zalinescu.  “Convex analysis in general vector spaces”.  World scientific. (2002).
  * [66] Panayotis Mertikopoulos, Bruno Lecouat, Houssam Zenati, Chuan-Sheng Foo, Vijay Chandrasekhar, and Georgios Piliouras.  “Optimistic mirror descent in saddle-point problems: Going the extra(-gradient) mile”.  In International Conference on Learning Representations.  (2019).  url: https://openreview.net/forum?id=Bkg8jjC9KQ.
  * [67] Sébastien Bubeck.  “Convex optimization: Algorithms and complexity”.  Found. Trends Mach. Learn. 8, 231–357 (2015).
  * [68] Yurii Nesterov.  “Primal-dual subgradient methods for convex problems”.  Mathematical Programming 120, 221–259 (2009).
  * [69] Andersen Ang.  “Projected gradient algorithm”.  https://angms.science/doc/CVX/CVX_PGD.pdf (2023).  Version: July 13, 2023. First draft: August 2, 2017.
  * [70] Satyen Kale.  “Efficient algorithms using the Multiplicative Weights Update method”.  PhD thesis.  Princeton University.  (2007).  url: https://www.proquest.com/dissertations-theses/efficient-algorithms-using-multiplicative-weights/docview/304824121/se-2.
  * [71] Anatoli Juditsky, Arkadi Nemirovski, and Claire Tauvel.  “Solving Variational Inequalities with Stochastic Mirror-Prox Algorithm”.  Stochastic Systems 1, 17–58 (2011).
  * [72] Yao-Liang Yu.  “The strong convexity of von neumann’s entropy”.  Lecture Notes for MSE 213: Optimization Theory (2020).  Available online: https://web.stanford.edu/~sidford/courses/20fa_opt_theory/sidford_mse213_2020fa_chap_5_extensions.pdf.
  * [73] Wayne Lin, Georgios Piliouras, Ryann Sim, and Antonios Varvitsiotis.  “Quantum Potential Games, Replicator Dynamics, and the Separability Problem” (2023).  arXiv:2302.04789.
  * [74] D. Petz.  “Bregman divergence as relative operator entropy”.  Acta Mathematica Hungarica 116, 127–131 (2007).
  * [75] Jean-Baptiste Hiriart-Urruty and Claude Lemaréchal.  “Fundamentals of convex analysis”.  Pages X, 259.  Grundlehren Text Editions. Springer Berlin, Heidelberg. (2001).  1 edition.

##  Appendix A Complex Matrix Differentiation and Gradients

In this appendix section, we review complex matrix differentiation and
gradients - an invaluable tool for framing our quantum game feedback and
analyzing our algorithm. Specifically, we will show why it is natural to
define the gradient of a real-valued function
$u:\mathbb{C}^{d}\rightarrow\mathbb{R}$ with respect to complex matrix
$Z\in\mathbb{C}^{d}$ as

$\nabla_{Z^{\top}}u(Z)$ (118)

and how this relates to a directional derivative. Concretely, our goal is to
compute the directional derivative of $u(Z(t))$ at $Z$ in the direction of
matrix $dZ/dt$.

In order to derive this derivative, we will start by introducing a few
important properties of $Z$ and $u$. To begin, note that $Z$ is a Hermitian
matrix of the form

$\displaystyle Z=X+iY,$ (119)

where $X$ is a real symmetric matrix and $Y$ is a real skew-symmetric matrix.
Thus, can use the chain rule to derive $du/dt$:

$\displaystyle\frac{du}{dt}=\sum_{i,j}\frac{\partial u}{\partial
X_{ij}}\frac{dX_{ij}}{dt}+\frac{\partial u}{\partial
Y_{ij}}\frac{dY_{ij}}{dt}.$ (120)

From the definition of the Wirtinger derivative, we also have that

$\displaystyle\frac{\partial u}{\partial Z_{ij}}=\frac{\partial u}{\partial
X_{ij}}-i\frac{\partial u}{\partial Y_{ij}}$ (121)

Finally, we note the following property of symmetric and skew-symmetric
matrices:

###### Proposition 6

For symmetric matrix $A$ and skew-symmetric matrix $B$,

$\displaystyle\textnormal{Tr}(AB)=\textnormal{Tr}(A^{\top}B)=\sum_{i,j}A_{ij}B_{ij}=0.$
(122)

#####  Proof.

$\textnormal{Tr}(AB)=\textnormal{Tr}(A^{\top}B)=\textnormal{Tr}(AB^{\top})=-\textnormal{Tr}(AB)\implies\textnormal{Tr}(AB)=0$
$\square$ Since $X$ is symmetric, $\partial u/\partial X$ and $dX/dt$ are
symmetric. Similarly, since $Y$ is skew-symmetric, $\partial u/\partial Y$ and
$dY/dt$ are symmetric. Therefore, by Proposition 6,

$\displaystyle\sum_{i,j}\frac{\partial u}{\partial X_{ij}}\frac{dY_{ij}}{dt}$
$\displaystyle=\textnormal{Tr}\left(\frac{\partial u}{\partial
X}\frac{dY}{dt}\right)=0$ (123) $\displaystyle\sum_{i,j}\frac{\partial
u}{\partial Y_{ij}}\frac{dX_{ij}}{dt}$
$\displaystyle=\textnormal{Tr}\left(\frac{\partial u}{\partial
Y}\frac{dX}{dt}\right)=0$ (124)

Now, we can combine the previous properties to obtain $du/dt$ expressed in
terms of the complex matrix $Z$. Starting from Equation 120, then leveraging
Equation 123 and Equation 124, we get that

$\displaystyle\frac{du}{dt}$ $\displaystyle=\sum_{i,j}\frac{\partial
u}{\partial X_{ij}}\frac{dX_{ij}}{dt}+\frac{\partial u}{\partial
Y_{ij}}\frac{dY_{ij}}{dt}$ (125) $\displaystyle=\sum_{i,j}\left(\frac{\partial
u}{\partial X_{ij}}\frac{dX_{ij}}{dt}+\frac{\partial u}{\partial
Y_{ij}}\frac{dY_{ij}}{dt}\right)+i\>\textnormal{Tr}\left(\frac{\partial
u}{\partial X}\frac{dY}{dt}\right)-i\>\textnormal{Tr}\left(\frac{\partial
u}{\partial Y}\frac{dX}{dt}\right)$ (126)
$\displaystyle=\sum_{i,j}\frac{\partial u}{\partial
X_{ij}}\frac{dX_{ij}}{dt}+i\frac{\partial u}{\partial
X_{ij}}\frac{dY_{ij}}{dt}-i\frac{\partial u}{\partial
Y_{ij}}\frac{dX_{ij}}{dt}+\frac{\partial u}{\partial
Y_{ij}}\frac{dY_{ij}}{dt}$ (127) $\displaystyle=\sum_{i,j}\left(\frac{\partial
u}{\partial X_{ij}}-i\frac{\partial u}{\partial
Y_{ij}}\right)\left(\frac{dX_{ij}}{dt}+i\frac{dY_{ij}}{dt}\right).$ (128)

By Equation 119 and Equation 121, this achieves the desired directional
derivative, which is expressed in terms of the gradient $\nabla_{Z^{\top}}u$:

$\displaystyle\frac{du}{dt}$ $\displaystyle=\sum_{i,j}\left(\frac{\partial
u}{\partial X_{ij}}-i\frac{\partial u}{\partial
Y_{ij}}\right)\left(\frac{dX_{ij}}{dt}+i\frac{dY_{ij}}{dt}\right)=\sum_{i,j}\frac{\partial
u}{\partial Z_{ij}}\frac{dZ_{ij}}{dt}$ (129)
$\displaystyle=\textnormal{Tr}\left(\left(\frac{\partial u}{\partial
Z}\right)^{\top}\frac{dZ}{dt}\right)=\textnormal{Tr}\left(\frac{\partial
u}{\partial
Z^{\top}}\frac{dZ}{dt}\right)=\textnormal{Tr}\left(\nabla_{Z^{\top}}u\cdot\frac{dZ}{dt}\right).$
(130)

##  Appendix B Mirror & Proximal Maps

In this appendix we offer further intuition for definitions and provide
intuitive proofs of key lemmas regarding the mirror and proximal maps/steps
introduced in Section 2.5.

See 6

#####  Proof.

Since $\mathcal{H}^{d}$ is the set of all complex matrices, this is an
unconstrained minimization problem. We can thus solve for the maximum by
setting the gradient of the minimization term to zero, as

$0=\nabla_{P^{\top}}\left(\textnormal{Tr}[D^{\dagger}P]-h(P)\right)=D-\partial
h(P),$ (131)

which implies the solution

$\partial h(P)=D\>\>\implies\>\>P=\partial^{-1}h(D).$ (132)

$\square$

See 4

#####  Proof.

The constrained minimization problem has Lagrangian

$\mathcal{L}=\textnormal{Tr}[\nabla
f(X_{t})^{\dagger}X_{t+1}]-\frac{1}{\eta}\left(h(X_{t+1})-h(X_{t})-\textnormal{Tr}[\nabla
h(X_{t})^{\dagger}(X_{t+1}-X_{t})]\right)-\lambda(\textnormal{Tr}[X_{t+1}-1]),$
(133)

with gradients

$\begin{cases}\nabla_{X_{t+1}^{\top}}\mathcal{L}&=\nabla
f(X_{t})-\frac{1}{\eta}\left(\nabla h(X_{t+1})-\nabla
h(X_{t})\right)-\lambda\mathbb{1}\\\
\nabla_{\lambda}\mathcal{L}&=\textnormal{Tr}[X_{t+1}]-1.\end{cases}.$ (134)

Setting both of these gradients to zero gives us the equations

$\begin{cases}\partial h(X_{t+1})&=\partial h(X_{t})+\eta\nabla
f(X_{t})-\eta\lambda\mathbb{1}\\\ \textnormal{Tr}[X_{t+1}]&=1,\end{cases},$
(135)

which implies the solution

$X_{t+1}=\partial^{-1}h(\partial h(X_{t})+\eta\nabla
f(X_{t})-\eta\lambda\mathbb{1}),$ (136)

with $\lambda$ chosen to ensure that $X_{t}$ lies in the feasible set
$\mathcal{X}^{d}$; i.e., it satisfies the constraint

$\textnormal{Tr}[X_{t+1}]=\textnormal{Tr}[\partial^{-1}h(\partial
h(X_{t})+\eta\nabla f(X_{t})-\eta\lambda\mathbb{1})]=1.$ (137)

$\square$

See 5

#####  Proof.

The constrained minimization problem has Lagrangian

$\displaystyle\mathcal{L}$ $\displaystyle=\textnormal{Tr}[\nabla
f(Y_{t+1})^{\dagger}(X_{t+1}-X_{t})]-\frac{1}{\eta}\mathcal{BD}_{h}(X_{t+1}||X)-\lambda(\textnormal{Tr}[X_{t+1}]-1),$
(138)

where

$\displaystyle\mathcal{BD}_{h}(X_{t+1}||X)=h(X_{t+1})-h(X_{t})-\textnormal{Tr}[\partial
h(X_{t})^{\dagger}(X_{t+1}-X_{t})].$ (139)

The Lagrangian has gradients

$\begin{cases}\nabla_{X_{t+1}^{\top}}\mathcal{L}&=\nabla
f(Y_{t+1})-\frac{1}{\eta}\left(\partial h(X_{t+1})-\partial
h(X_{t})\right)-\lambda\mathbb{1}\\\
\nabla_{\lambda}\mathcal{L}&=\textnormal{Tr}[X_{t+1}]-1.\end{cases}$ (140)

Setting both of these gradients to zero gives us the equations

$\begin{cases}\partial h(X_{t+1})&=\partial h(X_{t})+\eta\nabla
f(Y_{t+1})-\eta\lambda\mathbb{1}\\\ \textnormal{Tr}[X_{t+1}]&=1,\end{cases}$
(141)

which implies the solution

$X_{t+1}=\partial^{-1}h(\partial h(X_{t})+\eta\nabla
f(Y_{t+1})-\eta\lambda\mathbb{1}),$ (142)

with $\lambda$ chosen to ensure that $X_{t}$ lies in the feasible set
$\mathcal{X}^{d}$, i.e., it satisfies the constraint

$\textnormal{Tr}[X_{t}]=\textnormal{Tr}[\partial^{-1}h(\partial
h(X_{t})+\eta\nabla f(Y_{t+1})-\eta\lambda\mathbb{1})]=1.$ (143)

$\square$

See 6

#####  Proof.

We begin by noting that

$\displaystyle\nabla h(Y)=\frac{1}{2}\nabla||Y||_{F}^{2}=Y,$ (144)

which implies that the mirror map $\nabla h$ is the identity in this case.
Therefore,

$\displaystyle\nabla^{-1}h(Y)=\nabla h(Y)=Y.$ (145)

Furthermore, since

$\displaystyle||X-Y||_{F}^{2}$
$\displaystyle=||X||_{F}^{2}+||Y||_{F}^{2}-2\>\textnormal{Tr}[X^{\dagger}Y]$
(146)
$\displaystyle=||X||_{F}^{2}+(2||Y||_{F}^{2}-||Y||_{F}^{2})-2\>\textnormal{Tr}[X^{\dagger}Y]$
(147)
$\displaystyle=||X||_{F}^{2}-||Y||_{F}^{2}+2\>\textnormal{Tr}[Y^{\dagger}Y]-2\>\textnormal{Tr}[X^{\dagger}Y]$
(148)
$\displaystyle=||X||_{F}^{2}-||Y||_{F}^{2}-2\>\textnormal{Tr}[(X-Y)^{\dagger}Y],$
(149)

the Bregman divergence, as desired, is

$\displaystyle\mathcal{BD}_{h}(X||Y)$
$\displaystyle=h(X)-h(Y)-\textnormal{Tr}[(X-Y)^{\dagger}\nabla h(Y)]$ (150)
$\displaystyle=\frac{1}{2}||X||_{F}^{2}-\frac{1}{2}||Y||_{F}^{2}-\textnormal{Tr}[(X-Y)^{\dagger}\nabla
h(Y)]$ (151)
$\displaystyle=\frac{1}{2}||X||_{F}^{2}-\frac{1}{2}||Y||_{F}^{2}-\textnormal{Tr}[(X-Y)^{\dagger}Y]$
(152)
$\displaystyle=\frac{1}{2}\left(||X||_{F}^{2}-||Y||_{F}^{2}-2\textnormal{Tr}[(X-Y)^{\dagger}Y]\right)$
(153) $\displaystyle=\frac{1}{2}||X-Y||_{F}^{2}.$ (154)

The corresponding regularized best-response is

$\displaystyle\textnormal{Mir}\Pi_{\mathcal{C}}^{h}(Y)$
$\displaystyle=\underset{C\in\mathcal{C}}{\arg\min}\left\\{\textnormal{Tr}[Y^{\dagger}C]-\frac{1}{2}||C||_{F}^{2}\right\\},$

which has Lagrangian

$\displaystyle\mathcal{L}=\textnormal{Tr}[Y^{\dagger}C]-\frac{1}{2}||C||_{F}^{2}-\lambda(\textnormal{Tr}[C]-1),$
(155)

with gradients

$\displaystyle\nabla\mathcal{L}=\begin{cases}\nabla_{C^{\top}}\mathcal{L}&=Y-C-\lambda\mathbb{1}\\\
\nabla_{\lambda}\mathcal{L}&=-\textnormal{Tr}[C]+1.\end{cases}$ (156)

Note, however, that these are also gradients of the Lagrangian

$\displaystyle\mathcal{L}^{\prime}=\frac{1}{2}||Y-C||_{F}^{2}-\lambda(\textnormal{Tr}[C]-1),$
(157)

corresponding to the minimization of the orthogonal projection

$\displaystyle\textrm{Orth}\Pi_{\mathcal{C}}(Y)={\arg\min}_{C\in\mathcal{C}}||Y-C||_{F}^{2}.$

Therefore, as desired,

$\displaystyle\textnormal{Mir}\Pi_{\mathcal{C}}^{h}(Y)={\arg\min}_{C\in\mathcal{C}}||Y-C||_{F}^{2}=\textrm{Orth}\Pi_{\mathcal{C}}(Y).$

Finally, the corresponding proximal projection is given by

$\displaystyle\textnormal{Prox}\Pi^{h}_{\mathcal{C}}(X,Y)={\arg\min}_{C\in\mathcal{C}}\left\\{\textnormal{Tr}[Y^{\dagger}(X-C)]+\frac{1}{2}||C-X||^{2}_{F}\right\\},$
(158)

which has Lagrangian

$\displaystyle\mathcal{L}=\textnormal{Tr}[Y^{\dagger}(X-C)]+\frac{1}{2}||C-X||^{2}_{F}-\lambda(\textnormal{Tr}[C]-1),$
(159)

with gradients

$\displaystyle\nabla\mathcal{L}=\begin{cases}\nabla_{C^{\top}}\mathcal{L}&=X+Y-C-\lambda\mathbb{1}\\\
\nabla_{\lambda}\mathcal{L}&=-\textnormal{Tr}[C]+1.\end{cases}$ (160)

Note, however, that these are also gradients of the Lagrangian

$\displaystyle\mathcal{L}^{\prime}=\frac{1}{2}||X+Y-C||_{F}^{2}-\lambda(\textnormal{Tr}[C]-1),$
(161)

corresponding to the minimization of the orthogonal projection

$\displaystyle\textrm{Orth}\Pi_{\mathcal{C}}(X+Y)={\arg\min}_{C\in\mathcal{C}}||X+Y-C||_{F}^{2}.$

Therefore, as desired,

$\displaystyle\textnormal{Prox}\Pi^{h}_{\mathcal{C}}(X,Y)={\arg\min}_{C\in\mathcal{C}}||X+Y-C||_{F}^{2}=\textrm{Orth}\Pi_{\mathcal{C}}(X+Y).$

$\square$

See 7

#####  Proof.

We will begin with an informal proof of the Bregman divergence of the von
Neumann entropy of density matrices (for the original, rigorous proof refer to
[74]). In doing so, we will leverage the following facts about density
matrices $X$ and $Y$:

  1. 1.

$\nabla_{X^{\top}}\textnormal{Tr}[X\log X]=\mathbb{1}+\log X$

  2. 2.

$\textnormal{Tr}[X]=\textnormal{Tr}[Y]=1$

  3. 3.

Density matrices are Hermitian, which implies $X^{\dagger}=X$ and
$Y^{\dagger}=Y$.

Thus, the Bregman divergence is the relative quantum entropy as follows,

$\displaystyle\mathcal{BD}_{h}(X||Y)$