\subcaptionsetup

[figure]skip=-6mm,singlelinecheck=off \AtAppendix \AtAppendix \AtAppendix \AtAppendix

Quantum algorithm for large-scale market equilibrium computation

Po-Wei Huang¹ and Patrick Rebentrost^{1, 2}
¹Centre for Quantum Technologies, National University of Singapore
²Department of Computer Science, National University of Singapore
[email protected], [email protected]

Abstract

Classical algorithms for market equilibrium computation such as proportional response dynamics face scalability issues with Internet-based applications such as auctions, recommender systems, and fair division, despite having an almost linear runtime in terms of the product of buyers and goods. In this work, we provide the first quantum algorithm for market equilibrium computation with sub-linear performance. Our algorithm provides a polynomial runtime speedup in terms of the product of the number of buyers and goods while reaching the same optimization objective value as the classical algorithm. Numerical simulations of a system with 16384 buyers and goods support our theoretical results that our quantum algorithm provides a significant speedup.

1 Introduction

The balance of supply and demand is a fundamental and well-known law that determines the price of goods in a market. In a market with a set of $n$ buyers and $m$ goods, the competitive equilibrium [1, 2] determines the optimal price and allocation of goods such that the supply equals the demand in the given market. The computation of the competitive equilibrium is known as the market equilibrium computation problem, whose unique solution was shown to exist under a general model of the economics in the seminal work of Arrow and Debreu [3]. The relevance of such problems in algorithmic game theory [4, 5] is substantiated by the first welfare theorem, which implies that the competitive equilibria are Pareto-efficient [6], where no allocation is available that makes one agent better without making another one worse. In competitive equilibrium from equal income (CEEI) scenarios, such equilibria are further known to by envy-free [7, 8], where no agent would prefer an allocation received by another agent over their own.

The market equilibrium computation problem has, in recent years, been extended to various large-scale Internet-based markets [9], including auction markets [10], fair item allocation/fair division [11, 12, 13], scheduling problems [14] and recommender systems [15]. Such developments call for the need to further develop algorithmic theories for markets and the computation of market equilibria.

We focus on a particular type of market known as the Fisher market [16, 17], where there is a set of $n$ buyers interested in buying $m$ infinitely-divisible goods, and where each buyer has their monetary budget that has no intrinsic value apart from being used to purchase goods. We mainly consider Fisher markets with linear utilities, where the total utility gained by purchasing goods is strictly linear to the value and proportion of the goods obtained.

While combinatorial algorithms that can obtain exact and approximate solutions to such solutions have been discovered [18, 19, 20, 21], these algorithms tend to scale poorly against the growing number of buyers and goods. One can otherwise formulate the market equilibrium computation problem as an optimization problem that maximizes a convex objective function known as the Eisenberg-Gale (EG) convex program [22, 23]. Such optimization algorithms can produce approximate solutions much faster than that of combinatorial algorithms. One such example that is commonly used for the market equilibrium problem is the proportional response (PR) dynamics [24, 25]. The PR dynamics is an iterative algorithm that converges with a rate of $\nicefrac{{1}}{{T}}$ where $T$ is the number of iterations. Each iteration of the PR dynamics has a cost of $\mathcal{O}(mn)$ from proportionally updating individual bids that a buyer should make for different goods.

Given the high number of buyers and goods that can exist in Internet-based markets, the problem of further algorithmic speedups to the computation continues to be an active field of research. Gao and Kroer [26] discovered that by using projected gradient descent instead of PR dynamics, the market equilibrium can be found with linear convergence. Apart from the number of iterations, attempts to reduce the cost per iteration, such as using clustering to reduce the problem size [15], have also been made. However, it is not yet clear whether these methods can provide advantages beyond a constant-factor speedup.

In this work, we consider a Fisher market with $n$ buyers and $m$ goods, where the objective is to find an approximate market equilibrium whose EG objective function is within an additive error $\epsilon$ of the optimal EG objective value. We provide a method to reduce the cost per iteration by utilizing quantum norm estimation and quantum inner product estimation [27, 28] and provide the first quantum algorithm to achieve sublinear performance in terms of the product of the buyers and goods in market equilibrium computation. To arrive at the quantum algorithm, we show an alternate version of the PR dynamics with erroneous updates, which we term the faulty proportional response (FPR) dynamics. We then provide a quantum algorithm that provides a quadratic speedup in terms of the smaller dimension between buyers and goods, as well as less memory consumption, albeit being QRAM instead of classical RAM. We summarize our results in Table 1.

Algorithm	Iterations	Runtime	Memory	Result Prep.
PR dynamics [24]	$\displaystyle\frac{\log m}{\varepsilon}$	$\displaystyle\tilde{\mathcal{O}}\left(\frac{mn}{\varepsilon}\right)$	$\mathcal{O}(mn)$	NA, in RAM
Quantum alg.	$\displaystyle\frac{2\log m}{\varepsilon}$	$\displaystyle\tilde{\mathcal{O}}\left(\frac{\sqrt{mn\max(m,n)}}{\varepsilon^{2% }}\right)$	$\mathcal{O}(m+n)^{*}$	$\mathcal{O}(\operatorname{poly}\log mn)$

Table 1: Main result.

n

is the number of buyers,

m

is the number of goods, and

\varepsilon

indicates the additive error of the computed values with minimal-achievable EG objective value. The memory complexity for the quantum algorithm (annotated with *) refers to the use of quantum access to classical memory, which QRAM can achieve (see Definition 1), instead of classical RAM. As the computed competitive equilibrium takes

\mathcal{O}(mn)

memory to store, our quantum algorithm does not provide the entire matrix, but instead provide query access to the result. The result preparation column refers to the runtime cost of preparing query access to the competitive equilibrium.

2 Preliminaries

Notations.

Let $[n]:=\{0,1,\dots,n-1\}$ . We use $\odot$ to represent element-wise multiplication, as well as $\oplus$ for bit-wise XOR operation and $\otimes$ for tensor products. For vectors $u\in\mathbb{R}^{N}$ , we denote a vector’s $\ell_{p}$ norm by $\|v\|_{p}:=\sqrt[p]{\sum_{i=1}^{N}|v_{i}|^{p}}$ . Let $\mathcal{M}_{M\times N}(\mathbb{R})$ indicate the space of square matrices of size $M\times N$ over $\mathbb{R}$ . We denote the $i$ -th row vector of $A$ by $A_{i,*}$ and the $j$ -th column vector of $A$ by $A_{*,j}$ . We further define $\mathbb{I}$ as $[0,1]$ , and the $n$ -unit simplex as $\mathbb{S}^{n}$ , i.e. $\mathbb{S}^{n}=\{v\in\mathbb{I}^{n},\|v\|_{1}=1\}$ . For sets of numbers, we add the subscript $\cdot_{+}$ to indicate a constraint on positivity for elements in the set. We use $\ket{k}$ to denote a binary encoding of a real number $k$ up to arbitrary precision into a quantum state, and $\ket{\bar{0}}$ to denote a multi-qubit zero state whose number of qubits can be inferred from the context. Lastly, we use $\mathcal{\tilde{O}}(\cdot)$ to omit polylogarithmic factors in asymptotic runtime/memory analysis.

Quantum computation.

Quantum algorithms are shown to be able to provide asymptotic speedups over classical counterparts [29, 30, 31] by utilizing characteristics of quantum mechanics such as superposition to access data all at once. In this work, the main quantum speedup stems from the fast computation of $\ell_{1}$ norms and inner products [27, 28], which is in turn powered by a technique known as quantum amplitude estimation (QAE) [32]. Classical approximation algorithms that use Monte Carlo methods for probability estimation up to precision $\epsilon$ have runtime $\mathcal{O}(1/\epsilon^{2})$ due to the concentration of precision being correlated to the variance. On the other hand, QAE can achieve the same precision by reading input at once in superposition and repeatedly amplifying the precision of our estimation, which takes $\mathcal{O}(1/\epsilon)$ runtime and provides a quadratic speedup. Many subtle improvements to the QAE algorithm have since been made after its discovery, such as simplifying subroutines [33, 34, 35], restoring the initial state [36, 37] and compensating for bias [38].

Theorem 2.1 (Quantum amplitude estimation; Theorem 2, [39]).

Let $t\in\mathbb{N}$ . We are given one copy of a quantum state $\ket{\psi}$ as input, as well as a unitary transformation $U=I-2\lvert\psi\rangle\langle\psi\rvert$ , and a unitary transformation $V=I-2P$ for some projector $P$ . There exists a quantum algorithm that outputs $\tilde{a}$ , an estimate of $a=\|P\ket{\psi}\|^{2}$ , such that

|\tilde{a}-a|\leq 2\pi\frac{\sqrt{a(1-a)}}{M}+\frac{\pi^{2}}{M^{2}}

with probability at least $8/\pi^{2}$ , using $M$ applications of $U$ and $V$ each.

In this paper, we use QAE to estimate $\ell_{1}$ norms and inner products of vectors $v\in\mathbb{I}^{N}$ up to a multiplicative error in $M\in\mathcal{O}(\frac{\sqrt{N}}{\epsilon}\ln(\frac{1}{\delta}))$ runtime with probability $1-\delta$ , invoking a quadratic speedup in both the dimension and the error rate. We defer the formulation and details to Appendix A.

Apart from quantum subroutines that provide speedups, we also require the usage of arithmetic operations such as addition, subtraction, multiplication, and division on quantum computers. We assume the arithmetic model, which would allow us to ignore issues arising from the fixed point representation of numbers¹¹1If the fixed point representation with an additive error of $\mu$ is considered, the additional multiplicative cost required for operations is then $\mathcal{O}(\operatorname{poly}\log\nicefrac{{1}}{{\mu}})$ . Considering $\mu\in\Omega(1/\operatorname{poly}(m,n))$ , the additional cost is $\mathcal{O}(\operatorname{poly}\log(m,n))$ , which are polylogarithmic factors that we already omit in this paper.. We further assume that we have access to quantum arithmetic circuits [40, 41] that can perform such arithmetic operations in $\mathcal{O}(1)$ gates, and that by using such circuits, computation of the $n$ -th power of a number, where $n\in\mathbb{N}$ , can be achieved in $\mathcal{O}\left(\operatorname{poly}\log n\right)$ gates, using methods like binary exponentiation [42]. We note that quantum arithmetic circuits can be used to execute the same operation on multiple numbers in parallel if the numbers are held in superposition.

Lastly, we need to access the input matrices and intermediate vectors as a superposition of encoded quantum states. Such quantum access to the classical data in memory can be achieved by quantum random access memory (QRAM)²²2Our memory unit can be more precisely termed QRACM [43, 44] or QROM [45] as opposed to QRAQM [43, 44] or QRAG [46], whose memory registers store quantum states instead of classical numbers. However, both are more commonly and jointly referred to as QRAM in literature. as follows. We refer the reader to [44] for a survey on QRAM.

Definition 1 (Quantum random access memory; [47, 48]).

Let $n\in\mathbb{N}$ and $c\in\mathcal{O}(1)$ . Also let $w$ be a vector of bit strings such that $\forall i\in[n],w_{i}\in\{0,1\}^{c}$ . A quantum RAM provides access to $w_{i}$ in superposition after a one-time construction cost of $\tilde{\mathcal{O}}\left(n\right)$ , where each access costs $\mathcal{O}\left(\operatorname{poly}\log n\right)$ .

Fisher market equilibrium.

In the Fisher market model [16, 17], we are given a market of $m$ infinitely divisible goods to be divided among $n$ buyers. Without loss of generality, we assume a unit supply for each good. Each buyer $i\in[n]$ has a budget of $B_{i}>0$ that has no intrinsic value apart from being used to purchase goods where, again without loss of generality, we assume $B\in\mathbb{S}^{n}$ . Each buyer also has a utility function $u_{i}:\mathbb{R}^{m}\to\mathbb{R}_{+}$ that maps an allocation of portions of $m$ items to a utility value. We can then define the allocation matrix $x\in\mathcal{M}_{n\times m}(\mathbb{R}_{+})$ such that $x_{ij}$ is the portion of item $j$ allocated to buyer $i$ , where $x_{i}\in\mathbb{R}^{m}$ is the bundle of products allocated to buyer $i$ . In this paper, we consider linear utility functions such that $u_{i}(x_{i})=\sum_{j\in[m]}v_{ij}x_{ij}$ , where $v_{ij}>0$ is the value for a unit of item $j$ for buyer $i$ .

Given the Fisher market, we want to compute its competitive equilibrium, which consists of the price vector $p\in\mathbb{R}^{m}$ for each item $j$ and allocation matrix $x$ such that each buyer $i$ exhausts their entire budget $B_{i}$ to acquire a bundle of items $x_{i}$ that maximizes each of their utility $u_{i}(x_{i})$ .

