Fast and Interpretable Dynamics for Fisher Markets via Block-Coordinate Updates

TL;DR

This paper introduces scalable block-coordinate first-order methods for Fisher market equilibrium computation, offering interpretable dynamics and fast convergence suitable for large-scale problems.

Contribution

It develops novel block-coordinate first-order algorithms with interpretations as market dynamics, improving scalability and convergence analysis for Fisher market equilibria.

Findings

Methods have provably fast convergence rates.

Algorithms are interpretable as market dynamics.

Applicable to large-scale Fisher market problems.

Abstract

We consider the problem of large-scale Fisher market equilibrium computation through scalable first-order optimization methods. It is well-known that market equilibria can be captured using structured convex programs such as the Eisenberg-Gale and Shmyrev convex programs. Highly performant deterministic full-gradient first-order methods have been developed for these programs. In this paper, we develop new block-coordinate first-order methods for computing Fisher market equilibria, and show that these methods have interpretations as t\^atonnement-style or proportional response-style dynamics where either buyers or items show up one at a time. We reformulate these convex programs and solve them using proximal block coordinate descent methods, a class of methods that update only a small number of coordinates of the decision variable in each iteration. Leveraging recent advances in the convergence analysis of these methods and structures of the equilibrium-capturing convex programs, we establish fast convergence rates of these methods.

Equations

$${\mathcal{D}}_i(p)=\arg \underset{x_i}{\max }\{u_i=⟨v_i,x_i⟩:⟨p,x_i⟩\le B_i\}.$$

$$\begin{array}{rl}\underset{x}{\max }& {\sum }_i{B}_i\log u_i\\ \text{s.t.}& u_i\le ⟨v_i,x_i⟩\forall i\\ & {\sum }_i{x}_{ij}\le 1\forall j\\ & x\ge 0.\end{array}$$

$$\underset{x\in {\mathbb{R}}^{m\times n}}{\min }F(x):=f(x)+\psi (x)=f(x)+{\sum }_j{\psi }_j(x_{\cdot j})$$

$$\arg \underset{y_{\cdot j}}{\min }⟨{\nabla }_{\cdot j}f(x),y_{\cdot j}-x_{\cdot j}⟩+\frac{1}{{\eta }_j}\Vert y_{\cdot j}-x_{\cdot j}{\Vert }^2+{\psi }_j(y_{\cdot j}).$$

$$T_{j}(x)=\text{Proj}_{\Delta_{n}}\Big{(}x_{\cdot j}-\eta_{j^{\prime}}\nabla_{\cdot j}f(x)\Big{)},$$

$$L_j=\underset{i\in [n]}{\max }\frac{B_i{v}_{ij}^2}{{\underline{u}}_i^2}.$$

$$\tilde{f}(y)\le \tilde{f}(x)+⟨{\nabla }_{\cdot j}\tilde{f}(x),y-x⟩+\frac{L_j}{2}{\Vert y-x\Vert }^2.$$

$${\eta }_{j_k}\Vert {\nabla }_{\cdot j_k}\tilde{f}(x^k)-{\nabla }_{\cdot j_k}\tilde{f}(x^{k-1})\Vert \le \Vert x_{\cdot j_k}^k-x_{\cdot j_k}^{k-1}\Vert .$$

$$-2L\underset{y}{\min }⟨\nabla \tilde{f}(x),y-x⟩+\frac{L}{2}{\Vert y-x\Vert }^2+\psi (y)-\psi (x)$$

$$\frac{1}{2}{\mathcal{D}}_{\psi }(x,L)\ge \min \left\{\frac{\mu }{{\theta }^2(A,C)},L\right\}(F(x)-F^{\ast })$$

$$\mathbf{E}\big{[}F(x^{k+1})\big{]}-F^{*}\leq(1-\rho)^{k}\left(F(x^{0})-F^{*}\right),$$

$$x_{\cdot j}^k\leftarrow {\text{Proj}}_{{\Delta }^n}\left(x_{\cdot j}^{k-1}-{\eta }_j{g}^{k-1}\right).$$

$${\sum }_i(y_i-t{)}^{+}=1$$

$${\sum }_i(x_{ij}^{k-1}-{\eta }_j{g}_i^{k-1}-t{)}^{+}=1.$$

$${\sum }_i{D}_i^k(p_j)=1,D_i^k(p_j)=(x_{ij}^{k-1}-{\eta }_j{g}_i^{k-1}-{\eta }_j{p}_j{)}^{+}.$$

$$\begin{array}{rl}\underset{b}{\max }& {\sum }_{i,j}{a}_{ij}{b}_{ij}-{\sum }_j{p}_j(b)\log p_j(b)\\ \text{s.t.}& {\sum }_j{b}_{ij}=B_i\forall i\\ & b\ge 0.\end{array}$$

$$\underset{b\in {\mathbb{R}}^{m\times n}}{\min }\Phi (b):=\phi (b)+r(b)=\phi (b)+{\sum }_i{r}_i(b_i)$$

$$\arg \underset{a}{\min }⟨{\nabla }_i\phi (b),a-b_i⟩+\frac{1}{{\alpha }_i}{D}_{\mathrm{KL}}(a,b_i)+r_i(a)$$

$$b_{ij}^{+}=\frac{1}{Z_i}{b}_{ij}{\left(\frac{v_{ij}}{p_j}\right)}^{{\alpha }_i}\forall j=1,2,\dots ,m$$

$$\mathbf{E}\big{[}\Phi(b^{k})\big{]}-\Phi^{*}\leq\frac{n}{n+k}\Big{(}\Phi^{*}-\Phi(b^{0})+D(b^{*},b^{0})\Big{)}$$

$$\mathbf{E}\left[\Phi (b^k)-\Phi (b^{\ast })\right]$$

$$\le \frac{n}{n+k}\left(\Phi (b^{\ast })-\Phi (b^0)+\sum _i\frac{1}{{\alpha }_i^0}D(b_i^{\ast },b^0)\right)$$

$$+\frac{1}{k}\sum _{l=1}^k\mathbf{E}\left[\sum _i\frac{1}{{\alpha }_i^l}D(b_i^{\ast },b^l)-\sum _i\frac{1}{{\alpha }_i^{l-1}}D(b_i^{\ast },b^l)\right].$$

$$\phi (b^k)-{\phi }^{\ast }\le \frac{1}{{\rho }^{-}}\cdot \frac{D(b^{\ast },b^0)}{k}$$

$$L_j=\underset{i\in [n]}{\max }\frac{B_i{v}_{ij}{\underline{x}}_{ij}^{\rho -2}}{{\underline{u}}_i(\rho )}\text{for all }j\in [m]$$

$$\mathbf{E}\big{[}F(x^{k})\big{]}-F^{*}\leq\Big{(}1-\frac{\mu(L)}{m}\Big{)}^{k}\left(F(x^{0})-F^{*}\right)$$

$$k\ge \frac{2\log \frac{\sqrt{8D(b^{\ast },b^0)}{W}^{\frac{1}{1-\rho }}}{ϵ}}{\log \frac{n}{n-1+\rho }},W=\frac{n}{{\min }_{v_{ij}>0}{v}_{ij}\cdot {\min }_i{B}_i},$$

$${\nabla }_{ij}^2\tilde{f}(x)=\{\begin{array}{ll}{B}_i{v}_{ij}^2/u_i^2(x)& u_i(x)\ge {\underline{u}}_i\\ B_i{v}_{ij}^2/{\underline{u}}_i^2& u_i(x)<{\underline{u}}_i.\end{array}$$

$$\left|{\nabla }_{ij}\tilde{f}(y)-{\nabla }_{ij}\tilde{f}(x)\right|=\left|{\nabla }_{ij}\tilde{f}(y_i)-{\nabla }_{ij}\tilde{f}(x_i)\right|\le \left|\frac{B_i{v}_{ij}^2}{{\underline{u}}_i^2}\right||y_{ij}-x_{ij}|=L_j\Vert y_i-x_i\Vert ,$$

$${\nabla }_{ij}^2\tilde{f}(x)\le \frac{B_i{v}_{ij}^2}{{\underline{u}}_i^2}\le \underset{i\in [n]}{\max }\frac{B_i{v}_{ij}^2}{{\underline{u}}_i^2}=L_j.$$

$${\Vert {\nabla }_{\cdot j}\tilde{f}(y)-{\nabla }_{\cdot j}\tilde{f}(x)\Vert }^2=\sum _i{\left|{\nabla }_{ij}\tilde{f}(y)-{\nabla }_{ij}\tilde{f}(x)\right|}^2\le L_j^2\sum _i\Vert y_i-x_i{\Vert }^2=L_j^2\Vert y-x{\Vert }^2=L_j^2\Vert y_{\cdot j}-x_{\cdot j}{\Vert }^2,$$

Full text

1 Introduction

In a market equilibrium (ME) a set of items is allocated to a set of buyers via a set of prices for the items and an allocation of items to buyers such that each buyer spends their budget optimally, and all items are fully allocated. Due to its rich structural properties and strong fairness and efficiency guarantees, ME has long been used to develop fair division and online resource allocation mechanisms (Gao and Kroer 2021; Aziz and Ye 2014; Barman, Krishnamurthy, and Vaish 2018; Arnsperger 1994). Market model and corresponding equilibrium computation algorithms have been central research topics in market design and related areas in economics, computer science and operations research with practical impacts (Scarf et al. 1967; Kantorovich 1975; Othman, Sandholm, and Budish 2010; Daskalakis, Goldberg, and Papadimitriou 2009; Cole et al. 2017; Kroer et al. 2019). More specifically, many works in market design rely on the assumption that ME can be computed efficiently for large-scale market instances. For example, the well-known fair division mechanism without money competitive equilibrium from equal incomes (CEEI) requires computing a ME of a Fisher market under uniform buyer budgets (Varian et al. 1974). Recent work has also established close connections between market equilibria and important solution concepts in the context of large-scale Internet markets, such as pacing equilibria in repeated auctions (Conitzer et al. 2018, 2019; Kroer and Stier-Moses 2022). Motivated by the classical and emerging applications described above, we are interested in developing efficient equilibrium computation algorithms for large-scale market instances. In general, computing a ME is a hard problem (Chen and Teng 2009; Vazirani and Yannakakis 2011; Othman, Papadimitriou, and Rubinstein 2016). However, for the case of Fisher markets and certain classes of utility functions, efficient algorithms are known (Devanur et al. 2008; Zhang 2011; Gao and Kroer 2020), often based on solving a specific convex program—whose solutions are ME and vice versa—using an optimization algorithm. In this paper, we focus on two well-known such convex programs, namely, the Eisenberg-Gale (EG) (Eisenberg and Gale 1959; Eisenberg 1961) and Shmyrev convex programs (Shmyrev 2009; Cole et al. 2017).

Most existing equilibrium computation literature studies the case of a static market where all buyers and items are present in every time step of the equilibrium computation process. In contrast, we study a setting where only a random subset of buyer-item pairs show up at each time step. Such a setting is well-motivated from a computational perspective, since stochastic methods are typically more efficient for extremely large problems. Secondly, our model allows us to model new types of market dynamics. We make use of recent advances in stochastic first-order optimization, more specifically, block-coordinate-type methods, to design new equilibrium computation algorithms for this setting. The resulting equilibrium computation algorithms have strong convergence guarantees and consistently outperform deterministic full-information algorithms in numerical experiments. In addition, many of the optimization steps not only give efficient update formulas for market iterates, but also translate to interpretable market dynamics.

Summary of contribution.

We propose two stochastic algorithms for computing large-scale ME: (proximal) block-coordinate descent on EG (BCDEG) and block-coordinate proportional response (BCPR). These algorithms are derived by applying stochastic block-coordinate-type algorithms on reformulated equilibrium-capturing convex programs. More specifically, BCDEG is based on (proximal) stochastic block coordinate descent (BCD) and BCPR is based on a non-Euclidean (Bregman) version of BCD. We show that these algorithms enjoy attractive theoretical convergence guarantees, and discuss important details for efficient implementation in practice. Furthermore, we show that the Euclidean projection onto the simplex in BCDEG (and other preexisting projected-gradient-type methods) has a tâtonnement-style interpretation. We then demonstrate the practical efficiency of our algorithms via extensive numerical experiments on synthetic and real market instances, where we find that our algorithms are substantially faster than existing state-of-the-art methods such as proportional response dynamics.

Preliminaries and notation.

Unless otherwise stated, we consider a linear Fisher market with $n$ buyers and $m$ items. We use $i\in[n]$ to denote a buyer and $j\in[m]$ to denote an item. Each buyer $i$ has a budget $B_{i}>0$ and each item has supply one. An allocation (or bundle) for buyer $i$ is a vector $x_{i}\in\mathbb{R}^{m}_{+}$ specifying how much buyer $i$ gets of each item $j$. Given $x_{i}$, buyer $i$ gets utility $\langle v_{i},x_{i}\rangle=\sum_{j}v_{ij}x_{ij}$, where $v_{i}\in\mathbb{R}^{m}_{+}$ is their valuation vector. Given prices $p\in\mathbb{R}^{m}_{+}$ (i.e., price of item $j$ is $p_{j}$) on all items, buyer $i$ pays $\langle p,x_{i}\rangle=\sum_{j}p_{j}x_{ij}$ for $x_{i}$. Given prices $p$ and budget $B_{i}$, a bundle $x_{i}$ is budget feasible for buyer $i$ if $\langle p,x_{i}\rangle\leq B_{i}$. We also use $x_{\cdot j}\in\mathbb{R}^{n}_{+}$ to denote a vector of amounts of item $j$ allocated to all buyers. The demand set of buyer $i$ is the set of budget-feasible utility-maximizing allocations:

[TABLE]

A market equilibrium (ME) is an allocation-price pair $(x^{*},p^{*})$ such that $x^{*}_{i}\in\mathcal{D}_{i}(p^{*})$ for all $i$ and $\sum_{i}x^{*}_{ij}\leq 1$ for all $j$, with equality if $p^{*}_{j}>0$.

2 Related Work

Since the seminal works by Eisenberg and Gale (Eisenberg and Gale 1959; Eisenberg 1961), there has been an extensive literature on equilibrium computation for Fisher markets, often based on convex optimization characterizations of ME (Devanur et al. 2008; Zhang 2011; Birnbaum, Devanur, and Xiao 2011; Cole et al. 2017; Gao and Kroer 2020; Garg and Kapoor 2006). Equilibrium computation algorithms for more general market models with additional constraints—such as indivisible items and restrictions on the set of permissible bundles for each buyer—have also been extensively studied (Othman, Sandholm, and Budish 2010; Budish et al. 2016); these algorithms are often based on approximation algorithms, mixed-integer programming formulations, and local search heuristics.

In Gao and Kroer (2020), the authors considered three (deterministic) first-order optimization methods, namely projected gradient (PG), Frank-Wolfe (FW) and mirror descent (MD) for solving convex programs capturing ME. To the best of our knowledge there are no existing results on block-coordinate methods for ME. In the optimization literature there is an extensive and ongoing literature on new block-coordinate-type algorithms and their analysis (see, e.g., Tseng (2001); Wright (2015); Beck and Tetruashvili (2013); Hanzely and Richtárik (2019, 2021); Liu and Wright (2015); Nesterov (2012); Richtárik and Takáč (2014); Gao et al. (2020); Attouch, Bolte, and Svaiter (2013); Zhang (2020); Reddi et al. (2016)). As mentioned previously, our BCDEG algorithm is based on proximal block-coordinate descent (PBCD) applied to EG. Linear convergence of the mean-square error of the last iterate of PBCD for nonsmooth, finite-sum, composite optimization (with an objective function of the form $F(x)=\sum_{i}f_{i}(x)+\psi(x)$) has been established under different error bound conditions (Richtárik and Takáč 2014; Karimi, Nutini, and Schmidt 2016; Reddi et al. 2016). For BCPR, we adopt the analysis of a recently proposed non-Euclidean (Bregman) PBCD (Gao et al. 2020), which in turn made use of the convergence theory developed in Bauschke, Bolte, and Teboulle (2017).

3 Block Coordinate Descent Algorithm for the EG Program

(Proximal) Block coordinate descent methods (BCD) are often used to solve problems whose objective function consists of a smooth part and a (potentially) nonsmooth part. The second part is typically block-coordinate-wise separable. BCD algorithms update only a small block of coordinates at each iteration, which makes each iteration much cheaper than for deterministic methods, and this enables scaling to very large instances. The EG convex program, which captures market equilibria, can be written in a form amenable to proximal block coordinate descent method. The EG convex program for buyers with linear utility functions is

[TABLE]

Any optimal solution $x^{*}$ to (EG) and the (unique) optimal Lagrange multipliers $p^{*}\in\mathbb{R}^{m}_{+}$ associated with the constraints $\sum_{i}x_{ij}\leq 1$, $j\in[m]$ forms a market equilibrium. In fact, this holds more generally if the utility functions are concave, continuous, nonnegative, and homogeneous with degree $1$.

To apply BCD, note that (EG) is of the following form:

[TABLE]

where $f(x)=-\sum\nolimits_{i}{B_{i}\log{\langle v_{i},x_{i}\rangle}}$ and $\psi_{j}(x_{\cdot j})=0$ if $\sum_{i}x_{ij}\leq 1,x_{\cdot j}\geq 0$ and $+\infty$ otherwise.

Thus, the nonsmooth term $\psi(x)$ decomposes along the items $j$, and we can therefore treat $x$ as a set of blocks: each item $j$ has a block of allocation variables corresponding to how much of item $j$ is given to each buyer. We use $x_{\cdot j}=(x_{1j},\ldots,x_{nj})$ to denote the $j$’th block. Given the full gradient of $f$ at $x$ as $\nabla f(x)$, we will also need the partial gradient w.r.t. the $j$th block $x_{\cdot j}$, which we denote as $\nabla_{\cdot j}f(x)$, that is, $\nabla_{\cdot j}f(x)=\left(\nabla_{1j}f(x),\ldots,\nabla_{nj}f(x)\right)$.

In each iteration $t$ of the BCD method, we first choose an index $j^{\prime}\in[m]$ at random, with a corresponding stepsize $\eta_{j^{\prime}}$. Then, the next iterate $x^{+}$ is generated from $x$ via $x_{\cdot j}^{+}=T_{j}(x)$ if $j=j^{\prime}$ and $x_{\cdot j}$ otherwise, where $T_{j}(x)$ equals

[TABLE]

The above proximal mapping is equivalent to

[TABLE]

so we can generate $T_{j}(x)\in\mathbb{R}^{n}$ via a single projection onto the $n$-dimensional simplex. Since the projection is the most computationally expensive part, our method ends up being cheaper than full projected gradient descent by a factor of $m$.

To make sure $\nabla_{\cdot j}f(x)$ exists, we need to bound buyer utilities away from zero at every iteration. To that end, let $\underline{u}_{i}=\langle v_{i},(B_{i}/\sum_{l=1}^{n}B_{l})\vec{1}\rangle$ be the utility of the proportional allocation. Then we perform “quadratic extrapolation” where we replace the objective function $f(x)$ with the function $\tilde{f}(x)=\sum_{i}\tilde{f}_{i}(x_{i})=\sum_{i}\tilde{g}_{i}(u_{i})$ where $\tilde{g}_{i}(u_{i})=-B_{i}\log{u_{i}}$ if $u_{i}\geq\underline{u}_{i}$ and $\log{\underline{u}_{i}}-\frac{B_{i}}{\underline{u}_{i}}\left(u_{i}-\underline{u}_{i}\right)+\frac{B_{i}}{2\underline{u}_{i}^{2}}\left(u_{i}-\underline{u}_{i}\right)^{2}$ otherwise. For (EG), replacing $f$ with $\tilde{f}$ does not affect any optimal solution (Gao and Kroer 2020, Lemma 1).

We will also need the following Lipschitz bound which ensures that the iterates are descent steps, meaning that the expected objective value is non-increasing as long as the stepsizes are not too large. The upper bound on the allowed stepsizes (that ensure descent iterates) for a specific block is governed by the Lipschitz constant w.r.t. that block of coordinates. More details and proofs can be found in Appendix A.

Lemma 1**.**

For any $j\in[m]$, let

[TABLE]

Then, for all $x,y\in\mathcal{X}$ such that $x,y$ differ only in the $j$th block, we have

[TABLE]

Let $L$ be the “global” Lipschitz constant which determines the stepsize of (full) gradient descent. Let $L_{\max}=\max_{j}L_{j}$. Generally, it is easy to see that $1\leq L/L_{\max}\leq n$, but we can get a stronger bound than this using the gradient $(\nabla f)_{ij}=B_{i}v_{ij}/u_{i}$. When only the variables $x_{\cdot j}$ in block $j$ change, we have that $L_{j}$ is bounded by the maximal diagonal value of the $j$’th-block (sub) Hessian matrix. In contrast, for $L$ it depends on the maximal trace of the $i$’th-block (sub) Hessian matrix over all buyers $i$. From a market perspective this can be interpreted as follows: when we only adjust one block, each buyer’s utility fluctuates based on that one item. On the contrary, for full gradient methods, every item contributes to the change of $u_{i}$. This yields that the ratio of $L_{j}/L$ is roughly $\max_{i}v_{ij}^{2}/\max_{i}\lVert v_{i}\rVert^{2}$.

Algorithm 1 states the (proximal) BCD algorithm for the EG convex program. Each iteration only requires a single projection onto an $n$-dimensional simplex, as opposed to $m$ projections for the full projected gradient method. Moreover, we show in Appendix A that for linear utilities we can further reduce the computational cost per iteration of Algorithm 1 when the valuation matrix $v_{ij}$ is sparse.

Line search strategy.

As mentioned in previous work such as Richtárik and Takáč (2014), line search is often very helpful for BCD. If it can be performed cheaply, then it can greatly speed up numerical convergence. We show later that this occurs for our setting.

We incorporate line search in Algorithm 1 with a (coordinate-wise) $l_{2}$ smoothness condition. The line search modifies Algorithm 1 as follows: after computing $x_{\cdot j_{k}}^{k}$, we check whether

[TABLE]

If this check succeeds then we increase the stepsize by a small multiplicative factor and go to the next iteration. If it fails then we decrease by a small multiplicative factor and redo the calculation of $x_{\cdot j_{k}}^{k}$. This line search algorithm can be implemented in $O(n)$ cost per iteration, whereas a full gradient method requires $O(nm)$ time. A full specification of BCDEG with line search (BCDEG-LS) is given in Appendix A.

Convergence analysis.

Next we establish the linear convergence of Algorithm 1 under reasonably-large (fixed) stepsizes, as well as for the line search variant.

Following prior literature on the linear convergence of first-order methods for structured convex optimization problems under “relaxed strong convexity” conditions, we show that BCDEG generates iterates that have linear convergence of the expected objective value. For more details on these relaxed sufficient conditions that ensure linear convergence of first-order methods, see Karimi, Nutini, and Schmidt (2016) and references therein.

Gao and Kroer (2020) showed that (EG) and other equilibrium-capturing convex programs can be reformulated to satisfy these conditions. Hence, running first-order methods on these convex programs yields linearly-convergent equilibrium computation algorithms. Similar to the proof of Gao and Kroer (2020, Theorem 2) and Karimi, Nutini, and Schmidt (2016, (39)), we first establish a Proximal-PŁ inequality.

Lemma 2**.**

For any feasible $x$ and any $L>0$, define $\mathcal{D}_{\psi}(x,L)$ to be

[TABLE]

then we have the inequality

[TABLE]

where $\theta(A,C)$ is the Hoffman constant of the polyhedral set of optimal solutions of (EG), which is characterized by matrices $A$ and $C$ where $C$ is a matrix capturing optimality conditions, and $F^{*}$ is the optimal value.

In previous inequalities, they essentially showed that for any $L\geq\lambda={\mu}/{\theta^{2}(A,C)}>0$, $\frac{1}{2}\mathcal{D}_{\psi}(x,L)\geq\lambda(F(x)-F^{*})$. However, here we generate a Proximal-PŁ inequality giving a lower bound on $\frac{1}{2}\mathcal{D}_{\psi}(x,L)$ for any $L>0$.

Then, combining Lemma 2 with Richtárik and Takáč (2014, Lemmas 2 & 3), we can establish the main convergence theorem for BCDEG.

Theorem 1**.**

Given an initial iterate $x^{0}$ and stepsizes $\eta_{j}^{0}=1/L_{j},\,\forall\;j$ satisfying (4), let $x^{k}$ be the random iterates generated by Algorithm 1. Then,

[TABLE]

where $\rho=\min\left\{\frac{\mu}{mL_{\max}\theta^{2}(A,C)},\frac{1}{m}\right\}$ and $L_{\max}=\max_{j}L_{j}$.

Karimi, Nutini, and Schmidt (2016) also develop a block-coordinate method that applies to EG. Unlike their result, our result admits larger stepsizes that can vary per block, as well as a line search strategy, which is helpful for practical performance as we show in Section 8.

4 Economic Interpretation of Projected Gradient Steps

As Goktas, Viqueira, and Greenwald (2021) argue, one drawback of computing market equilibrium via projected gradient methods such as BCDEG is that these methods do not give a natural interpretation as market dynamics. To address this deficiency, in this section we show that Algorithm 1, and projected gradient descent more generally, can be interpreted as distributed pricing dynamics that balance supply and demand.

The projection step in Algorithm 1, for an individual buyer $i$ and a chosen item $j$, is as follows (where we use $j$ for $j_{k}$ and drop the time index for brevity):

[TABLE]

As is well-known, the projection of a vector $y\in\mathbb{R}^{n}$ onto the simplex $\Delta^{n}=\left\{x\in\mathbb{R}^{n}_{+}:\sum_{i}x_{i}=1\right\}$ can be found using an $O(n\log n)$ algorithm (the earliest discovery that we know of is Held, Wolfe, and Crowder (1974); see Appendix B for a discussion of more recent work on simplex projection). The key step is to find the (unique) number $t$ such that

[TABLE]

and compute the solution as $x=(y-t\cdot\mathbf{1})^{+}$ (component-wise). This can be done with a simple one-pass algorithm if the $y_{i}$ are sorted. In fact, the number $t$ corresponds to the (unique) optimal Lagrange multiplier of the constraint $\sum_{i}x_{i}=1$ in the KKT conditions.

In the projection step in Algorithm 1, (9) has the form

[TABLE]

Recall that we have $g^{k-1}<0$. Note that the left-hand side of (10) is non-increasing in $t$ and strictly decreasing around the solution (since some terms on the left must be positive for the sum to be $1$). Furthermore, setting $t=0$ gives a lower bound of $\sum_{i}(x_{ij}^{k-1}-\eta_{j}g^{k-1}_{i}-t)^{+}>\sum_{i}x^{k-1}_{ij}=1$. Hence, the unique solution $t^{*}$ must be positive. Now we rewrite $t$ as $\eta_{j}p_{j}$ for some “price” $p_{j}$. Then, (10) can be written as

[TABLE]

In other words, the projection step is equivalent to finding $p_{j}^{k}$ that solves (11). Here, $D^{k}_{i}$ can be viewed as the linear demand function of buyer $i$ at time $k$ given a prior allocation $x^{k-1}_{ij}$. The solution $p_{j}^{k}$ can be seen as a market-clearing price, since $\sum_{i}D^{k}_{i}(p_{j}^{k})$ exactly equals the unit supply of item $j$. After setting the price, the updated allocations can be computed easily just as in the simplex projection algorithm, that is, $x^{k}_{ij}=\left(x^{k-1}_{ij}-\eta_{j}g^{k-1}-\eta_{j}p^{k}_{j}\right)^{+}$ for all $i$. Equivalently, the new allocations are given by the current linear demand function: $x^{k}_{ij}=D^{k}_{i}(p^{k}_{j})$.

Summarizing the above, we can recast Algorithm 1 into the following dynamic pricing steps. At each time $k$, the following events occur.

• •

An item $j$ is sampled uniformly at random.

• •

For each buyer $i$, her demand function becomes $D_{i}^{k}(p_{j})\leftarrow(x^{k-1}_{ij}-\eta_{j}g^{k-1}-\eta_{j}p_{j})^{+}$.

• •

Find the unique price $p_{j}^{k}$ such that $\sum_{i}D^{k}_{i}(p^{k}_{j})=1$.

• •

Each buyer chooses their new allocation of item $j$ via $x^{k}_{ij}=D^{k}_{i}(p^{k}_{j})$.

To get some intuition for the linear demand function, note that when $u_{i}\geq\underline{u}_{i}$, we have $g^{k-1}_{i}=-{B_{i}v_{ij}}/{u^{k-1}_{i}}$, and therefore it holds that $D^{k}_{i}\left({B_{i}v_{ij}}/{u^{k-1}_{i}}\right)=x^{k-1}_{ij}$. In other words, the current demand of buyer $i$ is exactly the previous-round allocation $x^{k-1}_{ij}$ if the price of item $j$ is $(B_{i}/u^{k-1}_{i})v_{ij}$. This can be interpreted in terms familiar from the solution of EG: let $\beta^{k-1}_{i}=B_{i}/u^{k-1}_{i}$ be the utility price at time $k-1$ for buyer $i$, then we get that after seeing prices $p_{j}^{k}$, buyer $i$ increases their allocation of goods that beat their current utility price, and decreases their allocation on goods that are worse than their current utility price. The stepsize $\eta_{j}$ denotes buyer $i$’s responsiveness to price changes.

The fact that the prices are set in a way that equates the supply and demand is reminiscent of tâtonnement-style price setting. The difference here is that the buyers are the ones who slowly adapt to the changing environment, while the prices are set in order to achieve exact market clearing under the current buyer demand functions.

5 Relative Block Coordinate Descent Algorithm for PR dynamics

Birnbaum, Devanur, and Xiao (2011) showed that the Proportional Response dynamics (PR) for linear buyer utilities can be derived by applying mirror descent (MD) with the KL divergence on the Shmyrev convex program formulated in buyers’ bids (Shmyrev 2009; Cole et al. 2017). The authors derived an $O(1/k)$ last-iterate convergence guarantee of PR (MD) by exploiting a “relative smoothness” condition of the Shmyrev convex program. This has later been generalized and led to faster MD-type algorithms for more general relatively smooth problems (Hanzely and Richtárik 2021; Lu, Freund, and Nesterov 2018; Gao et al. 2020). In this section, we propose a randomized extension of PR dynamics, which we call block coordinate proportional response (BCPR). BCPR is based on a recent stochastic mirror descent algorithm (Gao et al. 2020). We provide stepsize rules and show that each iteration involves only a single buyer and can be performed in $O(m)$ time.

Let $b\in\mathbb{R}^{n\times m}$ denote the matrix of all buyers’ bids on all items and $p_{j}(b):=\sum_{i}b_{ij}$ denote the price of item $j$ given bids $b$. Denote $a_{ij}=\log v_{ij}$ if $v_{ij}>0$ and [math] otherwise. The Shmyrev convex program is

[TABLE]

It is known that an optimal solution $(b^{*},p^{*})$ of (S) gives equilibrium prices $p^{*}_{j}$. Corresponding equilibrium allocations can be constructed via $x_{ij}^{*}=b_{ij}^{*}/p_{j}^{*}$ for all $i,j$. (S) can be rewritten as minimization of a smooth finite-sum convex function with a potentially nonsmooth convex separable regularizer:

[TABLE]

where $\varphi(b)=-\sum\nolimits_{i,j}{b_{ij}\log\left(\frac{v_{ij}}{p_{j}(b)}\right)}$ and $r_{i}(b_{i})=0$ if $\sum_{j}b_{ij}=B_{i},\,b_{i}\geq 0$ and $+\infty$ otherwise.

Now we introduce the relative randomized block coordinate descent (RBCD) method for (12). We use the KL divergence as the Bregman distance in the proximal update of $b$. Let $D_{\rm KL}(q^{1},q^{2})=\sum_{i}{q^{1}_{i}\log{(q^{1}_{i}/q^{2}_{i})}}$ denote the KL divergence between $q^{1}$ and $q^{2}$ (assuming $\sum_{i}q^{1}_{i}=\sum_{i}q^{2}_{i}$ and $q^{1}_{i},q^{2}_{i}>0$, $\forall\,i$). In each iteration, given a current $b$, we select $i\in[n]$ uniformly at random and only the update $i$-th block of coordinates $b_{i}$. The next iterate $b^{+}$ is $b_{i}^{+}=T_{i}(b_{i})$ for $i$ and $b^{+}=b$ for the remaining blocks. Here, $T_{i}(b_{i})$ equals

[TABLE]

where $\alpha_{i}>0$ is the stepsize. It is well-known that (13) is equivalent to the following simple, explicit update formula:

[TABLE]

where $Z_{i}$ is a normalization constant such that $\sum_{j}b_{ij}^{+}=B_{i}$.

Algorithm 2 states the full block coordinate proportional response dynamics.

Convergence Analysis.

The objective function of (S) is relatively smooth with $L=1$ (Birnbaum, Devanur, and Xiao 2011, Lemma 7). This means that $\alpha_{i}=1$ is a safe lower bound for stepsizes. For Algorithm 2, a last-iterate sublinear convergence rate is given by Gao et al. (2020) for $0<\alpha_{i}<\frac{1+\theta_{i}}{L_{i}}$, where $\theta_{i}$ is the Bregman symmetry measure introduced by Bauschke, Bolte, and Teboulle (2017). For the KL divergence $\theta_{i}=0$. Their proof still goes through for $\alpha_{i}=1/L_{i}$, which yields the following

Theorem 2**.**

Let $b^{k}$ be random iterates generated by Algorithm 2 with $\alpha_{i}^{k}=1/L_{i}$ for all $k$, then

[TABLE]

where $\Phi^{*}$ is the optimal objective value.

Line search.

To speed up BCPR, we introduce a line search strategy and an adaptive stepsize strategy. BCPR with line search can be implemented by comparing $D_{\rm KL}(p^{+},p)$ and $D_{\rm KL}(b_{i}^{+},b_{i})$, which takes $O(m)$ time, and is much cheaper than computing the whole objective function value of (S). Beyond that, by storing $p$, we also avoid touching all variables in each iteration. Therefore, the amount of data accessed and computation needed is $O(m)$ per iteration (vs. $O(nm)$ for full gradient methods). In the adaptive strategy we compute larger stepsize based on Lipschitz estimates using a closed-form formula. The BCPR with Line Search (BCPR-LS) and Adaptive BCPR (A-BCPR) are formally stated in Appendix C.

As with BCDEG, our experiments demonstrate that larger stepsizes can accelerate Algorithm 2. When we consider a series of stepsizes $\{\alpha_{i}^{k}\}_{k\in\mathbb{N}}$ generated by line search or an adaptive strategy, we can show the inequality

[TABLE]

However, we cannot guarantee convergence, as we are unable to ensure the convergence of the last term above.

6 Proportional Response with Line Search

In this section we extend vanilla PR dynamics to PR dynamics with line search (PRLS), by developing a Mirror Descent with Line Search (MDLS) algorithm. Intuitively, the LS strategy is based on repeatedly incrementing the stepsize and checking the relative-smoothness condition, with decrements made when the condition fails. This is similar to the projected gradient method with line search (PGLS) in Gao and Kroer (2020, A.5), but replaces the $\ell_{2}$ norm with the Bregman divergence. This allows larger stepsizes while guaranteeing the same sublinear last-iterate convergence. The general MDLS algorithm is stated in Appendix D.

Convergence rate.

Birnbaum, Devanur, and Xiao (2011, Theorem 3) showed that a constant stepsize of $\alpha^{k}=1/L$ ensures sublinear convergence at a $1/k$ rate. One of the key steps in establishing the rate is the descent lemma, which also holds in the line search case:

Lemma 3**.**

Let $b^{k}$ be MDLS iterates. Then, $D(b^{*},b^{k+1})\leq D(b^{*},b^{k})$ for all $k$.

We have the following theorem.

Theorem 3**.**

Let $b^{k}$ be iterates generated by PRLS starting from any initial solution feasible $b^{0}$, then we have

[TABLE]

where $\rho^{-}$ is the shrinking factor of stepsize.

Unlike the result in BCPR-LS, the deterministic algorithm maintains its convergence guarantee with line search. The main issue in the block coordinate case is the lack of monotonicity of $\mathbf{E}[D_{i}(b^{*},b^{k})]$, which is avoided by the deterministic algorithm. In the proof, we also give a tighter, path-dependent bound.

7 Block Coordinate Descent Algorithm for CES utility function

In this section, we show that block coordinate descent algorithms also work for the case where buyers have constant elasticity of substitution (CES) utility functions (for $\rho\in(0,1)$).

Formally, a CES utility function for buyer $i$ with parameters $v_{i}\in\mathbb{R}^{m}_{+}$ and $\rho\in(0,1)$ is $u_{i}(x_{i})=(\sum_{j}{v_{ij}x_{ij}^{\rho}})^{\frac{1}{\rho}}$ where $v_{ij}$ denotes the valuation per unit of item $j$ for buyer $i$. CES utility functions are convex, continuous, nonnegative and homogeneous and hence the resulting ME can still be captured by the EG convex program.

BCDEG for CES utility function.

The resulting EG program is of similar form as (1) with $f(x)=-\sum\nolimits_{i}{\frac{B_{i}}{\rho}\log{\langle v_{i},x_{i}^{\rho}\rangle}}$, where $x_{i}^{\rho}=(x_{i1}^{\rho},\ldots,x_{im}^{\rho})$ and the same separable $\psi(x)$ as (1). Hence, as for linear Fisher markets, we can apply block-coordinate descent.

Since $u$ and $x$ may reach [math] at some iterates, we need smooth extrapolation techniques to ensure the existence of gradients. First, similar to Zhang (2011, Lemma 8), we lower bound $x^{*}$ for all $i,j:v_{ij}>0$, which ensures that extrapolation will not affect the equilibrium when $\min_{i,j}v_{ij}>0$. Our bounds are tighter than Zhang (2011).

Lemma 4**.**

For a market with CES utility functions with $\rho\in(0,1)$ and market equilibrium allocation $x^{*}$, for any $i,j$ such that $v_{ij}>0$, we have $x_{ij}^{*}\geq\underline{x}_{ij}^{*}:=\omega_{1}(i)^{\frac{1}{1-\rho}}\omega_{2}(i)^{\frac{\rho(\rho+1)}{1-\rho}}$, where $\omega_{1}(i)=\frac{B_{i}}{m\sum_{l}B_{l}}$, $\omega_{2}(i)=\min_{j:v_{ij}>0}v_{ij}/\max_{j}v_{ij}$.

In Appendix E, we show how to use this bound to perform safe extrapolation. This yields the following theorem for applying BCDEG to CES utilities. Due to its similarity to the linear case, we give the full algorithm in the appendix. The theorem is a direct consequence of Lemma 4 and Richtárik and Takáč (2014, Theorem 7).

Theorem 4**.**

Let $x^{k}$ be the random iterates generated by BCDEG  for CES utility function ($\rho\in(0,1)$) with stepsizes $\eta_{j}^{k}=1/L_{j}\;\forall\;k$, where

[TABLE]

*and $\underline{u}_{i}(\rho)=\sum_{j}v_{ij}{\underline{x}^{*}_{ij}}^{\rho}$ (defined in Lemma 4). Then, *

[TABLE]

where $\mu(L)$ is the strong-convexity modulus w.r.t. the weighted norm $\sum_{j}L_{j}\lVert\cdot\rVert^{2}_{\cdot j}$.

BCPR for CES utility function.

Unlike for linear utilities, the EG program for CES utility cannot be converted to a simple dual problem. Hence, we cannot view PR for $\rho\in(0,1)$ as a mirror-descent algorithm and analyze it with typical relative smoothness techniques. However, Zhang (2011) nonetheless showed convergence of PR for $\rho\in(0,1)$. We show that we can still extend their proof to show convergence of block coordinate PR for CES utility.

Theorem 5**.**

Let $b^{k}$ be the random iterates generated by BCPR  for CES utility function ($\rho\in(0,1)$). For any $\epsilon>0$, when

[TABLE]

we have $\mathbf{E}\left[\frac{\lvert b_{ij}-b_{ij}^{*}\rvert}{b_{ij}^{*}}\right]\leq\epsilon$ for all $i,j$ such that $v_{ij}>0$.

8 Numerical Experiments

We performed numerical experiments based on both simulated (for linear and CES ($\rho\in(0,1)$) utilities) and real data to test the scalability of our algorithms.

To measure the amount of work performed, we measure the number of accesses to cells in the valuation matrix. For the deterministic algorithms, each iteration costs $n\times m$, whereas for our BCDEG algorithms each iteration costs $n$, and for our BCPR algorithms each iteration costs $m$. For algorithms that employ line search, we count the valuation accesses required in order to perform the line search as well.

To measure the accuracy of a solution, we use the duality gap and average relative difference between $u$ and $u^{*}$. The instances are small enough that we can compute the equilibrium utilities $u^{*}$ using Mosek (2010).

Simulated low-rank instances.

To simulate market instances, we generate a set of valuations that mimic approximately low rank valuations, which are prevalent in real markets, and were previously studied in the market equilibrium context by Kroer et al. (2019). The valuation for item $j$ and buyer $i$ is generated as: $v_{ij}=v_{i}v_{j}+\epsilon_{ij}$, where $v_{i}\sim\mathcal{N}(1,1)$, $v_{j}\sim\mathcal{N}(1,1)$, and $\epsilon_{ij}\sim\emph{uniform}(0,1)$. Here, buyer $i$’s valuation for item $j$ consists of three parts: a value of item $j$ itself ($v_{j}$), buyers $i$’s average valuation ($v_{i}$), and a random term $\epsilon_{ij}$. We consider markets with $n=m=400$. All budgets and supplies are equal to one.

Movierating instances.

We generate a market instance using a movie rating dataset collected from twitter called Movietweetings (Dooms, De Pessemier, and Martens 2013). Here, users are viewed as the buyers, movies as items, and ratings as valuations. Each buyer is assigned a unit budget and each item unit supply. We use the “snapshots $200$K” data set and remove users and movies with too few entries. Using the matrix completion software fancyimpute (Rubinsteyn and Feldman 2016), we estimate missing valuations. The resulting instance has $n=691$ buyers and $m=632$ items.

First we compare each of our new algorithms in terms of the different stepsize strategies: BCDEG vs. BCDEG-LS, BCPR vs. A-BCPR vs. BCPR-LS, and PR vs PRLS. The results are shown in Fig. 1 and Fig. 2 in the upper left, upper right, and lower left corners. In all cases we see that our new line search variants perform the best.

Second, we then compare our block-coordinate algorithms to the best deterministic state-of-the-art market equilibrium algorithms: PGLS (Gao and Kroer 2020) and PRLS. The results are shown in Fig. 1 and Fig. 2 in the lower right corner.

For all markets either BCDEG or BCPR-LS is best on all metrics, followed by PRLS, and PGLS in order. In general, we see that the stochastic block-coordinate algorithms converge faster than their deterministic counterparts across the board, even after thousands of iterations of the deterministic methods. Thus, block-coordinate methods seem to be better even at high precision, while simultaneously achieving better early performance due to the many more updates performed per unit of work.

CES utilities.

Similar to linear utilities, we generate a $n=m=200$ scale market instance with the CES-utility parameters generated as $v_{ij}=v_{i}v_{j}+\epsilon_{ij}$, where $v_{i}\sim\mathcal{N}(1,0.2^{2})$, $v_{j}\sim\mathcal{N}(1,0.2^{2})$, and $\epsilon_{ij}\sim\emph{uniform}(0,0.2)$. For the CES instances we were not able to obtain high-accuracy solutions from existing conic solvers, and thus we only measure performance in terms of the dual gap. Somewhat surprisingly, we find that for CES utilities vanilla PR converges very fast; faster than BCPR and all versions of BCDEG. Moreover, BCDEG suffers from extremely small stepsizes, so we have to use BCDEG-LS. Here, we used specific small parameters for the distributions in the simulated utilities. More discussion and experiments on CES utilities are given in Appendix E.

9 Conclusion and Future Work

We proposed two stochastic block-coordinate algorithms for computing large-scale ME: (proximal) block-coordinate descent on EG (BCDEG) and block-coordinate proportional response (BCPR). For each algorithm we provided theoretical convergence guarantees and showed numerically that they outperform existing state-of-the-art algorithms. We also provided a new economic interpretation of the projected gradient update used in BCDEG. For future work, we are interested in deriving a sublinear convergence rate for BCPR with adaptive stepsizes, extending it to leverage distributed and parallel computing capabilities, and allowing more general dynamic settings and other buyer utility models.

Appendix A Block Coordinate Descent Algorithm for the EG Program

An efficient implementation of BCDEG

In order to implement Algorithm 1 efficiently, we temporarily store a vector $u\in\mathbb{R}^{n}$ in memory, which represents the current utilities for all buyers. Then, at each iteration $k$, we can dynamically update $u$ by only substituting $x_{\cdot j_{k}}$ while updating a particular block $j_{k}$, which can be done in $O(n)$ time. This avoids performing calculations of the form $u_{i}=\langle v_{i},x_{i}\rangle$, which require $O(nm)$ time when performed for each buyer. See Algorithm 3 for details.

(Proximal) Block Coordinate Descent for the EG Program with Line Search

We formally state BCDEG-LS here. This algorithm often outperforms others in our numerical experimental results. The difference to BCDEG is that we check a descent condition (the if condition in Algorithm 4) after performing a tentative proximal step. If the descent condition fails the check (the true case of the if-else statement), then we decrease the stepsize and perform another tentative step. If it passes the check then we commit to the tentative proximal step.

Proof of Eq. 5

Proof.

It can be verified by simple calculations that both $\tilde{f}_{i}(x)$ and $\nabla\tilde{f}_{i}(x)$ are continuous on their domains. Each entry of Hessian matrix of $\tilde{f}$ is of the form:

[TABLE]

Let $x$ and $y$ be two vectors which only differ in one coordinate $j$. Then,

[TABLE]

where the inequality is due to

[TABLE]

Therefore,

[TABLE]

that is,

[TABLE]

To derive (5), we can apply integration:

[TABLE]

∎

Remark.

As a comparison, we provide a proof for a standard (global) Lipschitz constant $L$ as follows. Comparing this with the above proof, we show that when only one block of the variables $x_{\cdot j}$ change, the coordinate-wise Lipschitz constant $L_{j}$ is bounded by the maximum (over all $i$) diagonal value of the Hessian matrix. In contrast, for $L$ it depends on the maximum (over all $i$) sub-trace of the Hessian matrix.

Lemma 5**.**

For any $j\in[m]$, let

[TABLE]

Then, for all $x,y\in\mathcal{X}$, we have

[TABLE]

Proof.

Let $x$ and $y$ be two vectors in the domain of $\tilde{f}$, then

[TABLE]

where $y^{j^{\prime}}_{i}$ is a vector whose first $(j^{\prime}+1)$’th components are copy of the first $(j^{\prime}+1)$’th components of $y_{i}$ and other components are copy of the last $(m-j^{\prime}-1)$’th components of $x_{i}$; $x^{j^{\prime}}_{i}$ is a vector whose first $j^{\prime}$’th components are copy of the first $j^{\prime}$’th components of $y_{i}$ and other components are copy of the last $(m-j^{\prime})$’th components of $x_{i}$. That is, each pair of $y^{j^{\prime}}_{i}$ and $x^{j^{\prime}}_{i}$ only differ in $(j^{\prime}+1)$’th component. Also, we have $|y_{i}^{j^{\prime}}-x_{i}^{j^{\prime}}|\leq\|y_{i}-x_{i}\|$ for all $i,j^{\prime}$.

[TABLE]

since

[TABLE]

Therefore,

[TABLE]

∎

Note that this is a tighter bound than Gao and Kroer (2020) because $\max_{i\in[n]}\frac{B_{i}\lVert v_{i}\rVert^{2}}{\underline{u}_{i}^{2}}\leq\max_{i\in[n]}\frac{B_{i}}{\underline{u}_{i}^{2}}\max_{i\in[n]}\lVert v_{i}\rVert^{2}$.

Proof of Lemma 2

Proof.

Since $g(u)$ is strongly convex, there exists a unique $u^{*}$ such that

[TABLE]

Thus, the set of optimal solutions $\mathcal{X}^{*}$ can be described by the following polyhedral set for some $C$:

[TABLE]

Assume that the optimal polyhedral set $\mathcal{X}^{*}$ is non-empty, then the Hoffman inequality (Hoffman 2003) tells us there exists some positive constant such that

[TABLE]

Let $\mathcal{X}$ be the set of feasible allocations. For any $x\in\mathcal{X}$ such that $Cx\leq d$, we have

[TABLE]

From strong convexity of $\tilde{g}$ and (29), we have

[TABLE]

If $L\geq\frac{\mu}{\theta^{2}(A,C)}$, then

[TABLE]

Hence,

[TABLE]

Since $\mathcal{D}(x,\cdot)$ is non-decreasing (Karimi, Nutini, and Schmidt 2016, Lemma 1), we have

[TABLE]

If $L<\frac{\mu}{\theta^{2}(A,C)}$, then

[TABLE]

Hence,

[TABLE]

Combining the two cases, we can conclude (6). ∎

Proof of Theorem 1

The following convergence analysis is typical for uniform (block) coordinate descent algorithms. Leveraging (6), it is easy to derive a linear convergence rate for both BCDEG and BCDEG-LS. This proof scheme was given in Karimi, Nutini, and Schmidt (2016, Appendix H), but here we give a variant with different stepsizes across blocks and iterations.

Proof.

To remind readers of notations, we define (again)

[TABLE]

for BCDEG and

[TABLE]

for BCDEG-LS.

Then, we have (e.g. for BCDEG)

[TABLE]

where we used Lipschitz smoothness (or line search condition for BCDEG-LS) of $\tilde{f}$, the definition of $T_{j}(x)$, and separability on each block $j$.

For BCDEG-LS, replacing (38) with the following inequality:

[TABLE]

and by $1/\eta_{j}^{k}\leq L_{j}\leq L_{\max}$ (used in the last inequality) we have (39).

By Lemma 2, we have

[TABLE]

Subtracting $F^{*}$ on the both sides, Rearranging the above inequality, and taking expectation w.r.t. $j_{1},\ldots,j_{k}$, we obtain (7) by induction. ∎

Appendix B Algorithm Details for Projection onto a Simplex

We show a key theorem for projection onto a simplex algorithm and a formally stated algorithm as background for readers who are not familiar with this algorithm.

Theorem 6**.**

(Chen and Ye 2011, Theorem 2.2.)* For any vector $y\in\mathbb{R}^{n}$, the projection of $y$ onto $\Delta^{n}$ is obtained by the positive part of $y-\hat{t}$:*

[TABLE]

where $\hat{t}$ is the only one in $\{t_{i}:i=0,\ldots,n-1\}$ that falls into the corresponding interval as follows,

[TABLE]

Hence, to find $x$, we only need to find the $t_{i}$ in (43) that falls into the corresponding interval, claim it as the optimal $\hat{t}$, and then compute $x$ based on (42). The following algorithm is implied by the above theorem.

Discussion on alternative algorithms for projection onto a simplex

In this paper, we assume that the above algorithm is an efficient method to project an array onto a simplex, and interpret projection-type ME-solving algorithms based on this algorithm. There are several variants of this type of projection algorithm, and some of them give improved complexity results or practical performance. Condat (2016) summarize most of these algorithms, which use different methods to find the correct threshold value. Some of these algorithms exploit problem structure to obtain $O(n)$ expected or practical complexity. However, in the worst case, there is still no better result than $O(n\log(n))$. Perez et al. (2020) proposed a method with a worst-case linear-time complexity result. Their result is analogous to bucket sort, it assumes that the set of possible values that you might encounter is constant.

Appendix C Block Coordinate Proportional Response

Block Coordinate Proportional Response with Line Search

We formally state Block Coordinate Proportional Response with Line Search (BCPR-LS) algorithm as follows. Note that $\delta\in[0,1)$ is a “conservative” factor - we can set $\delta>0$ to use more conservative stepsize strategy than standard line search ($\delta=0$). In $\delta>0$ case, we can guarantee some (weak) convergence property for BCPR-LS.

Proof of Eq. 14

Proof.

The original minimization problem can be solved by Lagrangian method by introducing a Lagrangian multiplier $\lambda$:

[TABLE]

Let $\nabla_{a_{1},a_{2},\ldots,a_{m},\lambda}\mathcal{L}=0$, we have

[TABLE]

The partial derivatives are equal to

[TABLE]

Let a_{j}=\frac{1}{Z}b_{ij}\Big{(}\frac{v_{ij}}{p_{j}}\Big{)}^{\alpha}, $\lambda=1+\frac{1}{\alpha}-\log{Z}$, where Z=\frac{1}{B_{i}}\sum_{j}{b_{ij}\Big{(}\frac{v_{ij}}{p_{j}}\Big{)}^{\alpha}}. We can verify that

[TABLE]

Therefore, $(a_{1},a_{2},\ldots,a_{m},\lambda)$ is a stationary point of the Lagrange function. $a=(a_{1},a_{2},\ldots,a_{m})$ is a solution to the original minimization problem due to the convexity of the problem. ∎

Coordinate-wise relative smoothness (in expectation)

In the rest of this section, except for A-BCPR, we use general Bregman distance $D(a,b)$ between $a,b\in\mathbb{R}^{n}$ and $D_{L}(a,b)=\sum_{i}L_{i}D_{i}(a_{i},b_{i})$ where $D_{i}(a_{i},b_{i})$ is coordinate-wise Bregman distance and $L\in\mathbb{R}^{n}$. We also introduce $T^{k}=\sum_{i}T_{i}^{k}$ as below (similar to $l_{2}$ norm block coordinate descent):

[TABLE]

Then, it is standard to have the following lemma.

Lemma 6**.**

If $b^{k}$ is the random iterates generated by Algorithm 6, then

[TABLE]

Note that (48) is a special case (the cardinality of sampling equals to $1$, uniform sampling) of $h$-ESO assumption defined in Hanzely and Richtárik (2021). Essentially, this is coordinate-wise relative smoothness in expectation.

Proof.

Given any coordinate $i$, we have (let $b^{+}$ and $p^{+}$ be the next iterates of $b$ and $p$, respectively)

[TABLE]

and $D_{\rm KL}(b_{i}^{+},b_{i})=\sum_{j}{b_{ij}^{+}\log{\frac{b_{ij}^{+}}{b_{ij}}}}$. The line search condition ensures

[TABLE]

Note that if we let $a\in\mathbb{R}$ and $f(x)=x\log{\frac{x}{x-a}}$, then $f^{\prime}(x)=\log{(1+\frac{a}{x-a})}-\frac{a}{x-a}\leq\frac{a}{x-a}-\frac{a}{x-a}=0$. Since $p_{j}^{+}-p_{j}=b_{ij}^{+}-b_{ij}$, let $a=b_{ij}^{+}-b_{ij}$, we have

[TABLE]

where we use $p_{j}^{+}\geq b_{ij}^{+}$ and monotonity of $f(x)$. Note that (51) always holds when $\alpha_{i}=1$ for all $i$.

Hence, by separability on each block of coordinates we have

[TABLE]

As we know

[TABLE]

and

[TABLE]

(48) and (49) are equivalent.

∎

A useful descent lemma

Next, we will consider BCPR with line search (BCPR-LS). We can directly see that BCPR (with fixed stepsizes)’s sublinear convergence rate is a special case of our convergence property.

The following lemma can be derived from Lemma 6 and three point property (e.g. Hanzely and Richtárik (2021, Lemma 3.5)). Beyond existing results, we prove a version for line search and $\alpha_{i}^{k}=(1-\delta)/L_{i}^{k}$ stepsize strategy.

Lemma 7**.**

If $b^{k}$ is the random iterates generated by Algorithm 6, then

[TABLE]

where $T^{k}$ defined in (47) and $u\in{\rm dom}\,\varphi\cup{\rm dom}\,h$.

Proof.

Given any $b^{k}$ (iterate at the $k$th iteration), we have

[TABLE]

which follows from Lemma 6 and separability of $r$.

By the optimality of (47), we have

[TABLE]

for some $v^{k}\in\partial r(T^{k})$.

By convexity of $\varphi$ and $r$, for any $u\in{\rm dom}\,\varphi\cup{\rm dom}\,h$,

[TABLE]

[TABLE]

By the definition of Bregman distance we have so-called three point identity:

[TABLE]

for all $u_{i}\in{\rm dom}\,h_{i}$ and all $i\in[n]$.

Then, by combining all above inequalities and equalities we have

[TABLE]

for any $u$. Here, we use (57), (59) with (58), (60), (61) and $\alpha_{i}^{k}=(1-\delta)/L_{i}^{k}$ for $\delta\in[0,1)$. (56) follows from identity relations similar to (54) and (55). ∎

The following corollaries provide some (weak) convergence properties for BCPR-LS.

Corollary 1**.**

If $b^{k}$ is the random iterates generated by Algorithm 6, then

1. (i)

$\mathbf{E}[\Phi(b^{k})]$* is non-increasing;* 2. (ii)

$\sum_{k=0}^{\infty}\mathbf{E}[D_{1/\alpha^{k}}(b^{k+1},b^{k})]<\infty$, $\mathbf{E}[D_{i}(b_{i}^{k+1},b_{i}^{k})]\rightarrow 0$ as $k\rightarrow\infty$ for all $i$; 3. (iii)

$\sum_{k=0}^{\infty}\mathbf{E}[D_{1/\alpha^{k}}(b^{k},b^{k+1})]<\infty$, $\mathbf{E}[D_{i}(b_{i}^{k},b_{i}^{k+1})]\rightarrow 0$ as $k\rightarrow\infty$ for all $i$; 4. (iv)

$\mathbf{E}[D_{1/\alpha^{k}}(b^{*},b^{k})]\leq\mathbf{E}[D_{1/\alpha^{k}}(b^{*},b^{k+1})]$* as $k\rightarrow\infty$.*

Proof.

Let $u=b^{k}$ in (56), we have

[TABLE]

by using non-negativity of Bregman distance. This shows objective descent in expectation.

Also, we have

[TABLE]

As we know $\mathbf{E}[D_{\alpha^{k}}(b^{k+1},b^{k})|b^{k}]=\frac{1}{n}D_{\alpha^{k}}(T^{k},b^{k})$, (64) is equivalent to

[TABLE]

Summing up (65) over $0,1,\ldots,k-1$ and taking expectation w.r.t. $i_{0},i_{1},\ldots,i_{k-1}$, we have

[TABLE]

Let $k\rightarrow\infty$, we have $\sum_{l=0}^{\infty}\mathbf{E}[D_{\alpha^{l}}(b^{l+1},b^{l})]<\infty$, which implies $\mathbf{E}[D_{\alpha^{k}}(b^{k+1},b^{k})]\rightarrow 0$ as $k\rightarrow\infty$.

Similar to (64), we have

[TABLE]

Then, (iii) follows from the same argument as the proof of (ii).

Because $\mathbf{E}[\Phi(b^{k})]$ is non-increasing and $\mathbf{E}[\Phi(b^{k})]\geq\Phi^{*}>-\infty$, $\mathbf{E}[\Phi(b^{k+1})-\Phi(b^{k})]\rightarrow 0$.

Let $u=b^{*}$ in (62), we have

[TABLE]

due to the non-negativity of Bregman distance. Rearranging (68), summing it up over $0,\ldots,k$ and using the tower property, we have

[TABLE]

When $k\rightarrow\infty$, the left hand side of (69) goes to be non-negative, which imples (iv). ∎

Proof of Theorem 2

Proof.

For any $u\in{\rm dom}\,h$, we have

[TABLE]

By Lemma 7 and (70), let $u=b^{*}$, then we get

[TABLE]

Rearranging (71), we attain

[TABLE]

We take the expectation of (72) with respect to $\{i_{0},i_{1},\ldots,i_{k}\}$ and sum this inequality over $l=0,\ldots,k$ to obtain

[TABLE]

It is easy to see that when $\alpha_{i}^{k}=1/L_{i}$, the last term in (73) is zero and sublinear convergence rate follows after re-arrangement. For line search case, we can derive (16).

Remark.

Even though there are many equivalent formulas of the last term of (73) (see the following for examples), we cannot guarantee its convergence.

[TABLE]

∎

Adaptive Block Coordinate descent algorithm for PR dynamics (A-BCPR)

We also give an adaptive stepsize strategy. To do this, we find a series of $\{L_{i}^{k}\}_{k\in\mathbb{N}}$ with each we have coordinate-wise relative smoothness with $L_{i}^{k}<1$ at the $k$th iteration. We state the following lemma to support this algorithm.

Lemma 8**.**

Let $b_{ij}^{k+1}$ be the next iterates of $b_{ij}^{k}$ in BCPR with $1\leq\alpha_{i}^{k}\leq\bar{\alpha}_{i}^{k},\;\forall\;i,k$, where $1$ and $\bar{\alpha}_{i}^{k}$ are lower and upper bounds of stepsizes, respectively. Let

[TABLE]

and

[TABLE]

When $1\leq\beta_{i}^{k}\leq\sqrt{2}$, we have

[TABLE]

where

[TABLE]

Proof.

Let $c_{j}=p_{j}^{+}-b_{ij}^{+}=p_{j}-b_{ij}$. The Taylor series of $D_{\rm KL}(c_{j}+b_{ij}^{+},c_{j}+b_{ij})$ with the Lagrange form of the remainder is

[TABLE]

for some $\lambda\in[0,1]$.

We assume the unit budgets, i.e., $\sum_{j}b_{ij}=1,\forall i$. Since $1\leq\alpha_{i}\leq\bar{\alpha}_{i}$, we can derive

[TABLE]

Similarly, we have

[TABLE]

Then, we obtain

[TABLE]

On the other side, we have when $1<\beta_{i}\leq\sqrt{2}$,

[TABLE]

To verify (83), we can re-arrange it and check this inequality is equivalent to

[TABLE]

The above inequality can be verified by considering two cases. If $b_{ij}^{+}\geq b_{ij}$,

[TABLE]

where we use $b_{ij}^{+}-b_{ij}\leq 2\left(c_{j}+b_{ij}+\lambda(b_{ij}^{+}-b_{ij})\right)$ and $b_{ij}^{+}-b_{ij}\leq b_{ij}^{+}-\frac{1}{2\beta_{i}}b_{ij}$ when $1<\beta_{i}\leq 3$.

If $b_{ij}^{+}<b_{ij}$,

[TABLE]

where we use $b_{ij}-b_{ij}^{+}\leq b_{ij}^{+}$ holds when $1<\beta_{i}\leq 2$, and $b_{ij}-b_{ij}^{+}\leq\frac{2}{\beta_{i}^{2}}\left(b_{ij}^{+}-\frac{1}{2\beta_{i}}b_{ij}\right)$ holds when $1<\beta_{i}\leq\sqrt{2}$.

Combining (79), (82) and (83), we attain

[TABLE]

where the second inequality follows since all terms are nonnegative in the two summations. Lemma 8 follows because the right-hand side of (87) is increasing in $\frac{b_{ij}}{p_{j}}$. ∎

Note that by the buyer optimality condition, when $b_{ij}^{*}>0$, $v_{ij}/p_{j}^{*}=\max_{j}v_{ij}/p_{j}^{*}$. This implies $\beta_{i}^{*}=1$. Hence, Lemma 8 is meaningful.

We state adaptive block coordinate descent algorithm for PR dynamics as follows.

Appendix D Proportional Response with Line Search

We state general Mirror Descent with Line Search (MD-LS) here.

Proof of Theorem 3

Proof.

We prove sublinear convergence for general Mirror Descent with Line Search, which immediately implies Theorem 3 as we know $L=1$ for (S).

Same as Birnbaum, Devanur, and Xiao (2011, Lemma 5), it can be easily seen that Algorithm 9 is still a descent objective algorithm, that is,

[TABLE]

By the convexity of $f$ and the break condition in Algorithm 8, we know that

[TABLE]

Using Birnbaum, Devanur, and Xiao (2011, Lemma 9) (originally in Chen and Teboulle (1993)), we have, for any $u\in\mathcal{X}$,

[TABLE]

and then

[TABLE]

The above is similar to Birnbaum, Devanur, and Xiao (2011, Corollary 10). Let $b^{*}$ be an optimal solution. Then, similar to Birnbaum, Devanur, and Xiao (2011, Lemma 8), we have

[TABLE]

This implies Lemma 3 because

[TABLE]

Summing up (92) across $k=0,1,\ldots$, we have

[TABLE]

where the first inequality follows from (88). Rearranging the above gives the following path-dependent bound:

[TABLE]

By Lemma 3 and $\alpha^{l}\geq\frac{\rho_{-}}{L}$, (17) follows from

[TABLE]

∎

Appendix E Block Coordinate Descent Algorithm for CES Utility Function

Details of quadratic extrapolation for CES utility function

In this section, we use $u_{i}^{\rho}$ to substitute $u_{i}(\rho)$ to avoid confusion as the decision variable of $u_{i}(\rho)$ is $x_{i}$.

Denote the objective function of EG convex program for CES utility function as

[TABLE]

where

[TABLE]

Then, we replace $g$ and $u^{\rho}$ with $\tilde{g}$ and $\tilde{u}^{\rho}$ (defined as follows), respectively. Similar to the linear case,

[TABLE]

and

[TABLE]

where $\underline{u}_{i}$ and $\underline{x}_{ij}$ are specific lower bounds for $u^{*}_{i}$ and $x^{*}_{ij}$ for some $i,j$. By similar discussion as Gao and Kroer (2020, Lemma 1), we know $\underline{u}_{i}$ exists for all $i$. By Lemma 4, we know $\underline{x}_{ij}$ exists for all $i,j$.

Proof of Lemma 4

Here, we provide a proof for Lemma 4, which leverages fixed point property at equilibrium like Zhang (2011, Lemma 8), but generates a tighter bound.

Proof.

Consider a buyer $i$ and an item $j$ with $v_{ij}>0$. If for the buyer $i$, all other $k\neq j$ is associated with $v_{ik}=0$, obviously $b_{ij}^{*}=B_{i}$ and $x_{ij}^{*}=b_{ij}^{*}/p_{j}^{*}\geq B_{i}/\sum_{l}B_{l}$.

Otherwise, if there is $v_{ik}>0$ for some $k\neq j$, since a market equilibrium is a fixed point of the proportional response dynamics (and this can be extended to block coordinate case naturally), we have

[TABLE]

Rearrange (101) and then we get

[TABLE]

and

[TABLE]

To be specific, we consider $b_{ij}^{*}<\frac{B_{i}}{m}$ because otherwise the lower bound is directly given. For these $i,j$, we can always choose another $k$ (there is at least one as discussed above) such that $b_{ik}^{*}\geq\frac{B_{i}}{m}$. Hence, we can bound $p_{j}^{*}$ by $b_{ij}^{*}\leq p_{j}^{*}$, (101) and $b_{ik}^{*}\leq p_{k}^{*}$:

[TABLE]

where $\omega_{2}(i)=\min_{j:v_{ij}>0}v_{ij}/\max_{j}v_{ij}$. Combining (103) with (104), we obtain

[TABLE]

where $\omega_{1}(i)=\frac{B_{i}}{m\sum_{l}B_{l}}$. All other cases give greater lower bounds than this. ∎

Proof of Theorem 4

We use the following lemma to establish a smoothness guarantee for $\tilde{f}(x)=\sum_{i}\tilde{g}(\sum_{j}\tilde{u}_{ij}^{\rho}(x_{ij}))$.

Lemma 9**.**

For any $j\in[m]$, let

[TABLE]

Then, for all feasible $x,y$ such that $x,y$ differ only in the $j$th block, we have

[TABLE]

Proof.

Let $\tilde{u}_{i}$ stand for $\tilde{u}_{i}(\rho)$ for conciseness. We can obtain

[TABLE]

Then, (107) follows from the same derivation as Eq. 5. ∎

Theorem 4 can be obtained by Richtárik and Takáč (2014, Theorem 7) and the fact that there exists a strong-convexity modulus $\mu(L)$ for $\tilde{f}$ w.r.t. the norm $\sum_{j}L_{j}{\lVert\cdot\rVert}_{\cdot j}^{2}$ when $\tilde{f}$ is of the form described above.

Remark.

Similar to Eq. 107, we can establish a standard Lipschitz smoothness with constant $L=\max_{i\in[n]}\frac{\rho B_{i}\lVert v_{i}\rVert^{2}}{\underline{u}_{i}^{2}}\underline{x}_{i}^{2\rho-2}$. This can be seen a generalized result from Lemma 5 as we have

[TABLE]

Proof of Theorem 5

We make use of an auxiliary function introduced in Zhang (2011), which they use to show convergence for deterministic PR.

Proof.

In this proof, $D_{i}(\cdot,\cdot)$ denotes KL divergence of the $i$th block of the decision variables. Define

[TABLE]

thus $T_{ij}^{k}=B_{i}u_{ij}^{k}(\rho)/u_{i}^{k}(\rho)$.

Let $m(k,i)=D_{i}(b_{i}^{*},b_{i}^{k})$ and $m(k)=\sum_{i}m(k,i)$, then we have

[TABLE]

Similar to Zhang (2011), we obtain

[TABLE]

Hence,

[TABLE]

Since $D_{\rm KL}(p^{*},p^{k})\geq 0$ and $\sum_{i}\frac{B_{i}}{\rho}\log{\frac{\langle v_{i},{x_{i}^{*}}^{\rho}\rangle}{\langle v_{i},{x_{i}^{k}}^{\rho}\rangle}}\geq 0$ (by the optimality of EG program), we have

[TABLE]

After taking expectation w.r.t. $i_{1},i_{2},\ldots,i_{k}$, we attain $\mathbf{E}[m(k)]\leq(\rho^{\prime})^{k}m(0)$.

Note that we have the following lemma.

Lemma 10**.**

(Zhang 2011, Lemma 9)* For two positive sequences $\{b^{*}_{ij}\}_{i\in[n],j\in[m]}$ and $\{b_{ij}\}_{i\in[n],j\in[m]}$ such that $\sum_{i,j}b^{*}_{ij}=\sum_{i,j}b_{ij}$, let $\eta=\max_{i,j}\frac{\lvert b_{ij}-b^{*}_{ij}\rvert}{b^{*}_{ij}}$. Then,*

[TABLE]

We consider the case where $\eta<1$. From Zhang (2011, Lemma 8), we know $b^{*}_{ij}\geq\left(\frac{1}{W^{2}}\right)^{\frac{1}{1-\rho}}$ for any $i,j$ such that $v_{ij}>0$ where $W=\frac{n}{\min_{v_{ij}>0}v_{ij}\cdot\min_{i}B_{i}}$. Then, (114) tells us that

[TABLE]

We take expectation on both sides of the above inequality, and use Jensen’s inequality, then we have

[TABLE]

Hence, to guarantee $\mathbf{E}\frac{\lvert b_{ij}-b_{ij}^{*}\rvert}{b_{ij}^{*}}\leq\epsilon<1$, we need $8W^{\frac{2}{1-\rho}}\mathbf{E}\left[m(k)\right]\leq\epsilon^{2}$. Since $\mathbf{E}[m(k)]\leq(\rho^{\prime})^{k}m(0)$, it suffices to have

[TABLE]

which is equivalent to k\geq{2\log{\frac{\sqrt{8m(0)}W^{\frac{1}{1-\rho}}}{\epsilon}}}\Big{/}{\log{\frac{n}{n-1+\rho}}}.

∎

Discussion about BCD methods for CES utilities and more experimental results

We showed that for linear utilities, the EG program can be solved efficiently by applying block coordinate descent algorithms. In contrast to linear utilities, there are some difficulties in computing market equilibria with CES ($\rho\in(0,1)$) utilities by BCD-type methods.

The key issue is that, while we can obtain a lower bound for $x_{ij}^{*}$ (for those $i,j:v_{ij}>0$), this bound is too small to generate a reasonable (coordinate-wise) Lipschitz constant, since the lower bound appears in the the bound for $L_{j}$ (see (106)).

Even though our algorithms are equipped with line search strategies, we find that they still struggle to find efficient stepsizes. We even find it hard to select an appropriate initial stepsize. Often, the stepsizes are too small to reach equilibrium, especially when the number of iterations is not very large. This makes BCDEG-LS and PGLS perform worse than BCPR and PR, because the latter do not need to choose stepsizes; they have (sublinear) convergence under stepsize $1$.

This difficulty can be overcome when volatility of valuations is small - in this case, the lower bound of $x$ in Eq. 106 is larger.

See Fig. 4 for more experimental results on CES utilities. We used the model $v_{ij}=v_{i}v_{j}+\epsilon$ where $v_{i}\sim\mathcal{N}(1,\hat{v})$, $v_{j}\sim\mathcal{N}(1,\hat{v})$ and $\epsilon\sim\emph{uniform}(0,\hat{v})$, where $\hat{v}$ is a parameter to control volatility. Here, we set $n=m=200$.

Appendix F Discussion on wall-clock running times

In this paper, we compare different algorithms based on the number of times they query for a buyer-item valuation (i.e. the number of accesses to cells of the valuation matrix). In this sense, if we count one query as one unit cost, then item-block algorithms (BCDEG, BCDEG-LS) cost $n$ units in each iteration, buyer-block algorithms (BCPR, A-BCPR, BCPR-LS) cost $m$ units in each iteration, and full gradient algorithms (PGLS, PR) cost $nm$ units in each iteration. Another way to compare these algorithms would be to compare wall-clock running times. In this section we justify our choice of measuring valuation accesses instead.

For the projected-gradient style methods (the full deterministic method and BCDEG variants), they each perform one $n$-dimensional projection for every $n$ units of work, and thus all these methods are directly comparable in terms of our measure. Similarly, all proportional-response-style algorithms perform the same amount of work for every $m$ queries made. Thus they are also directly comparable in terms of our measure. When comparing projected-gradient-style methods to PR-style methods, there is a difference in that projection is required in the former. Yet in practice it is known that efficient implementations of the best projection-style algorithms tend to run in linear time as well  (Condat 2016). Thus cross comparisons should also be fair. The main reason why wall-clock might differ from our measure is due to efficiency of the implementation. All our experiments are implemented in python. For that reason, deterministic methods could have an advantage over block-coordinate methods in terms of wall clock, since NumPy and MKL enable highly-optimized matrix/vector operations. In turn, this means that deterministic methods can spend more time ‘in C’ as opposed to ‘in python.’ If all algorithms are implemented without highly optimized matrix operations that exploit multi-core CPUs, experimental results on wall-clock running times will be consistent with our query-based complexity cost. Similarly, if they are all optimized fully to take advantage of hardware, parallelization, etc., we can still expect similar results.