The market equilibrium of Fisher markets can be captured by solving the Eisenberg-Gale (EG) convex program [22, 23]. The program is derived from maximizing the budget-weighted geometric mean of the buyers’ utilities, which satisfies natural properties such as invariance of the optimal solution to rescaling and splitting [49]. Using the $\log$ on the geometric mean, the EG program is as follows:

\max_{x\geq 0}\sum_{i\in[n]}B_{i}\log u_{i}(x_{i})\text{ s.t. }\sum_{i\in[n]}x% _{ij}\leq 1,\forall j\in[m].

(2.1)

Such convex programs (maximization of a concave function subject to constraints) can be solved by interior point methods [50], but may not scale to large markets. We discuss this further in Section 6.

For the linear Fisher market, an alternative convex program that obtains the same market equilibrium was shown by Shmyrev [51]. Supposing that each buyer $i$ submits a bid $b_{ij}$ for item $j$ such that the sum of the bid of the buyer matches their budget $B_{i}$ such that each buyer $i$ is allocated $x_{ij}=b_{ij}/p_{j}$ of item $j$ , we have the following convex program:

\max_{b\geq 0}\sum_{ij}b_{ij}\log\frac{v_{ij}}{p_{j}}\text{ s.t. }\sum_{i\in[n% ]}b_{ij}=p_{j},\forall j\in[m];\sum_{j\in[m]}b_{ij}=B_{i},\forall i\in[n].

(2.2)

As the allocation matrix and price vector can be directly computed from and be used to compute the bid matrix, the bid matrix can be used as a direct representation of the market equilibrium itself, and hence, is the output of the algorithms we discuss in our paper.

Proportional response dynamics.

The proportional response (PR) dynamics is an iterative algorithm [24, 25, 52] that obtains the Fisher market equilibrium computation by updating the bids $b_{ij}$ submitted by buyer $i$ for item $j$ . For each time step, the elements of the price vector $p_{j}$ are computed by summing the bids for item $j$ such that $p_{j}=\sum_{i}b_{ij}$ . The allocation $x_{ij}$ is then obtained by taking $x_{ij}=b_{ij}/p_{j}$ . The buyers then update the bids such that the new bid is proportional to the utility $u_{i}=\sum_{j}v_{ij}x_{ij}$ gained in the current time step such that $b_{ij}^{(t+1)}=B_{i}v_{ij}x_{ij}^{(t)}/u_{i}^{(t)}$ . It was shown by Birnbaum et al. [53] that the PR dynamics is equivalent to mirror descent [54, 55] to a Bregman divergence [56] of the Shmyrev convex program.

For ease of discussion, we write the objective function of the EG and Shmyrev convex programs as functions of the bid matrix $b$ , obtaining the EG objective function $\Phi(b)=-\sum_{i\in[n]}B_{i}\log u_{i}$ and Shmyrev objective function $\Psi(b)=\sum_{i\in[n],j\in[m]}b_{ij}\log\frac{p_{j}}{v_{ij}}$ . We denote the optimal bid $b^{*}=\operatorname{arg\,min}_{b\in\mathcal{S}}\Phi(b^{*})$ , where $\mathcal{S}=\left\{b\in\mathcal{M}_{n\times m}(\mathbb{I}):\sum_{j\in[m]}b_{i,% *}=B_{i}\right\}$ .

The convergence bounds of the PR dynamics regarding the EG and Shmyrev objective functions for linear Fisher markets were found as follows:

Theorem 2.2 (Convergence of PR dynamics; [53]).

Considering a linear Fisher market, for $b_{ij}^{(t)}$ as iteratively defined by the proportional response dynamics where $b_{ij}^{(0)}=\frac{B_{i}}{m}$ , we have

\Psi(b^{(T)})-\Psi(b^{*})\leq\frac{\log m}{T},\quad\Phi(b^{(T-1)})-\Phi(b^{*})% \leq\frac{\log m}{T}.

(2.3)

An alternate end-to-end proof of the convergence of both convex programs that varies from Birnbaum et al. [53]’s approach and centered around the EG function can be found in Appendix B, elements of which we use in the proof of later sections. Two notable results that we prove and utilize are: 1) $\Psi(b^{(t+1)})\leq\Phi(b^{(t)})+\sum_{i\in[n]}B_{i}\log B_{i}\leq\Psi(b^{(t)})$ , and 2) the telescoping sum of the difference of the KL divergence of the optimal bid and the iterating bids can be lower bounded by the difference of the current EG objective function and the optimal EG function.

3 Faulty proportional response dynamics

Before moving on to our quantum algorithm, we propose the faulty proportional response (FPR) dynamics, which computes an erroneous update to compute a sequence of bids $\hat{b}_{ij}^{(t)}$ , which still retains a convergence guarantee, serving as a counterpart to Theorem 2.2. We first define a faulty update we use for the FPR dynamics:

Definition 2 (Faulty proportional response update).

Let $t\geq 0$ and $\hat{b}^{(t)}\in\mathcal{M}_{n\times m}(\mathbb{R}_{+})$ . Given $\epsilon_{p}\in(0,0.5)$ such that $\forall j,t,\|\tilde{p}_{j}^{(t)}-\hat{p}_{j}^{(t)}\|\leq\hat{p}_{j}^{(t)}% \epsilon_{p}$ where $\hat{p}_{j}^{(t)}=\sum_{i\in[n]}\hat{b}_{ij}^{(t)}$ . Further, given $\epsilon_{\nu}\in(0,0.5)$ such that $\forall i,t,\|\tilde{\nu}_{i}^{(t)}-\hat{\nu}_{i}^{(t)}\|\leq\hat{\nu}_{j}^{(t% )}\epsilon_{\nu}$ where $\hat{\nu}_{i}^{(t)}=\sum_{j\in[m]}v_{ij}\hat{b}_{ij}^{(t)}/\tilde{p}_{j}^{(t)}$ . A faulty proportional response update of the bids from timestep $t$ to $t+1$ is then expressed as follows:

\hat{x}_{ij}^{(t)}=\frac{\hat{b}_{ij}^{(t)}}{\tilde{p}_{j}^{(t)}},\quad\hat{b}% _{ij}^{(t+1)}=B_{i}\frac{v_{ij}x_{ij}^{(t)}}{\tilde{\nu}_{i}^{(t)}}.

Note that $\tilde{p}_{j}$ provides an estimation to the price $\hat{p}_{j}=\sum_{i\in[n]}\hat{b}_{ij}$ , $\tilde{\nu}_{i}$ does not provide an estimation to the exact utility $\hat{u}_{i}=\sum_{j\in[m]}v_{ij}\hat{b}_{ij}/\hat{p}_{j}$ , but $\nu_{i}$ , which replaces $\hat{p}_{j}$ in the computation of $u_{i}$ with $\tilde{p}_{j}$ .

The convergence bounds of the FPR dynamics regarding the EG objective function for linear Fisher markets were found as follows:

Theorem 3.1 (Convergence of the FPR dynamics).

Considering a linear Fisher market, for $b_{ij}^{(t)}$ as iteratively defined by the faulty proportional response dynamics where $\hat{b}_{ij}^{(0)}=\frac{B_{i}}{m}$ , we have

\min_{t\in[T]}\Phi(\hat{b}^{(t)})-\Phi(b^{*})\leq\frac{2\log m}{T}

when $\epsilon_{\nu}\leq\frac{\log m}{8T}$ and $\epsilon_{p}\leq\frac{\log m}{6T}$ .

A high-level idea of the proof follows from the telescoping sum trick to upper bound the EG objective functions with KL divergence from our proof of PR dynamics but with the consideration of error. We show an end-to-end proof of the convergence of the EG objective function in Appendix C.

Notice that in the FPR dynamics, we do not enforce the monotonicity of the iterations, but instead simply take the minimum value over all iterations. The error terms $\epsilon_{p}$ and $\epsilon_{\nu}$ in the FPR dynamics are only upper bounded such that the total sum of objective values over $T$ iterations (plus the original iteration) can be upper bounded by $\log m$ plus an accumulated error over $T$ iterations also within $\log m$ . If we enforce the monotonicity of the iterations to take the last iteration, the error would require $\mathcal{O}(1/T^{2})$ precision and would produce a $\mathcal{O}(T^{3})$ algorithmic dependency instead of $\mathcal{O}(T^{2})$ .

However, given the formulation of a faulty update, a problem that comes into question is whether the computation of the exact value of the function $\Phi(b)$ is supported, as we do not compute $u_{i}$ in the process of updating. Without computation of $\Phi(b)$ , one can not be sure which iteration of $\hat{b}^{(t)}$ is the minimum. However, we use the computed value of $\tilde{\nu}_{i}^{(t)}$ as an estimator for the function $\Phi(b)$ . The following result is then obtained.

Theorem 3.2.

Considering a linear Fisher market, for $b_{ij}^{(t)}$ as iteratively defined by the faulty proportional response dynamics where $\hat{b}_{ij}^{(0)}=\frac{B_{i}}{m}$ . Let $t^{*}=\operatorname{arg\,max}_{t\in[T]}B_{i}\log\tilde{\nu}_{i}^{(t)}$ . Then

\Phi(\hat{b}^{(t^{*})})-\Phi(b^{*})\leq\frac{2\log m}{T}

when $\epsilon_{\nu}\leq\frac{\log m}{8T}$ and $\epsilon_{p}\leq\frac{\log m}{6T}$ .

The proof of this theorem can similarly be found in Appendix C, which has the same proof idea as Theorem 3.1 apart from some slight differences in error handling.

4 Quantum algorithm

We present our quantum algorithm for solving linear Fisher market equilibrium computation based on the FPR dynamics. Our quantum algorithm does not aim to provide speedups in terms of the number of iterations but provides speedups on the iteration cost of the PR dynamics algorithm. Our algorithm, while reducing the runtime in terms of the number of buyers $n$ or goods $m$ , increases runtime in terms of the number of iterations $T$ but as the $T$ is logarithmically dependent on $m$ , there is an overall quadratic speedup provided in the smaller of the two dimensions.

In this section, we further assume that $v\in\mathcal{M}_{n\times m}(\mathbb{I})$ . We note that the multiplicative scaling of $v_{ij}$ does not affect the bid matrix $b$ generated in the FPR dynamics as errors are multiplicative. Hence if the values are larger than $1$ , we scale down the values by dividing the queried $v_{ij}$ by a number that is larger than $\max_{ij}v_{ij}$ .

To compute the market equilibrium for the Fisher market by the FPR dynamics in the quantum setting, we require the data input of both the budget vector $B\in\mathbb{S}^{n}$ and the value matrix $v\in\mathcal{M}_{n\times m}(\mathbb{I})$ . We assume quantum query access to the budget and vector and value matrix by the index is readily given to us as part of the problem input without having to load classical data into a quantum system. That is, given an index state and ancilla quantum registers we can store the value of the budget and value according to the index in the ancilla register. Note that these operations can be performed in superposition, such that $\frac{1}{\sqrt{mn}}\sum_{i,j}\boldsymbol{\ket{i}}\boldsymbol{\ket{j}}\ket{\bar% {0}}\rightarrow\frac{1}{\sqrt{mn}}\sum_{i,j}\boldsymbol{\ket{i}}\boldsymbol{% \ket{j}}\ket{v_{ij}}$ .

We do not explicitly state how the data input of the budget and value entries are generated; they could be extracted from entries of a matrix already preloaded in QRAM, or generated/reconstructed from a low-rank approximation of the matrix [15], which takes $\mathcal{O}\left(k\operatorname{poly}\log mn\right)$ cost to access $k$ -rank approximations using quantum arithmetic circuits, but with much lower memory consumption³³3With low-rank approximations, the loading of classical data into QRAM would only take $\tilde{\mathcal{O}}(k(m+n))$ runtime.. We note that the low-rank approximation assumption of the value matrix has not yet been utilized to produce reductions in resource consumption in classical methods as the PR dynamics and other methods to compute market equilibrium [26] require all $\mathcal{O}(mn)$ entries of the full value matrix.

Storing the results of the computed bids $\hat{b}^{(t)}\in\mathcal{M}_{n\times m}(\mathbb{I})$ in QRAM would require a cost of $\tilde{\mathcal{O}}(mn)$ which would remove all possibility of potential speedups. The same applies to the allocation matrix $\hat{x}^{(t)}$ . Hence, every time we require the usage of $\hat{b}^{(t)}$ oder $\hat{x}^{(t)}$ , we compute them on-the-fly as follows:

\hat{b}_{ij}^{(T)}=\frac{B_{i}^{T+1}v_{ij}^{T}}{m\prod_{t=0}^{T-1}\tilde{p}_{j% }^{(t)}\prod_{k=0}^{T-1}\tilde{\nu}_{i}^{(t)}},\quad\hat{x}_{ij}^{(T)}=\frac{B% _{i}^{T}v_{ij}^{T-1}}{m\prod_{t=0}^{T}\tilde{p}_{j}^{(t)}\prod_{k=0}^{T-1}% \tilde{\nu}_{i}^{(t)}}

(4.1)

Given quantum access to the values of $\Pi_{p}^{(T)}:=\prod_{t=0}^{T}\tilde{p}_{j}^{(t)}$ and $\Pi_{\nu}^{(T)}:=\prod_{t=0}^{T}\tilde{\nu}_{i}^{(t)}$ , one can encode the values of $\hat{b}^{(T)}$ and $\hat{x}^{(T)}$ into a quantum state in superposition via quantum arithmetic circuits in runtime of $\mathcal{O}\left(\operatorname{poly}\log Tmn\right)$ . The quantum access $\Pi_{p}^{(T)}$ and $\Pi_{\nu}^{(T)}$ cost $\mathcal{O}\left(\operatorname{poly}\log mn\right)$ as they are obtained from QRAM, and the operation of taking the $T$ -th power of the budget and value cost $\mathcal{O}\left(\operatorname{poly}\log T\right)$ by binary exponentiation [42].

The remaining steps are to compute the price vector $\tilde{p}$ and utility vector $\tilde{u}$ in each iteration. Each entry $\tilde{p}_{j}$ is the estimation of the $\ell_{1}$ norm of $\hat{b}_{*,j}$ and each entry $\tilde{\nu}_{i}$ is the estimation of the inner product between $\hat{x}_{i,*}$ and $v_{i,*}$ , which can both be obtained using amplitude estimation. $\Pi_{p}^{(t)}$ and $\Pi_{\nu}^{(t)}$ can then be iteratively updated by multiplying by the values of $\tilde{p}_{j}$ and $\tilde{\nu}_{i}$ each iteration. The full algorithm is shown in Algorithm 1.

Input: Quantum access to

B

and

v

, Timestep

T

, Price error

\epsilon_{p}

, Utility error

\epsilon_{\nu}

Output: Query access to the values of a bid matrix estimator

\hat{b}^{(t^{*})}

constructed using values in QRAM and budget and value access

\texttt{maxEGVal}=-\inf

b_{ij}^{(0)}=\frac{B_{i}}{m}

2 for $t=0$ to $T$ do

3 for $j=0$ to $m$ do

\tilde{p}_{j}^{(t)}=(1\pm\epsilon_{p})\left\|\hat{b}_{*,j}^{(t)}\right\|_{1}

via q norm estimation with success prob.

1-\frac{\delta}{2mT}

6 Store vector

\Pi_{p}^{(t)}=\tilde{p}^{(t)}\odot\Pi_{p}^{(t-1)}

into QRAM

7 Gain access to

\hat{x}_{ij}^{(t)}

via

\Pi_{p}^{(t)}

and

\Pi_{\nu}^{(t-1)}

in QRAM

8 for $i=0$ to $n$ do

\tilde{\nu}_{i}^{(t)}=(1\pm\epsilon_{\nu})\left\langle x_{i,*}^{(t)},v_{i,*}\right\rangle

via q inner product estimation with success prob.

1-\frac{\delta}{2nT}

11 Store vector

\Pi_{\nu}^{(t)}=\tilde{\nu}^{(t)}\odot\Pi_{\nu}^{(t-1)}

into QRAM

12 Gain access to

\hat{b}_{ij}^{(t+1)}

via

\Pi_{p}^{(t)}

and

\Pi_{\nu}^{(t)}

in QRAM

13 Classically compute

\tilde{\Phi}^{(t)}=\sum_{i\in[n]}B_{i}\log(\nu_{i}^{(t)})

14 if $\tilde{\Phi}^{(t)}>\texttt{maxEGVal}$ then

\texttt{maxEGVal}=\tilde{\Phi},\,\texttt{bestPiP}=\Pi_{p}^{(t-1)},\,\texttt{% bestPiNu}=\Pi_{\nu}^{(t-1)}

18return bestPiP and bestPiNu in QRAM

Algorithm 1 Quantum faulty proportional response dynamics

Theorem 4.1 (Quantum algorithm for faulty proportional response dynamics).

Let $\delta\in(0,0.5),n,m,T\in\mathbb{N}$ , $\epsilon_{p}\leq\frac{\log m}{8T}$ , and $\epsilon_{\nu}\leq\frac{1}{6T}$ . Given quantum access to $B$ and $v$ , and access to QRAM, with success probability $1-\delta$ , Algorithm 1 produces values stored in QRAM such that query access to the values of $\hat{b}^{(t^{*})}$ can be constructed, where

\Phi(\hat{b}^{(t^{*})})-\Phi(b^{*})\leq\frac{2\log m}{T},

with $\tilde{\mathcal{O}}(T^{2}\sqrt{mn\max(m,n)}\log\nicefrac{{1}}{{\delta}})$ runtime and $\tilde{\mathcal{O}}(m+n)$ QRAM space. To provide query access to $\hat{b}^{(t^{*})}$ , an cost of $\mathcal{O}\left(\operatorname{poly}\log Tmn\right)$ is incurred from accessing $\Pi_{p}^{(t^{*}-1)}$ and $\Pi_{\nu}^{(t^{*}-1)}$ in QRAM.

Proof.

Per union bound [57], we find that the total success probability is at least $1-\delta$ . Note that the output of Algorithm 1 of bestPiP and bestPiNu corresponds to the values of $\Pi_{p}^{(t^{*}-1)}$ and $\Pi_{\nu}^{(t^{*}-1)}$ that can be used to construct $b^{(t^{*})}$ per Equation 4.1. This gives us the guarantee of convergence shown in Theorem 3.2.

Moving to the runtime analysis, the quantum norm estimation subroutine takes $\mathcal{O}(T^{2}\sqrt{n}\log\frac{mT}{\delta})$ for $mT$ iterations, while quantum inner product estimation takes $\mathcal{O}(T^{2}\sqrt{m}\log\frac{nT}{\delta})$ for $nT$ iterations, resulting in a runtime of $\tilde{\mathcal{O}}(T^{2}\sqrt{mn\max(m,n)}\log\frac{1}{\delta})$ . For uses of QRAM, the construction on Lines 5 and 9, is a one-time cost of $\tilde{\mathcal{O}}(n)$ and $\tilde{\mathcal{O}}(m)$ , respectively, with a total runtime of $\tilde{\mathcal{O}}(T(m+n))$ . The classical computation of the EG value in Line 11 costs $\mathcal{O}(Tn)$ . We note that the quantum norm and inner product estimation subroutine is the main bottleneck of the algorithm, and hence the total runtime is then $\tilde{\mathcal{O}}(T^{2}\sqrt{mn\max(m,n)}\log\frac{1}{\delta})$ .

For the memory complexity, for the $t$ -th iteration, we require 6 vectors in QRAM: the current iteration $\Pi_{p}^{(t)}$ and $\Pi_{\nu}^{(t)}$ , the best iteration bestPiP and bestPiNu and the previous iteration $\Pi_{p}^{(t-1)}$ and $\Pi_{\nu}^{(t-1)}$ , in case we need to update bestPiP and bestPiNu. Note that to update the best iteration, we simply reroute the register of the previous iteration to being the best iteration. There is no need to copy data or reconstruct a new QRAM as the data from the previous iteration is no longer needed in the next iteration. Therefore, the memory is $\tilde{\mathcal{O}}(m+n)$ for storing the 6 vectors. ∎

5 Numerical simulations

We simulate the market equilibrium computation under PR dynamics and our quantum algorithm. To showcase the effects of quantum speedups, we fixed the number of queries to all bid matrices $b^{(t)}$ and observed the reduction of the objective value over the number of queries.

As an actual simulation of amplitude estimation using quantum gates over multiple qubits is costly, we directly compute the probability vector of $\Pr[Z=z]$ for $z\in[M]$ that one would obtain by amplitude estimation [32] for a target value $a$ ,

\Pr[Z=z]=\frac{\sin^{2}(M\Delta_{z}\pi)}{M^{2}\sin^{2}(\Delta_{z}\pi)}

(5.1)

where $\Delta_{z}=\min(\lvert z-\sin^{-1}(\sqrt{a})/\pi\rvert,\lvert 1-z+\sin^{-1}(% \sqrt{a})/\pi\rvert)$ , and $M$ is the number of times that call the unitaries $U$ and $V$ in QAE (see Theorem 2.1), and is linearly correlated to the runtime We then sample the output according to the computed probabilities to obtain an estimator $\tilde{a}=\sin^{2}(\pi\frac{z}{M})$ .

For our experiments, we generate data the input data $v$ where the value $v$ is sampled from a uniform distribution with range $[0,1)$ and a normal distribution $\mathcal{N}(0.5,0.25)$ , where we resample values that fall outside the range of $[0,1]$ . For the budget $B$ , we either sample from the same distribution as the value matrix or set the same budget for all buyers to simulate competitive equilibrium from equal income (CEEI) applications. Our simulation includes $n=16384$ buyers, $m=16384$ goods, and iterate for $T=16$ iterations for the PR dynamics. For the quantum algorithm, note that the queries per iteration would be reduced by $\sqrt{n}$ if we use an actual quantum computer, hence increasing the number of iterations to fix the number of queries. For amplitude estimation, we run for $\sqrt{T\sqrt{n}}=512$ iterations and set $M\in\mathcal{O}(\sqrt{T\sqrt{n}})$ . As the classical algorithms are deterministic, we rerun our quantum algorithm over $15$ times with the same sample of $B$ and $v$ to observe the variance of convergence progress. Experimental results are shown in Figure 1(a). Details on implementation and further experimental setup are found in Appendix D.

From the plots of Figure 1, we note that the results fit our theoretical results in that the quantum algorithm converges much faster than that of the PR dynamics [24]. Further, we also compare against the convergence of projected gradient descent, which supports empirical results by Gao and Kroer [26] that in the regime of mid-level accuracy and low iterations, PR dynamics-related algorithms, both classical and quantum, converge faster than projected gradient descent.

6 Discussion

Quasi-linear utilities.

For the bulk of our paper, we focus on the setting of linear utilities for Fisher markets. However, applications of market equilibrium computation in large-scale Fisher markets involve mostly quasi-linear utilities [9]. An approach for using PR dynamics for quasi-linear utilities proposed by Gao and Kroer [26]⁴⁴4There is another method proposed by Cheung et al. [58], which we find difficult to convert to quantum due to its use of thresholding, which would cause problems with faulty updates from the FPR dynamics. includes the usage of slack variables $\delta=(\delta_{1},\cdots,\delta_{m})$ that represent the buyers’ leftover budgets. The PR updates are then modified as follows:

b_{ij}^{(t+1)}=B_{i}\frac{v_{ij}x_{ij}^{(t)}}{\sum_{j^{\prime}}v_{ij^{\prime}}% x_{ij^{\prime}}^{(t)}+\delta_{i}^{(t)}},\;\delta_{i}^{(t+1)}=B_{i}\frac{\delta% _{i}^{(t)}}{\sum_{j^{\prime}}v_{ij^{\prime}}x_{ij^{\prime}}^{(t)}+\delta_{i}^{% (t)}}.

(6.1)

Further, PR dynamics for quasi-linear utilities exhibit a convergence rate of $\mathcal{O}(\log(m+1)/T)$ . Using the methods discussed in previous sections, the quasi-linear version of PR dynamics can then be readily adapted to its quantum version by employing the same techniques of computing and storing in QRAM the values of $\Pi_{p}^{(t)}$ and $\Pi_{\nu}^{(t)}$ in conjunction with on-the-fly computation of $\hat{b}^{(t)}$ , $\hat{x}^{(t)}$ and $\hat{\delta}^{(t)}$ .

Constant number of buyers.

Notice that our quantum algorithm provides a quadratic speedup on the smaller value in regards to the number of buyers $n$ and number of goods $m$ . Therefore given extreme cases where the number of buyers $n\in\mathcal{O}(1)$ , our algorithm does not provide a speedup. However, in such cases, quantum speedups may still be obtained simply by removing the amplitude estimation step for estimating the price for each item and replacing it with using quantum arithmetic circuits to compute the exact sum. We use a total of $\mathcal{O}(nT\operatorname{poly}\log(T,m,n))$ qubits to compute the values of $b_{i,*}$ separately on the $i$ -th set of qubits, and only conduct amplitude estimation when estimating the utility value for each buyer. Given that in this setting, $n\in\mathcal{O}(1)$ , the total runtime would then be $\mathcal{O}\left(T^{2}\sqrt{m}\log\nicefrac{{1}}{{\delta}}\right)$ , gaining a quadratic speedup over the number of goods $m$ .

Dequantization.

Given the work in recent years towards the development of quantum-inspired classical algorithms [59, 60, 61, 62] that achieve similar performances as quantum algorithms using sampling-based techniques, a natural question that arises is whether our algorithm can be “de-quantized”. The main speedup in our algorithm stems from the usage of estimation of $\ell_{1}$ norms and inner products. While the use of sampling techniques can indeed provide inner product estimations, they retain the same $\mathcal{O}(1/\epsilon^{2})$ dependency instead of the $\mathcal{O}(1/\epsilon)$ dependency of QAE. Hence, our algorithm performance may be hard to replicate in classical settings.

On the other hand, while it has been suggested that the computation of market equilibrium may benefit from low-rank approximations [15], methods of using such properties to accelerate the computation of gradients have not been proposed, given that the update of the PR dynamics rely on element-wise multiplication of matrices instead matrix multiplication. This would suggest that using sampling techniques to accelerate updates would be similarly difficult.

Potential and limitations for further quantum speedups.

Our quantum algorithm shares similarities to other quantum algorithms that are based on the multiplicative weight update (MWU) method [63, 64]. Such methods have found success in obtaining quantum speedups for LPs [65] and SDPs [66, 67, 68, 69], which have been extended to applications such as zero-sum games [65, 70], quadratic binary optimization [71], and financial applications [28, 72]. Apart from the MWU-esque PR dynamics, various other methods for computing market equilibrium have also been proposed. Can quantum speedups obtained from these methods exceed those of our quantum algorithm?

Tracing back to the roots of the EG convex program [22, 23] and Shmyrev convex program [51], it is well known that such programs can be solved in polynomial time with interior-point methods (IPM) [50]. However, as IPMs require using linear solvers as subroutines, and as there is no guarantee of well-conditioned systems, the quantum linear systems solver [73, 74] may not provide significant speedup. Therefore, it may be unlikely that quantum IPMs [75] can provide significant speedups.

First-order methods such as the Frank-Wolfe (FW) algorithm [76] and projected gradient descent (PGD) have also been discussed as candidates for solving market equilibrium [26], with PGD achieving linear convergence classically. While PGD obtains a superior asymptotic convergence rate in terms of the error $\epsilon$ compared to PR dynamics, as our quantum speedups stem from faster computations of results within a single iteration, it may be harder to find such speedups for PGD as there has been no evidence for quantum speedups in projections onto a simplex [77, 78] as required.

On the other hand, the FW algorithm has been shown to provide quantum speedups for regression [79, 80]. However, convergence results of FW [81, 82] show that $\Phi(b^{(T)})-\Phi(b^{*})\leq C_{\Phi}/(T+2)$ , where $C_{\Phi}$ can be shown to be $\mathcal{O}(n)$ by computing relevant values. The number of iterations $T$ required for convergence to additive error $\varepsilon$ is then $\mathcal{O}(n/\varepsilon)$ as compared to $\mathcal{O}(\log m/\varepsilon)$ of PR dynamics. This matches the results of Gao and Kroer [26], which show that FW has slow convergence empirically for market equilibrium computation. Prior no-go results suggest that quantum algorithms cannot provide speedups for the number of iterations $T$ when $T$ is independent of the problem dimension [83, 84]. Assuming no quantum speedups in the number of iterations, given the $\mathcal{O}(n)$ upper bound in the FW algorithm, the quantum algorithm based on FW can potentially have a higher dependency on $n$ than the classical PR dynamics.

Lastly, we ask whether random sampling of buyers or goods can provide further speedups. Classical results in first-order updates for randomly sampled buyers and goods [85] indicate that the number of iterations would increase multiplicatively by $m$ and $n$ , respectively, such that the total runtime cost of the algorithms remain at $\tilde{\mathcal{O}}(mn)$ . Assuming no quantum speedups on the number of iterations, further quantum speedups by incorporating sampling may be difficult.

Acknowledgments and Disclosure of Funding

The authors thank Gregory Kang Ruey Lau for discussions. This work is supported by the National Research Foundation, Singapore, and A*STAR under its CQT Bridging Grant and its Quantum Engineering Programme under grant NRF2021-QEP2-02-P05.

References

Arrow [1951] K. J. Arrow, An extension of the basic theorems of classical welfare economics, in Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, Berkeley Symp. on Math. Statist. and Prob (University of California Press, Berkeley, California, USA, 1951) pp. 507–532.
Debreu [1951] G. Debreu, The coefficient of resource utilization, Econometrica 19, 273 (1951).
Arrow and Debreu [1954] K. J. Arrow and G. Debreu, Existence of an equilibrium for a competitive economy, Econometrica 22, 265 (1954).
Vazirani [2007] V. V. Vazirani, Combinatorial algorithms for market equilibria, in Algorithmic Game Theory, edited by N. Nisan, T. Roughgarden, E. Tardos, and V. V. Vazirani (Cambridge University Press, 2007) p. 103–134.
Codenotti and Varadarajan [2007] B. Codenotti and K. Varadarajan, Computation of market equilibria by convex programming, in Algorithmic Game Theory, edited by N. Nisan, T. Roughgarden, E. Tardos, and V. V. Vazirani (Cambridge University Press, 2007) p. 135–158.
Mas-Colell et al. [1995] A. Mas-Colell, M. D. Whinston, and J. R. Green, Equilibrium and its basic welfare properties, in Microeconomic theory (Oxford Univ. Press, 1995) Chap. 16.
Foley [1967] D. K. Foley, Resource allocation and the public sector, Ph.D. thesis, Yale University (1967).
Varian [1974] H. R. Varian, Equity, envy, and efficiency, Journal of Economic Theory 9, 63–91 (1974).
Kroer and Stier-Moses [2021] C. Kroer and N. E. Stier-Moses, Market equilibrium models in large-scale internet markets, in Springer Series in Supply Chain Management (Springer International Publishing, 2021) p. 147–189.
Conitzer et al. [2022] V. Conitzer, C. Kroer, D. Panigrahi, O. Schrijvers, N. E. Stier-Moses, E. Sodomka, and C. A. Wilkens, Pacing equilibrium in first price auction markets, Management Science 68, 8515–8535 (2022).
Othman et al. [2010] A. Othman, T. Sandholm, and E. Budish, Finding approximate competitive equilibria: efficient and fair course allocation, in Proceedings of the 9th International Conference on Autonomous Agents and Multiagent Systems: Volume 1 - Volume 1, AAMAS ’10 (International Foundation for Autonomous Agents and Multiagent Systems, Richland, SC, 2010) p. 873–880.
Budish et al. [2017] E. Budish, G. P. Cachon, J. B. Kessler, and A. Othman, Course match: A large-scale implementation of approximate competitive equilibrium from equal incomes for combinatorial allocation, Operations Research 65, 314–336 (2017).
Babaioff et al. [2019] M. Babaioff, N. Nisan, and I. Talgam-Cohen, Fair allocation through competitive equilibrium from generic incomes, in Proceedings of the Conference on Fairness, Accountability, and Transparency, FAT* ’19 (Association for Computing Machinery, New York, NY, USA, 2019) p. 180.
Im et al. [2017] S. Im, J. Kulkarni, and K. Munagala, Competitive algorithms from competitive equilibria: Non-clairvoyant scheduling under polyhedral constraints, J. ACM 65 (2017).
Kroer et al. [2022] C. Kroer, A. Peysakhovich, E. Sodomka, and N. E. Stier-Moses, Computing large market equilibria using abstractions, Operations Research 70, 329–351 (2022).
Fisher [1891] I. Fisher, Mathematical investigations in the theory of value and prices, Ph.D. thesis, Yale University (1891).
Brainard and Scarf [2005] W. C. Brainard and H. E. Scarf, How to compute equilibrium prices in 1891, The American Journal of Economics and Sociology 64, 57–83 (2005).
Scarf [1967] H. E. Scarf, The core of an $n$ person game, Econometrica 35, 50 (1967).
Devanur et al. [2008] N. R. Devanur, C. H. Papadimitriou, A. Saberi, and V. V. Vazirani, Market equilibrium via a primal-dual algorithm for a convex program, J. ACM 55 (2008).
Orlin [2010] J. B. Orlin, Improved algorithms for computing Fisher’s market clearing prices, in Proceedings of the Forty-Second ACM Symposium on Theory of Computing, STOC ’10 (Association for Computing Machinery, New York, NY, USA, 2010) p. 291–300.
Végh [2012] L. A. Végh, Strongly polynomial algorithm for a class of minimum-cost flow problems with separable convex objectives, in Proceedings of the Forty-Fourth Annual ACM Symposium on Theory of Computing, STOC ’12 (Association for Computing Machinery, New York, NY, USA, 2012) p. 27–40.
Eisenberg and Gale [1959] E. Eisenberg and D. Gale, Consensus of subjective probabilities: The pari-mutuel method, The Annals of Mathematical Statistics 30, 165–168 (1959).
Eisenberg [1961] E. Eisenberg, Aggregation of utility functions, Management Science 7, 337–350 (1961).
Wu and Zhang [2007] F. Wu and L. Zhang, Proportional response dynamics leads to market equilibrium, in Proceedings of the Thirty-Ninth Annual ACM Symposium on Theory of Computing, STOC ’07 (Association for Computing Machinery, New York, NY, USA, 2007) p. 354–363.
Zhang [2011] L. Zhang, Proportional response dynamics in the Fisher market, Theoretical Computer Science 412, 2691–2698 (2011).
Gao and Kroer [2020] Y. Gao and C. Kroer, First-order methods for large-scale market equilibrium computation, in Advances in Neural Information Processing Systems, Vol. 33, edited by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Curran Associates, Inc., 2020) pp. 21738–21750.
Li et al. [2019] T. Li, S. Chakrabarti, and X. Wu, Sublinear quantum algorithms for training linear and kernel-based classifiers, in Proceedings of the 36th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 97, edited by K. Chaudhuri and R. Salakhutdinov (PMLR, 2019) pp. 3815–3824.
Rebentrost et al. [2021] P. Rebentrost, Y. Hamoudi, M. Ray, X. Wang, S. Yang, and M. Santha, Quantum algorithms for hedging and the learning of Ising models, Phys. Rev. A 103, 012418 (2021).
Montanaro [2016] A. Montanaro, Quantum algorithms: an overview, npj Quantum Information 2 (2016).
Dalzell et al. [2023] A. M. Dalzell, S. McArdle, M. Berta, P. Bienias, C.-F. Chen, A. Gilyén, C. T. Hann, M. J. Kastoryano, E. T. Khabiboulline, A. Kubica, G. Salton, S. Wang, and F. G. S. L. Brandão, Quantum algorithms: A survey of applications and end-to-end complexities (2023), arXiv:2310.03011 [quant-ph] .
Abbas et al. [2023] A. Abbas, A. Ambainis, B. Augustino, A. Baertschi, H. Buhrman, C. Coffrin, G. Cortiana, V. Dunjko, D. Egger, B. Elmegreen, N. Franco, F. Fratini, B. Fuller, J. Gacon, C. Gonciulea, S. Gribling, S. Gupta, S. Hadfield, R. Heese, G. Kircher, T. Kleinert, T. Koch, G. Korpas, S. Lenk, J. Marecek, V. Markov, G. Mazzola, S. Mensa, N. Mohseni, G. Nannicini, C. O’Meara, E. Peña Tapia, S. Pokutta, M. Proissl, P. Rebentrost, E. Sahin, B. Symons, S. Tornow, V. Valls, S. Woerner, M. Wolf-Bauwens, J. Yard, S. Yarkoni, D. Zechiel, S. Zhuk, and C. Zoufal, Quantum Optimization: Potential, Challenges, and the Path Forward, Tech. Rep. (Office of Scientific and Technical Information (OSTI), 2023).
Brassard et al. [2002] G. Brassard, P. Høyer, M. Mosca, and A. Tapp, Quantum amplitude amplification and estimation, in Quantum computation and information, Contemporary Mathematics, Vol. 305 (American Mathematical Society, Providence, RI, USA, 2002) pp. 53–74.
Suzuki et al. [2020] Y. Suzuki, S. Uno, R. Raymond, T. Tanaka, T. Onodera, and N. Yamamoto, Amplitude estimation without phase estimation, Quantum Information Processing 19 (2020).
Grinko et al. [2021] D. Grinko, J. Gacon, C. Zoufal, and S. Woerner, Iterative quantum amplitude estimation, npj Quantum Information 7 (2021).
Nakaji [2020] K. Nakaji, Faster amplitude estimation, Quantum Information and Computation 20, 1109–1123 (2020).
Harrow and Wei [2020] A. W. Harrow and A. Y. Wei, Adaptive quantum simulated annealing for Bayesian inference and estimating partition functions, in Proceedings of the Thirty-First Annual ACM-SIAM Symposium on Discrete Algorithms, SODA ’20 (Society for Industrial and Applied Mathematics, USA, 2020) p. 193–212.
Rall and Fuller [2023] P. Rall and B. Fuller, Amplitude estimation from quantum signal processing, Quantum 7, 937 (2023).
Cornelissen and Hamoudi [2023] A. Cornelissen and Y. Hamoudi, A sublinear-time quantum algorithm for approximating partition functions, in Proceedings of the 2023 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA) (Society for Industrial and Applied Mathematics, 2023) pp. 1245–1264.
Montanaro [2015] A. Montanaro, Quantum speedup of Monte Carlo methods, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471, 20150301 (2015).
Vedral et al. [1996] V. Vedral, A. Barenco, and A. Ekert, Quantum networks for elementary arithmetic operations, Phys. Rev. A 54, 147–153 (1996).
Takahashi [2009] Y. Takahashi, Quantum arithmetic circuits: A survey, IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences E92-A, 1276–1283 (2009).
Montgomery [1987] P. L. Montgomery, Speeding the Pollard and elliptic curve methods of factorization, Mathematics of Computation 48, 243–264 (1987).
Kuperberg [2013] G. Kuperberg, Another subexponential-time quantum algorithm for the dihedral hidden subgroup problem, in 8th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2013), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 22 (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2013) pp. 20–34.
Jaques and Rattew [2023] S. Jaques and A. G. Rattew, QRAM: A survey and critique (2023), arXiv:2305.10310 [quant-ph] .
Babbush et al. [2018] R. Babbush, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, A. Paler, A. Fowler, and H. Neven, Encoding electronic spectra in quantum circuits with linear t complexity, Phys. Rev. X 8, 041015 (2018).
Ambainis [2007] A. Ambainis, Quantum walk algorithm for element distinctness, SIAM Journal on Computing 37, 210–239 (2007).
Giovannetti et al. [2008a] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum random access memory, Phys. Rev. Lett. 100, 160501 (2008a).
Giovannetti et al. [2008b] V. Giovannetti, S. Lloyd, and L. Maccone, Architectures for a quantum random access memory, Phys. Rev. A 78, 052310 (2008b).
Jain and Vazirani [2010] K. Jain and V. V. Vazirani, Eisenberg–gale markets: Algorithms and game-theoretic properties, Games and Economic Behavior 70, 84–106 (2010).
Boyd and Vandenberghe [2004] S. Boyd and L. Vandenberghe, Interior-point methods, in Convex Optimization (Cambridge University Press, 2004) p. 561–630.
Shmyrev [2009] V. I. Shmyrev, An algorithm for finding equilibrium in the linear exchange model with fixed budgets, Journal of Applied and Industrial Mathematics 3, 505–518 (2009).
Levin et al. [2008] D. Levin, K. LaCurts, N. Spring, and B. Bhattacharjee, Bittorrent is an auction: Analyzing and improving bittorrent’s incentives, SIGCOMM Comput. Commun. Rev. 38, 243–254 (2008).
Birnbaum et al. [2011] B. Birnbaum, N. R. Devanur, and L. Xiao, Distributed algorithms via gradient descent for Fisher markets, in Proceedings of the 12th ACM Conference on Electronic Commerce, EC ’11 (Association for Computing Machinery, New York, NY, USA, 2011) p. 127–136.
Nemirovsky and Yudin [1983] A. S. Nemirovsky and D. B. Yudin, Problem complexity and method efficiency in optimization (Wiley, 1983).
Beck and Teboulle [2003] A. Beck and M. Teboulle, Mirror descent and nonlinear projected subgradient methods for convex optimization, Operations Research Letters 31, 167–175 (2003).
Bregman [1967] L. Bregman, The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming, USSR Computational Mathematics and Mathematical Physics 7, 200–217 (1967).
Boole [1847] G. Boole, The Mathematical Analysis of Logic (Cambridge University Press, 1847).
Cheung et al. [2021] Y. K. Cheung, S. Leonardos, and G. Piliouras, Learning in markets: Greed leads to chaos but following the price is right, in Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21, edited by Z.-H. Zhou (International Joint Conferences on Artificial Intelligence Organization, 2021) pp. 111–117.
Tang [2019] E. Tang, A quantum-inspired classical algorithm for recommendation systems, in Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, STOC 2019 (Association for Computing Machinery, New York, NY, USA, 2019) p. 217–228.
Arrazola et al. [2020] J. M. Arrazola, A. Delgado, B. R. Bardhan, and S. Lloyd, Quantum-inspired algorithms in practice, Quantum 4, 307 (2020).
Tang [2021] E. Tang, Quantum principal component analysis only achieves an exponential speedup because of its state preparation assumptions, Phys. Rev. Lett. 127, 060503 (2021).
Chia et al. [2022] N.-H. Chia, A. P. Gilyén, T. Li, H.-H. Lin, E. Tang, and C. Wang, Sampling-based sublinear low-rank matrix arithmetic framework for dequantizing quantum machine learning, J. ACM 69 (2022).
Arora et al. [2005] S. Arora, E. Hazan, and S. Kale, Fast algorithms for approximate semidefinite programming using the multiplicative weights update method, in 46th Annual IEEE Symposium on Foundations of Computer Science (FOCS’05) (2005) pp. 339–348.
Arora et al. [2012] S. Arora, E. Hazan, and S. Kale, The multiplicative weights update method: a meta-algorithm and applications, Theory of Computing 8, 121–164 (2012).
van Apeldoorn and Gilyén [2019] J. van Apeldoorn and A. Gilyén, Quantum algorithms for zero-sum games (2019), arXiv:1904.03180 [quant-ph] .
Brandão and Svore [2017] F. G. S. L. Brandão and K. M. Svore, Quantum speed-ups for solving semidefinite programs, in 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS) (IEEE Computer Society, Los Alamitos, CA, USA, 2017) pp. 415–426.
Brandão et al. [2019] F. G. S. L. Brandão, A. Kalev, T. Li, C. Y.-Y. Lin, K. M. Svore, and X. Wu, Quantum sdp solvers: Large speed-ups, optimality, and applications to quantum learning, in 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 132, edited by C. Baier, I. Chatzigiannakis, P. Flocchini, and S. Leonardi (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2019) pp. 27:1–27:14.
van Apeldoorn and Gilyén [2019] J. van Apeldoorn and A. Gilyén, Improvements in quantum sdp-solving with applications, in 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 132, edited by C. Baier, I. Chatzigiannakis, P. Flocchini, and S. Leonardi (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2019) pp. 99:1–99:15.
van Apeldoorn et al. [2020] J. van Apeldoorn, A. Gilyén, S. Gribling, and R. de Wolf, Quantum SDP-Solvers: Better upper and lower bounds, Quantum 4, 230 (2020).
Jain et al. [2022] R. Jain, G. Piliouras, and R. Sim, Matrix multiplicative weights updates in quantum zero-sum games: Conservation laws & recurrence, in Advances in Neural Information Processing Systems, Vol. 35, edited by S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (Curran Associates, Inc., 2022) pp. 4123–4135.
G.S L. Brandão et al. [2022] F. G.S L. Brandão, R. Kueng, and D. Stilck França, Faster quantum and classical SDP approximations for quadratic binary optimization, Quantum 6, 625 (2022).
Lim and Rebentrost [2024] D. Lim and P. Rebentrost, A quantum online portfolio optimization algorithm, Quantum Information Processing 23 (2024).
Harrow et al. [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Phys. Rev. Lett. 103, 150502 (2009).
Childs et al. [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46, 1920–1950 (2017).
Kerenidis and Prakash [2020] I. Kerenidis and A. Prakash, A quantum interior point method for LPs and SDPs, ACM Transactions on Quantum Computing 1 (2020).
Frank and Wolfe [1956] M. Frank and P. Wolfe, An algorithm for quadratic programming, Naval Research Logistics Quarterly 3, 95–110 (1956).
Duchi et al. [2008] J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra, Efficient projections onto the $\ell_{1}$ -ball for learning in high dimensions, in Proceedings of the 25th International Conference on Machine Learning, ICML ’08 (Association for Computing Machinery, New York, NY, USA, 2008) p. 272–279.
Condat [2015] L. Condat, Fast projection onto the simplex and the $\ell_{1}$ ball, Mathematical Programming 158, 575–585 (2015).
Du et al. [2022] Y. Du, M.-H. Hsieh, T. Liu, S. You, and D. Tao, Quantum differentially private sparse regression learning, IEEE Transactions on Information Theory 68, 5217–5233 (2022).
Chen and de Wolf [2023] Y. Chen and R. de Wolf, Quantum algorithms and lower bounds for linear regression with norm constraints, in 50th International Colloquium on Automata, Languages, and Programming (ICALP 2023), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 261, edited by K. Etessami, U. Feige, and G. Puppis (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2023) pp. 38:1–38:21.
Clarkson [2010] K. L. Clarkson, Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm, ACM Trans. Algorithms 6 (2010).
Jaggi [2013] M. Jaggi, Revisiting Frank-Wolfe: Projection-free sparse convex optimization, in Proceedings of the 30th International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 28, edited by S. Dasgupta and D. McAllester (PMLR, Atlanta, Georgia, USA, 2013) pp. 427–435.
Garg et al. [2021a] A. Garg, R. Kothari, P. Netrapalli, and S. Sherif, No quantum speedup over gradient descent for non-smooth convex optimization, in 12th Innovations in Theoretical Computer Science Conference (ITCS 2021), Leibniz International Proceedings in Informatics (LIPIcs), Vol. 185, edited by J. R. Lee (Schloss Dagstuhl – Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2021) pp. 53:1–53:20.
Garg et al. [2021b] A. Garg, R. Kothari, P. Netrapalli, and S. Sherif, Near-optimal lower bounds for convex optimization for all orders of smoothness, in Advances in Neural Information Processing Systems, Vol. 34, edited by M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan (Curran Associates, Inc., 2021) pp. 29874–29884.
Nan et al. [2023] T. Nan, Y. Gao, and C. Kroer, Fast and interpretable dynamics for Fisher markets via block-coordinate updates, Proceedings of the AAAI Conference on Artificial Intelligence 37, 5832–5840 (2023).
Jerrum et al. [1986] M. R. Jerrum, L. G. Valiant, and V. V. Vazirani, Random generation of combinatorial structures from a uniform distribution, Theoretical Computer Science 43, 169–188 (1986).
Durr and Høyer [1996] C. Durr and P. Høyer, A quantum algorithm for finding the minimum (1996), arXiv:quant-ph/9607014 [quant-ph] .
Paszke et al. [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, PyTorch: An imperative style, high-performance deep learning library, in Advances in Neural Information Processing Systems, Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).

Appendix A Quantum subroutines

In this section, we show prior results that obtain $\ell_{1}$ norms and inner products with quadratic speedups with amplitude estimation.

Lemma A.1 (Quantum state preparation and norm estimation; Lemma 5, [28]).

Let $n\in\mathbb{N}$ . We are given quantum query access to non-zero vector $w\in\mathbb{I}^{n}$ , with $\max_{j}w_{j}=1$ .

1.

There exists a quantum circuit that prepares the state $\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\boldsymbol{\ket{j}}\left(\sqrt{w_{j}}\ket{0}+% \sqrt{1-w_{j}}\ket{1}\right)$ with two queries and $\mathcal{O}\left(\log n\right)$ gates.
2.

Let $\epsilon>0$ and $\delta\in(0,1)$ . There exists a quantum algorithm that provides an estimate $\Gamma_{w}$ of the $\ell_{1}$ -norm $\|w\|_{1}$ such that $\left|\|w\|_{1}-\Gamma_{w}\right|\leq\epsilon\|w\|_{1}$ , with probability at least $1-\delta$ . The algorithm requires $\mathcal{O}\left(\frac{\sqrt{n}}{\epsilon}\log\frac{1}{\delta}\right)$ queries and $\tilde{\mathcal{O}}\left(\frac{\sqrt{n}}{\epsilon}\log\frac{1}{\delta}\right)$ quantum gates.

Proof.

We reiterate the proof of Lemma 5 in [28] for the convenience of the reader.

1.

First, using $\mathcal{O}(\log n)$ Hadamard gates, prepare the state $\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\boldsymbol{\ket{j}}$ . Then, by quantum query access to $w$ , obtain $\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\boldsymbol{\ket{j}}\ket{w_{j}}$ . By controlled rotation gates, we can then obtain $\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\boldsymbol{\ket{j}}\ket{w_{j}}(\sqrt{w_{j}}% \ket{0}+\sqrt{1-w_{j}}\ket{1})$ . By another quantum query access to $w$ , we can uncompute the intermediate registers and obtain $\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\boldsymbol{\ket{j}}(\sqrt{w_{j}}\ket{0}+\sqrt% {1-w_{j}}\ket{1})$ .

First observe that that with projector $P=I_{n}\otimes\lvert 0\rangle\langle 0\rvert$ and $\ket{\psi}=\frac{1}{\sqrt{n}}\sum_{j=1}^{n}\boldsymbol{\ket{j}}(\sqrt{w_{j}}% \ket{0}+\sqrt{1-w_{j}}\ket{1})$ , one can obtain $a=\|P\ket{\psi}\|_{2}^{2}=\frac{\|w\|_{1}}{n}$ . Setting $M\geq\frac{6\pi}{\epsilon}\sqrt{N}$ , we obtain an estimate

|\tilde{a}_{\mathrm{est}}-a|\leq 2\pi\frac{\sqrt{a(1-a)}}{M}+\frac{\pi^{2}}{M^% {2}}\leq\frac{\epsilon}{6\sqrt{N}}\left(2\sqrt{a}+\frac{\epsilon}{12}\right)% \leq\frac{3\sqrt{a}\epsilon}{6\sqrt{N}}\leq\frac{\sqrt{\|w\|_{1}}\cdot\epsilon% }{2N}\leq\frac{a}{2}\cdot\epsilon

(A.1)

with probability at least $\frac{8}{\pi^{2}}$ . Using the powering lemma [86], we can boost the success probability to $1-\delta$ by taking the median of $O(\log 1/\delta)$ runs of the QAE algorithm.

∎

Remark A.1.1.

Note that Lemma A.1 has the requirement that $\max_{j}w_{j}=1$ . For cases where this is not the case, we can use a maximum finding algorithm to divide all entries by the largest value. Such can be achieved by the following quantum minimum/maximum finding algorithm in $\mathcal{O}(\sqrt{n})$ runtime, which we introduce below. Recall that division takes $\mathcal{O}(1)$ runtime with quantum arithmetic circuits.

Lemma A.2 (Quantum minimum finding; Theorem 1, [87]).

Let $n\in\mathbb{N}$ . Given quantum query access to non-zero vector $w\in\mathbb{I}^{n}$ , we can find the minimum $w_{\min}=\min_{j\in[n]}w_{j}$ with success probability $1-\delta$ with $\mathcal{O}\left(\sqrt{n}\log\frac{1}{\delta}\right)$ queries and $\tilde{\mathcal{O}}\left(\sqrt{n}\log\frac{1}{\delta}\right)$ quantum gates.

Corollary A.2.1 (Quantum maximum finding).

Let $n\in\mathbb{N}$ . Given quantum query access to non-zero vector $w\in\mathbb{I}^{n}$ , we can find the maximum $w_{\max}=\min_{j\in[n]}w_{j}$ with success probability $1-\delta$ with $\mathcal{O}\left(\sqrt{n}\log\frac{1}{\delta}\right)$ queries and $\tilde{\mathcal{O}}\left(\sqrt{n}\log\frac{1}{\delta}\right)$ quantum gates.

Below we present a quantum inner product estimation algorithm simplified from Lemma 6 of [28].

Lemma A.3 (Quantum inner product estimation with relative accuracy).

Let $n\in\mathbb{N}$ , $\epsilon<0$ and $\delta\in(0,1)$ . We are given quantum query access to two vectors $u,v\in\mathbb{I}^{n}$ . An estimate $\Gamma_{u,v}$ for the inner product can be provided such that $|\Gamma_{u,v}-u\cdot v|\leq\epsilon\ u\cdot v$ with success probability $1-\delta$ . This estimate is obtained with $\mathcal{O}\left(\frac{\sqrt{n}}{\epsilon}\log\frac{1}{\delta}\right)$ queries and $\tilde{\mathcal{O}}\left(\frac{\sqrt{n}}{\epsilon}\log\frac{1}{\delta}\right)$ quantum gates.

Proof.

Using quantum arithmetic circuits, we can obtain $z_{j}=u_{j}v_{j}$ , i.e., $z=u\odot v$ , by the following:

\boldsymbol{\ket{j}}\to\boldsymbol{\ket{j}}\ket{u_{j}}\ket{v_{j}}\to% \boldsymbol{\ket{j}}\ket{z_{j}}\ket{v_{j}}\to\boldsymbol{\ket{j}}\ket{z_{j}}% \ket{\bar{0}}

(A.2)

Using quantum maximum finding in Corollary A.2.1 to find $z_{\max}$ up to probability $1-\delta/2$ , we can then obtain $z_{j}/z_{\max}$ . Lastly, using Lemma A.1, we can obtain $\Gamma_{u,v}=\Gamma_{z}$ such that $|\Gamma_{u,v}-u\cdot v|=|\Gamma_{z}-\|z\|_{1}|\leq\epsilon u\cdot v$ up to probability $1-\delta/2$ . Using a union bound [57], we find the total success probability of the entire process is $1-\delta$ . ∎

Appendix B Convergence guarantees for the PR dynamics

We show the convergence guarantee of the proportional response (PR) dynamics in regards to the Eisenberg-Gale convex program by Zhang [25] and improved upon by Birnbaum et al. [53], and the convergence in regards to the Shmyrev convex program, first shown also by Birnbaum et al. [53]. Recall that the negative target function from the Eisenberg-Gale convex program is

\Phi(b)=-\sum_{i\in[n]}B_{i}\log u_{i},

(B.1)

and the negative target function from the Shmyrev convex program is

\Psi(b)=-\sum_{i\in[n],j\in[m]}b_{ij}\log\frac{v_{ij}}{p_{j}}.

(B.2)

We first set up the following convex set:

\mathcal{B}=\left\{b\in\mathcal{M}_{n\times m}(\mathbb{R}_{+}):\sum_{j}b_{ij}=% B_{i}\right\}

(B.3)

To show convergence of the PR dynamics, we first need the following inequalities:

Lemma B.1.

Let $b^{*}=\operatorname{arg\,min}_{b\in\mathcal{B}}\Phi(b)$ . Then $\Phi(b^{*})=\Psi(b^{*})-\sum_{i\in[n]}B_{i}\log B_{i}$ .

Proof.

By KKT optimality constraints of the Eisenberg-Gale convex program, we see that

\frac{B_{i}}{u_{i}^{*}}=\frac{p_{j}^{*}}{v_{ij}},\quad\forall i,j,

(B.4)

which we use to show that $\Phi(b^{*})=\Psi(b^{*})-\sum_{i\in[n]}B_{i}\log B_{i}$ .

$\displaystyle\Phi(b^{*})$	$\displaystyle=-\sum_{i\in[n]}B_{i}\log u_{i}^{}=-\sum_{i\in[n],j\in[m]}b_{ij}% ^{}\log u_{i}^{*}$	(B.5)
	$\displaystyle=-\sum_{i\in[n],j\in[m]}b_{ij}^{}\log B_{i}\frac{v_{ij}}{p_{j}^{% }}=-\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{v_{ij}}{p_{j}^{}}-\sum_{i\in[n% ]}B_{i}\log B_{i}$	(B.6)
	$\displaystyle=\Psi(b^{*})-\sum_{i\in[n]}B_{i}\log B_{i}$	(B.7)

∎

Lemma B.2 (Lemma 19, [53]).

$\forall b\in\mathcal{B},\Phi(b)\leq\Psi(b)-\sum_{i\in[n]}B_{i}\log B_{i}.$

Proof.

We reiterate the proof of Lemma 19 in [53] for the convenience of the reader. By convexity of $-\log$ , we see

$\displaystyle\Phi(b)$	$\displaystyle=-\sum_{i\in[n]}B_{i}\log u_{i}=-\sum_{i\in[n]}B_{i}\log\sum_{j% \in[m]}\frac{b_{ij}}{p_{j}}v_{ij}$	(B.8)
	$\displaystyle=-\sum_{i\in[n]}B_{i}\log\sum_{j\in[m]}\frac{b_{ij}}{B_{i}}\frac{% v_{ij}}{p_{j}}-\sum_{i\in[n]}B_{i}\log B_{i}$	(B.9)
	$\displaystyle\leq-\sum_{i\in[n],j\in[m]}\frac{b_{ij}}{B_{i}}B_{i}\log\frac{v_{% ij}}{p_{j}}-\sum_{i\in[n]}B_{i}\log B_{i}$	(B.10)
	$\displaystyle=\Psi(b)-\sum_{i\in[n]}B_{i}\log B_{i}$	(B.11)

∎

Lemma B.3.

Let $\displaystyle b_{ij}^{\prime}=B_{i}\frac{b_{ij}v_{ij}/p_{j}(b)}{u_{i}(b)}$ . Then $\forall b\in\mathcal{B},\Psi(b^{\prime})\leq\Phi(b)+\sum_{i\in[n]}B_{i}\log B_% {i}.$

Proof.

Let $p_{j}^{\prime}=\sum_{i}b_{ij}^{\prime}$ . By concavity of $\log$ :

$\displaystyle\Psi(b^{\prime})$	$\displaystyle=\sum_{i\in[n],j\in[m]}b_{ij}^{\prime}\log\frac{p_{j}^{\prime}}{v% _{ij}}=\sum_{i\in[n],j\in[m]}B_{i}\frac{b_{ij}v_{ij}/p_{j}}{u_{i}}\log\frac{p_% {j}^{(t+1)}}{v_{ij}}$	(B.12)
	$\displaystyle\leq\sum_{i\in[n]}B_{i}\log\sum_{j\in[m]}\frac{b_{ij}v_{ij}/p_{j}% }{u_{i}}\frac{p_{j}^{\prime}}{v_{ij}}=\sum_{i\in[n]}B_{i}\log\left(\frac{1}{u_% {i}}\sum_{j\in[m]}\frac{b_{ij}p_{j}^{\prime}}{p_{j}}\right)$	(B.13)
	$\displaystyle=\sum_{i\in[n]}B_{i}\log\sum_{j\in[m]}\frac{b_{ij}p_{j}^{\prime}}% {B_{i}p_{j}}-\sum_{i\in[n]}B_{i}\log u_{i}+\sum_{i\in[n]}B_{i}\log B_{i}$	(B.14)
	$\displaystyle\leq\log\sum_{i\in[n],j\in[m]}\frac{b_{ij}p_{j}^{\prime}}{p_{j}^{% (t)}}-\sum_{i\in[n]}B_{i}\log u_{i}+\sum_{i\in[n]}B_{i}\log B_{i}$	(B.15)
	$\displaystyle=\log\sum_{j\in[m]}p_{j}^{\prime}-\sum_{i\in[n]}B_{i}\log u_{i}+% \sum_{i\in[n]}B_{i}\log B_{i}$	(B.16)
	$\displaystyle=\Phi(b)+\sum_{i\in[n]}B_{i}\log B_{i}.$	(B.17)

∎

From the above two lemmas, we gain the monotonically decreasing properties of iteratively updating $b$ via the PR dynamics on the negative target functions of the Eisenberg-Gale and Shmyrev convex programs:

Lemma B.4.

$\forall t\geq 0,\Phi(b^{(t+1)})\leq\Phi(b^{(t)})$ .

Proof.

Apply Lemma B.2 and Lemma B.3 consequently. ∎

Corollary B.4.1 (Lemma 5, [53]).

$\forall t\geq 0,\Psi(b^{(t+1)})\leq\Psi(b^{(t)})$ .

We now use the following lemmas to construct an end-to-end proof of the convergence of the PR dynamics. In a slight abuse of notation, we adapt the definition of KL divergence to matrices such that for $u,v\in\mathcal{M}(\mathbb{R}_{+})_{p\times q}$ , let $D(u\|v):=\sum_{i\in[p],j\in[q]}u_{ij}\log\frac{u_{ij}}{v_{ij}}$ . The following can then be shown:

Lemma B.5.

$\forall t\geq 0,\sum_{t=0}^{T}\left(\Phi(b^{(t)})-\Phi(b^{*})\right)\leq D(b^{% *}\|b^{(0)})$ .

Proof.

Similar by the proof of Theorem 3 of [25], we first lower bound $\Delta_{t}=D(b^{*}\|b^{(t)})-D(b^{*}\|b^{(t+1)})$ as follows:

$\displaystyle\Delta_{t}$	$\displaystyle=D(b^{}\\|b^{(t)})-D(b^{}\\|b^{(t+1)})$	(B.18)
	$\displaystyle=\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{b^{(t+1)}_{ij}}{b^{(t)% }_{ij}}=\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{B_{i}v_{ij}}{p_{j}^{(t)}u_{i% }^{(t)}}$	(B.19)
	$\displaystyle=\sum_{i\in[n],j\in[m]}\left(b_{ij}^{}\log\frac{v_{ij}}{p_{j}^{% }}+b_{ij}^{}\log\frac{p_{j}^{}}{p_{j}^{(t)}}-b_{ij}^{}\log u_{i}^{(t)}+b_{% ij}^{}\log B_{i}\right)$	(B.20)
	$\displaystyle=\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{v_{ij}}{p_{j}^{}}+% \sum_{j\in[m]}p_{j}^{}\log\frac{p_{j}^{}}{p_{j}^{(t)}}-\sum_{i\in[n]}B_{i}% \log u_{i}^{(t)}+\sum_{i\in[n]}B_{i}\log B_{i}$	(B.21)
	$\displaystyle=-\Psi(b^{})+D(p_{j}^{}\\|p_{j}^{(t)})+\Phi(b^{(t)})+\sum_{i\in[% n]}B_{i}\log B_{i}$	(B.22)
	$\displaystyle=D(p_{j}^{}\\|p_{j}^{(t)})+\Phi(b^{(t)})-\Phi(b^{})$	(B.23)
	$\displaystyle\geq\Phi(b^{(t)})-\Phi(b^{*})$	(B.24)

where the second-to-last equality is by Lemma B.1 and the inequality is by the positivity of KL divergence. Taking the telescoping sum of $\Delta_{t}$ , we see that

\sum_{t=0}^{T}\Delta_{t}=\sum_{t=0}^{T}D(b^{*}\|b^{(t)})-D(b^{*}\|b^{(t+1)})=D% (b^{*}\|b^{(0)})-D(b^{*}\|b^{(T+1)})\geq\sum_{t=0}^{T}(\Phi(b^{(t)})-\Phi(b^{*% }))

(B.25)

Hence, we obtain $\sum_{t=0}^{T}\left(\Phi(b^{(t)})-\Phi(b^{*})\right)\leq D(b^{*}\|b^{(0)})$ . ∎

Proposition B.6.

$\forall t\geq 0,\Phi(b^{(T-1)})-\Phi(b^{*})\leq\frac{D(b^{(0)}\|b^{*})}{T}$ .

Proof.

Combining Lemma B.4 and Lemma B.5, we can write

\Phi(b^{(T-1)})-\Phi(b^{*})\leq\frac{1}{T}\sum_{t=0}^{T-1}\Phi(b^{(t)})-\Phi(b% ^{*})\leq\frac{D(b^{*}\|b^{(0)})}{T}.

(B.26)

∎

Corollary B.6.1 (Lemma 3, [53]).

$\forall t\geq 0,\Psi(b^{(t)})-\Psi(b^{*})\leq\frac{D(b^{(0)}\|b^{*})}{T}$ .

Proof.

Apply Lemma B.3 to Proposition B.6. ∎

Lastly, we can upper bound the value $D(b^{*}\|b^{(0)})$ in terms of dimensions $m$ and $n$ given that each buyer initially divides the budget equally between all items such that $b^{(0)}_{ij}=\frac{B_{i}}{m}$ .

Lemma B.7 (Lemma 13, [53]; Theorem 7, [26]).

If $b_{ij}^{(0)}=\frac{B_{i}}{m}$ for all $i$ and $j$ , then $D(b^{*}\|b^{(0)})\leq\log m$ .

Proof.

Evaluating $D(b^{*}\|b^{(0)})$ , we have

D(b^{*}\|b^{(0)})=\sum_{ij}b_{ij}^{*}\log\frac{b_{ij}^{*}}{b_{ij}^{(0)}}=\sum_% {ij}b_{ij}^{*}\log\frac{mb_{ij}^{*}}{B_{i}}=\log m+\sum_{ij}b_{ij}^{*}\log% \frac{b_{ij}^{*}}{B_{i}}\leq\log m

(B.27)

∎

Plugging Lemma B.7 into Proposition B.6 and Corollary B.6.1, we obtain the convergence guarantee of Theorem 2.2.

Appendix C Convergence guarantees for the FPR dynamics

In this section, we prove the convergence guarantee of the faulty proportional response (FPR) dynamics. We first examine the immediate effects of allowing erroneous estimations of $u$ and $p$ in the FPR dynamics. Let $\hat{B}_{i}^{(t)}=\sum_{j\in[m]}\hat{b}_{ij}^{(t)}$ . Note that $\hat{B}_{i}^{(t)}\neq B_{i}$ as the normalization step of constructing $\hat{b}_{ij}$ is erroneous. By the construction of $\hat{b}^{(t)}$ by the FPR dynamics,

\hat{b}_{ij}^{(t)}=\frac{B_{i}}{\tilde{\nu}_{i}^{(t-1)}}v_{ij}\frac{\hat{b}_{% ij}^{(t-1)}}{\tilde{p}_{j}^{(t-1)}},

(C.1)

we can find that

\hat{B}_{i}^{(t)}=\sum_{j\in[m]}\hat{b}_{ij}^{(t)}=\sum_{j\in[m]}\frac{B_{i}}{% \tilde{\nu}_{i}^{(t-1)}}v_{ij}\frac{\hat{b}_{ij}^{(t-1)}}{\tilde{p}_{j}^{(t-1)% }}=\frac{B_{i}}{\tilde{\nu}_{i}^{(t-1)}}\sum_{j\in[m]}v_{ij}\frac{\hat{b}_{ij}% ^{(t-1)}}{\tilde{p}_{j}^{(t-1)}}=B_{i}\frac{\hat{\nu}_{i}^{(t-1)}}{\tilde{\nu}% _{i}^{(t-1)}},

(C.2)

where we can obtain the following inequality by definition of $\tilde{\nu}_{i}$ :

\frac{B_{i}}{1+\epsilon_{\nu}}\leq\hat{B}_{i}^{(t)}\leq\frac{B_{i}}{1-\epsilon% _{\nu}}

(C.3)

By summing $\hat{B_{i}}$ , we find that

\frac{1}{1+\epsilon_{\nu}}\leq\sum_{i\in[n]}\hat{B}_{i}^{(t)}=\sum_{i\in[n],j% \in[m]}\hat{b}_{ij}^{(t)}=\sum_{j\in[m]}\hat{p}_{j}^{(t)}\leq\frac{1}{1-% \epsilon_{\nu}}.

(C.4)

We now prove Theorem 3.1.

See 3.1

Proof.

Similar to the proof of Lemma B.5, we first lower bound $\hat{\Delta}_{t}=\sum_{i\in[n],j\in[m]}b_{ij}^{*}\log\frac{b_{ij}^{*}}{\hat{b}% ^{(t)}_{ij}}-\sum_{i\in[n],j\in[m]}b_{ij}^{*}\log\frac{b_{ij}^{*}}{\hat{b}^{(t% +1)}_{ij}}$ , where we use as follows:

$\displaystyle\hat{\Delta}_{t}$	$\displaystyle=\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{\hat{b}^{(t+1)}_{ij}}{% \hat{b}^{(t)}_{ij}}=\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{B_{i}v_{ij}}{% \tilde{p}_{j}^{(t)}\tilde{\nu}_{i}^{(t)}}$	(C.5)
	$\displaystyle=\sum_{i\in[n],j\in[m]}\left(b_{ij}^{}\log\frac{v_{ij}}{p_{j}^{% }}+b_{ij}^{}\log\frac{p_{j}^{}}{\tilde{p}_{j}^{(t)}}-b_{ij}^{}\log\tilde{% \nu}_{i}^{(t)}+b_{ij}^{}\log B_{i}\right)$	(C.6)
	$\displaystyle=\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{v_{ij}}{p_{j}^{}}+% \sum_{j\in[m]}p_{j}^{}\log\frac{p_{j}^{}}{\tilde{p}_{j}^{(t)}}-\sum_{i\in[n]% }B_{i}\log\tilde{\nu}_{i}^{(t)}+\sum_{i\in[n]}B_{i}\log B_{i}$	(C.7)
	$\displaystyle=-\Phi(b^{})+\sum_{j\in[m]}p_{j}^{}\log\frac{p_{j}^{*}}{\tilde{% p}_{j}^{(t)}}-\sum_{i\in[n]}B_{i}\log\tilde{\nu}_{i}^{(t)}.$	(C.8)

We now lower bound the second and third terms from the above individually as follows. Starting with the second term,

$\displaystyle\sum_{j\in[m]}p_{j}^{}\log\frac{p_{j}^{}}{\tilde{p}_{j}^{(t)}}$	$\displaystyle=\sum_{j\in[m]}p_{j}^{}\log\frac{p_{j}^{}}{\hat{p}_{j}^{(t)}/% \sum_{j^{\prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}}+\sum_{j\in[m]}p_{j}^{}\log% \frac{\hat{p}_{j}^{(t)}}{\tilde{p}_{j}^{(t)}}-\sum_{j\in[m]}p_{j}^{}\log\sum_% {j^{\prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}$	(C.9)
	$\displaystyle=D\left(p_{j}^{}\middle\\|\frac{\hat{p}_{j}^{(t)}}{\sum_{j^{% \prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}}\right)+\sum_{j\in[m]}p_{j}^{}\log% \frac{\hat{p}_{j}^{(t)}}{\tilde{p}_{j}^{(t)}}-\log\sum_{j\in[m]}\hat{p}_{j}^{(% t)}$	(C.10)
	$\displaystyle\geq D\left(p_{j}^{}\middle\\|\frac{\hat{p}_{j}^{(t)}}{\sum_{j^{% \prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}}\right)+\sum_{j\in[m]}p_{j}^{}\log% \frac{1}{1+\epsilon_{p}}-\log\frac{1}{1-\epsilon_{\nu}}$	(C.11)
	$\displaystyle\geq D\left(p_{j}^{*}\middle\\|\frac{\hat{p}_{j}^{(t)}}{\sum_{j^{% \prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}}\right)-\epsilon_{p}-2\epsilon_{\nu}% \geq-\epsilon_{p}-2\epsilon_{\nu}$	(C.12)

Moving on the the third term,

$\displaystyle-\sum_{i\in[n]}B_{i}\log\tilde{\nu}_{i}^{(t)}$	$\displaystyle=-\sum_{i\in[n]}B_{i}\log\hat{\nu}_{i}^{(t)}-\sum_{i\in[n]}B_{i}% \log\frac{\tilde{\nu}_{i}^{(t)}}{\hat{\nu}_{i}^{(t)}}$	(C.13)
	$\displaystyle\geq-\sum_{i\in[n]}B_{i}\log\hat{\nu}_{i}^{(t)}-\sum_{i\in[n]}B_{% i}\log(1+\epsilon_{\nu})$	(C.14)
	$\displaystyle=-\sum_{i\in[n]}B_{i}\log{\sum_{j\in[m]}\frac{\hat{p}_{j}^{(t)}}{% \tilde{p}_{j}^{(t)}}}\frac{v_{ij}\hat{b}_{ij}^{(t)}}{\hat{p}_{j}^{(t)}}-\log(1% +\epsilon_{\nu})$	(C.15)
	$\displaystyle\geq-\sum_{i\in[n]}B_{i}\log\sum_{j\in[m]}\frac{1}{1-\epsilon_{p}% }\frac{v_{ij}\hat{b}_{ij}^{(t)}}{\hat{p}_{j}^{(t)}}-\log(1+\epsilon_{\nu})$	(C.16)
	$\displaystyle=-\sum_{i\in[n]}B_{i}\log\sum_{j\in[m]}\frac{v_{ij}\hat{b}_{ij}^{% (t)}}{\hat{p}_{j}^{(t)}}+\log(1-\epsilon_{p})-\log(1+\epsilon_{\nu})$	(C.17)
	$\displaystyle\geq\Phi(\hat{b}^{(t)})-2\epsilon_{p}-\epsilon_{\nu}$	(C.18)

Hence, in total, we find that

\hat{\Delta}_{t}=\sum_{i\in[n],j\in[m]}b_{ij}^{*}\log\frac{\hat{b}^{(t+1)}_{ij% }}{\hat{b}^{(t)}_{ij}}\geq\Phi(\hat{b}^{(t)})-\Phi(b^{*})-3\epsilon_{p}-3% \epsilon_{\nu}.

(C.19)

Taking the telescoping sum of $\hat{\Delta}_{t}$ , we see that

\sum_{t=0}^{T}\hat{\Delta}_{t}=\sum_{i\in[n],j\in[m]}b_{ij}^{*}\log\frac{b_{ij% }^{*}}{b^{(0)}_{ij}}-\sum_{i\in[n],j\in[m]}b_{ij}^{*}\log\frac{b_{ij}^{*}}{% \hat{b}^{(t+1)}_{ij}}\geq\sum_{t=0}^{T}\left(\Phi(\hat{b}^{(t)})-\Phi(b^{*})-3% \epsilon_{p}-3\epsilon_{\nu}\right).

(C.20)

Taking the upper bound of $\sum_{t=0}^{T}\hat{\Delta}_{t}$ , we obtain

$\displaystyle\sum_{t=0}^{T}\hat{\Delta}_{t}$	$\displaystyle=\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{b_{ij}^{}}{b^{(0)}_{% ij}}-\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{\hat{b}_{ij}^{}}{\hat{b}^{(t+1% )}_{ij}}$	(C.21)
	$\displaystyle\leq\begin{multlined}\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{b_% {ij}^{}}{b^{(0)}_{ij}}-\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\frac{b_{ij}^{}}{% \hat{b}^{(t+1)}_{ij}/\sum_{i^{\prime}\in[n],j^{\prime}\in[m]}\hat{b}^{(t+1)}_{% i^{\prime}j^{\prime}}}\\ +\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\sum_{i^{\prime}\in[n],j^{\prime}\in[m]}% \hat{b}^{(t+1)}_{i^{\prime}j^{\prime}}\end{multlined}\sum_{i\in[n],j\in[m]}b_{% ij}^{}\log\frac{b_{ij}^{}}{b^{(0)}_{ij}}-\sum_{i\in[n],j\in[m]}b_{ij}^{}% \log\frac{b_{ij}^{}}{\hat{b}^{(t+1)}_{ij}/\sum_{i^{\prime}\in[n],j^{\prime}% \in[m]}\hat{b}^{(t+1)}_{i^{\prime}j^{\prime}}}\\ +\sum_{i\in[n],j\in[m]}b_{ij}^{}\log\sum_{i^{\prime}\in[n],j^{\prime}\in[m]}% \hat{b}^{(t+1)}_{i^{\prime}j^{\prime}}$	(C.24)
	$\displaystyle=D(b^{}\\|b^{(0)})-D\left(b^{}\middle\\|\frac{\hat{b}^{(t+1)}}{% \sum_{i^{\prime}\in[n],j^{\prime}\in[m]}\hat{b}^{(t+1)}_{i^{\prime}j^{\prime}}% }\right)+\log\sum_{i^{\prime}\in[n],j^{\prime}\in[m]}\hat{b}^{(t+1)}_{i^{% \prime}j^{\prime}}$	(C.25)
	$\displaystyle\leq D(b^{}\\|b^{(0)})-D\left(b^{}\middle\\|\frac{\hat{b}^{(t+1)}% }{\sum_{i^{\prime}\in[n],j^{\prime}\in[m]}\hat{b}^{(t+1)}_{i^{\prime}j^{\prime% }}}\right)+\log(1+\epsilon_{\nu})$	(C.26)
	$\displaystyle\leq D(b^{*}\\|b^{(0)})+\epsilon_{\nu}$	(C.27)

Hence, we obtain $\sum_{t=0}^{T}\left(\Phi(\hat{b}^{(T)})-\Phi(b^{*})\right)\leq D(b^{*}\|b^{(0)% })+(3T+4)\epsilon_{\nu}+(3T+3)\epsilon_{p}$ . Instead of $T$ , we plug in $T-1$ to obtain

\sum_{t=0}^{T-1}\left(\Phi(\hat{b}^{(t)})-\Phi(b^{*})\right)\leq D(b^{*}\|b^{(% 0)})+(3T+1)\epsilon_{\nu}+3T\epsilon_{p}.

(C.28)

With a simple observation that

T\cdot\min_{t\in[T]}\left(\Phi(\hat{b}^{(t)})-\Phi(b^{*})\right)\leq\sum_{t=0}% ^{T-1}\left(\Phi(\hat{b}^{(t)})-\Phi(b^{*})\right),

(C.29)

we find

\min_{t\in[T]}\left(\Phi(\hat{b}^{(t)})-\Phi(b^{*})\right)\leq\frac{D(b^{*}\|b% ^{(0)})}{T}+4\epsilon_{\nu}+3\epsilon_{p}.

(C.30)

To upper bound $D(b^{*}\|b^{(0)})$ , we use the result of Lemma 13 of [53] and Theorem 7 of [26] as follows:

D(b^{*}\|b^{(0)})=\sum_{ij}b_{ij}^{*}\log\frac{b_{ij}^{*}}{b_{ij}^{(0)}}=\sum_% {ij}b_{ij}^{*}\log\frac{mb_{ij}^{*}}{B_{i}}=\log m+\sum_{ij}b_{ij}^{*}\log% \frac{b_{ij}^{*}}{B_{i}}\leq\log m

(C.31)

where the last inequality is due to $\frac{b_{ij}^{*}}{B_{i}}\leq 1$ .

Then by setting $\epsilon_{\nu}=\frac{1}{8T}$ and $\epsilon_{\nu}=\frac{1}{6T}$ , we obtain

\min_{t\in[T]}\Phi(\hat{b}^{(t)})-\Phi(b^{*})\leq\frac{2\log m}{T}.

(C.32)

∎

Next, we prove Theorem 3.2.

See 3.2

Proof.

We slightly modify the proof of Theorem 3.1, and note that by Equation C.8 and Equation C.12, we have

\displaystyle\hat{\Delta}_{t}

\displaystyle=-\Phi(b^{*})+\sum_{j\in[m]}p_{j}^{*}\log\frac{p_{j}^{*}}{\tilde{% p}_{j}^{(t)}}-\sum_{i\in[n]}B_{i}\log\tilde{\nu}_{i}^{(t)}

\displaystyle\geq-\Phi(b^{*})-\epsilon_{p}-2\epsilon_{\nu}-\sum_{i\in[n]}B_{i}% \log\tilde{\nu}_{i}^{(t)}

(C.33)

Taking the telescoping sum and the upper bound from Equation C.27, we obtain

\sum_{t=0}^{T-1}\left(-\sum_{i\in[n]}B_{i}\log\tilde{\nu}_{i}^{(t)}-\Phi(b^{*}% )\right)\leq D(b^{*}\|b^{(0)})+(2T+1)\epsilon_{\nu}+T\epsilon_{p},

(C.34)

where we can note

\min_{t\in[T]}\left(-\sum_{i\in[n]}B_{i}\log\tilde{\nu}_{i}^{(t)}-\Phi(b^{*})% \right)\leq\frac{D(b^{*}\|b^{(0)})}{T}+3\epsilon_{\nu}+\epsilon_{p}.

(C.35)

Let $t^{*}=\operatorname{arg\,min}_{t\in[T]}\left(-\sum_{i\in[n]}B_{i}\log\tilde{% \nu}_{i}^{(t)}-\Phi(b^{*})\right)$ . Then by Equation C.18, we have the following:

-\sum_{i\in[n]}B_{i}\log\tilde{\nu}_{i}^{(t)}\geq\Phi(\hat{b}^{(t)})-2\epsilon% _{p}-\epsilon_{\nu}

(C.36)

Then we can obtain

\Phi(\hat{b}^{(t^{*})})-\Phi(b^{*})\leq\frac{D(b^{*}\|b^{(0)})}{T}+4\epsilon_{% \nu}+3\epsilon_{p}.

(C.37)

Lastly by setting $\epsilon_{\nu}=\frac{1}{8T}$ and $\epsilon_{\nu}=\frac{1}{6T}$ , we obtain

\Phi(\hat{b}^{(t^{*})})-\Phi(b^{*})\leq\frac{2\log m}{T}.

(C.38)

∎

Appendix D Experimental and implementation details

Our experiments are conducted on a single NVIDIA P100 GPU and written with the PyTorch library [88]. The optimal objective value is approximately computed by taking the results of the $1000$ -th iteration of the PR dynamics.

For the projected gradient descent (PGD) algorithm, our implementation is unlike Gao and Kroer [26], whose task is based on the CEEI scenario where agents are given a unit of fake money and whose end goal is only the allocation. We require information on both the allocation $x$ and price $p$ , hence our algorithm output should be the bids $b$ . Therefore, instead of formulating the problem after the EG objective function, we mirror⁵⁵5Pun intended. the PR dynamics in its equivalence to mirror descent [53] on the Shmyrev objective function and perform PGD on the latter (see Algorithm 2).

Input: Budget

B

, Value

v

Learning rate

\gamma

, Iterations

T

Output: Bids

b

b_{ij}^{(0)}=\frac{B_{i}}{m}

2 for $t=0$ to $T$ do

r_{ij}^{(t)}=b_{ij}^{(t)}-\gamma\cdot(1-\log v_{ij}/p_{j}^{(t)})

// Gradient step

3 for $i=0$ to $n$ do

b_{i,*}^{(t+1)}=\operatorname{Proj}(r_{i,*}^{(t)}\to\{x\in\mathbb{R}_{+}^{n},% \sum_{k}x_{k}=B_{i}\})

// Projection step onto a B_i-simplex

return

b^{(T)}

Algorithm 2 Projected Gradient Descent

We formulate the Shmyrev objective function into the following form to obtain convergence guarantees and the step size:

f(x)=h(Ax)+\langle q,x\rangle

(D.1)

where $x\in\mathbb{R}^{n},A\in\mathcal{M}_{d\times n}(\mathbb{R}),h:\mathbb{R}^{d}\to% \mathbb{R},q\in\mathbb{R}^{n}$ . Considering a flattened vector of the bids $b$ , we note that if

A=n\left\{\begin{pmatrix}1\\ 1\\ \vdots\\ 1\end{pmatrix}\right.\otimes\underbrace{\begin{pmatrix}1&0&\cdots&0\\ 0&1&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&1\end{pmatrix}}_{m},\quad q=\begin{pmatrix}-\log v_{11}\\ -\log v_{12}\\ \vdots\\ -\log v_{mn}\end{pmatrix},\quad h(x)=\sum_{i}x_{i}\log x_{i}

(D.2)

then $f=\Psi$ . Then by Theorem 3 of [26], by setting a learning rate of $\gamma=1/L\|A\|^{2}$ , where $L=1/\min_{j,t}p_{j}^{(t)}$ , we get linear convergence. Note that $\|A\|^{2}=n$ . Gao and Kroer [26] further provide a line search procedure to set the constant multiplier in the learning rate as well as provide sharper convergence guarantees, but as we only run for $16$ iterations, we do not perform the line search and fix the learning rate to the initial learning rate that Gao and Kroer [26] use in their empirical studies, which is $1000/L\|A\|^{2}$ .

For amplitude estimation, we set $M=\sqrt{T\sqrt{n}}/16=32$ . We scale down $M$ by the constant factor of $16$ to save memory consumption on the GPU, as we simulate amplitude estimation by computing the full probability distribution over $[M]$ . We compensate for the loss in accuracy of the estimation by employing the median-of-means estimator [54], where we take the median of $3$ estimators constructed from the mean of $7$ samples from the amplitude estimation subroutine. We also assume that the maximum finding algorithm is always successful in our algorithm.

$\displaystyle\sum_{j\in[m]}p_{j}^{}\log\frac{p_{j}^{}}{\tilde{p}_{j}^{(t)}}$	$\displaystyle=\sum_{j\in[m]}p_{j}^{}\log\frac{p_{j}^{}}{\hat{p}_{j}^{(t)}/% \sum_{j^{\prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}}+\sum_{j\in[m]}p_{j}^{}\log% \frac{\hat{p}_{j}^{(t)}}{\tilde{p}_{j}^{(t)}}-\sum_{j\in[m]}p_{j}^{}\log\sum_% {j^{\prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}$	(C.9)
	$\displaystyle=D\left(p_{j}^{}\middle\\|\frac{\hat{p}_{j}^{(t)}}{\sum_{j^{% \prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}}\right)+\sum_{j\in[m]}p_{j}^{}\log% \frac{\hat{p}_{j}^{(t)}}{\tilde{p}_{j}^{(t)}}-\log\sum_{j\in[m]}\hat{p}_{j}^{(% t)}$	(C.10)
	$\displaystyle\geq D\left(p_{j}^{}\middle\\|\frac{\hat{p}_{j}^{(t)}}{\sum_{j^{% \prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}}\right)+\sum_{j\in[m]}p_{j}^{}\log% \frac{1}{1+\epsilon_{p}}-\log\frac{1}{1-\epsilon_{\nu}}$	(C.11)
	$\displaystyle\geq D\left(p_{j}^{*}\middle\\|\frac{\hat{p}_{j}^{(t)}}{\sum_{j^{% \prime}\in[m]}\hat{p}_{j^{\prime}}^{(t)}}\right)-\epsilon_{p}-2\epsilon_{\nu}% \geq-\epsilon_{p}-2\epsilon_{\nu}$	(C.12)