\newpage

\section{Preliminaries}
We list and prove a few elementary lemmas used in subsequent proofs and discussions.

\begin{lemma}
	Let $n \geq m$, $A\in \RR^{m \times n}$ and $b\in \RR^m $ be such that $\mathcal{S} = \{ x\in \RR^n: Ax=b \} \neq \emptyset$. Let the singular value decomposition (SVD) of $A^\top$ be $A^\top = U\Sigma V^\top$, where $U\in \RR^{n\times r}$ and $V\in \RR^{m \times r}$ have orthonormal columns, $\Sigma = \diag(\sigma_1, \dots, \sigma_r)$, $\sigma_1\geq \dots \sigma_r >0$, $r = \rank(A)$. For any $x\in \RR^n$, the projection of $x$ onto $\mathcal{S}$ can be expressed as
	\[ \Pi_\cS(x) = (I - UU^\top) x + U \Sigma^{-1} V^\top b. \]
	In particular, when $A$ has full rank, 
	\[ \Pi_\cS(x) = (I - A^\top (AA^\top)^{-1}A)x + A^\top (AA^\top)^{-1}b. \]
	\label{lemma:proj-subspace-expression}
\end{lemma}
\textit{Proof.} The case of full rank $A$ is well-known, see, e.g., \cite[Eq. (1)]{plesnik2007finding}. For general $A$, $Ax = b \Leftrightarrow V\Sigma U^\top x = b$. Since $V\Sigma$ has full column rank, there exists a unique $w\in \RR^{r}$ such that $V\Sigma w = b$. In fact, $w = \Sigma^{-1}V^\top b$. Therefore, \[Ax = b \Leftrightarrow U^\top x = \Sigma^{-1} V^\top b.\] 
Since $U^\top$ has full row rank, the formula follows directly from the full rank case.
\qed

\begin{lemma}
	Under the same assumptions as Lemma \ref{lemma:proj-subspace-expression}, for any $x\in \RR^n$, it holds that \[ \sigma_r \cdot \| x - \Pi_\mathcal{S}(x) \| \leq \|Ax-b\|. \] \label{lemma:hoffman-const-subspace}
\end{lemma}
\textit{Proof.} 
By Lemma \eqref{lemma:proj-subspace-expression} and the fact that $U$ and $V$ have orthonormal columns, 
\begin{align*}
	\|x - \Pi_\cS(x) \| &= \| U(U^\top x - \Sigma^{-1}V^\top  b) \| = \| U^\top x - \Sigma^{-1}V^\top  b \| \\
	& = \|  \Sigma^{-1}V^\top  (Ax -b) \| \leq \|\Sigma^{-1}\| \|Ax-b\|,
\end{align*}
where $\|\Sigma^{-1}\| = \frac{1}{\sigma_r}$.\qed


\begin{lemma}
	For $B>0$ and $g\in \RR^d$, let
	\[ x^* = \argmin \left\{ \langle g, x\rangle + \sum_{i=1}^d x_i \log x_i : x\geq 0,\, \ones^\top x = B \right\}. \] \label{lemma:md-subproblem}
	Then, $x^*_i = B\cdot \frac{e^{-g_i}}{\sum_\ell e^{-g_\ell}}$, $i\in [d]$. 
\end{lemma}
\textit{Proof.} It can be easily verified via KKT optimality conditions. The Lagrangian is 
\[ L(x, \lambda) = \langle g, x \rangle + \sum_i x_i \log x_i - \lambda(\ones^\top x - B). \]
By the first-order condition, for any $\lambda$, the minimizer $x(\lambda)$ of $L(x, \lambda)$ has $x_i(\lambda) = e^{\lambda - 1 - g_i}$. Primal feasibility implies $\sum_i x_i(\lambda) = B \Rightarrow e^\lambda = \frac{B}{\sum_i e^{-(1+g_i)}}$. Therefore, 
\[ x^*_i = B\cdot \frac{e^{-g_i}}{\sum_\ell e^{-g_\ell}}. \]
\qed 

\subsection{Proof of Theorem \ref{thm:eg-gives-me-for-certain-ui}} \label{app:thm-eg-and-mkt-equi}
Existence of an optimal solution of \eqref{eq:eisenberg-gale-primal} can be shown via a standard argument. Specifically, since $u_i$ are CCNH, the objective function $ F(x) = \sum_i B_i \log u_i(x_i)$ is continuous and concave on the interior of the feasible region $\mathcal{X} = \left\{ x\in \RR_+^{n\times m}: \sum_i x_{ij} \leq 1,\, \forall\, j \right\}$. Since $u_i(\ones) >0$ for all $i$, by homogeneity, $u_i(\ones/n) >0$ for all $i$ for the feasible allocation $x^0_i = \ones/n$, $i\in [n]$. Hence, the optimal objective value $>-\infty$. Meanwhile, for any feasible $x$, 
\[ F(x) \leq \sum_i B_i \log u_i(\ones) < \infty.  \]
Hence, $F^* = \sup_{x\in \mathcal{X}} F(x)$ is finite. Furthermore, 
\[ \mathcal{X}' = \left\{ x\in \mathcal{X},\, F(x) \geq F(x^0) \right\} \]
is a compact set, on which $F$ is continuous. Therefore, $F^* = F(x^*)$ for some $X^* \in \mathcal{X}$. 

Next, we show that an optimal solution of \eqref{eq:eisenberg-gale-primal} gives a ME. 
Consider the minimization form of \eqref{eq:eisenberg-gale-primal}. Let $p_j\geq 0$ be the dual variable associated with constraint $\sum_i x_{ij} \leq 1$. The Lagrangian is
\[\mathcal{L}(x,p) = \left[ -\sum_i B_i \log u_i(x) + \left \langle p, \sum_i x_i\right \rangle - \sum_j p_j \right]. \]
The Lagrangian dual is 
\begin{align}
\max_{p\geq 0} \, g(p) := \min_{x\geq 0}\, \mathcal{L}(x,p). \label{eq:general-eg-dual}
\end{align}

Since $u_i$ are homogeneous, \eqref{eq:eisenberg-gale-primal} has a strictly feasible solution with finite objective value ($\sum_i x_{ij} < 1$ for all $j$), e.g., $x_i = \ones/(n+1)$ for all $i$. 
The dual \eqref{eq:general-eg-dual} clearly has a strictly feasible solution $p>0$ with finite $g(p) = \min_{x\geq 0} \mathcal{L}(x, p)$. Therefore, strong duality holds by Slater's condition and the KKT conditions are necessary and sufficient for (primal and dual) optimality of a solution pair $(x,p)$ (see, e.g., \cite[Appendix D]{ben2019lectures}). Let $x^*$ and $p^*$ be optimal solutions to the primal \eqref{eq:eisenberg-gale-primal} and dual \eqref{eq:general-eg-dual}, respectively. Clearly, $u_i(x^*_i)>0$ for all $i$ (since $F^*$ is finite). 
By Lagrange duality, we have
\[x^* \in \argmin_{x\geq 0} \mathcal{L}(x, p^*).\]
In other words, each $x^*_i$ maximizes \[ r_i(x_i, p^*) := B_i \log u_i(x_i) - \langle p^*, x_i\rangle\] 
on $x_i \geq 0$. We show that $(x^*, p^*)$ is a market equilibrium. 

\paragraph{Buyer optimality} First, we verify that $\langle p^*, x^*_i\rangle = B_i$ for all $i$. Assume that $\langle p^*, x^*_i\rangle > B_i$ for some $i$. Let $\tilde{x}_i = (1-\epsilon) x^*_i$, where $0\leq \epsilon<1$. Consider the smooth function
\[\phi(\epsilon) = B_i \log u_i((1-\epsilon)x^*_i) - \langle p^*, (1-\epsilon)x^*_i\rangle = B_i \log u_i(x^*_i) - \langle p^*, x^*_i\rangle + B_i \log (1-\epsilon) + \epsilon \langle p^*, x^*_i\rangle. \]
Note that $\phi$ is differentiable on $(0,1)$ and \[\phi'(\epsilon) = - \frac{B_i}{1-\epsilon} + \langle p, x^*_i\rangle.\] 
Since $\langle p^*, x^*_i\rangle >B_i$, we have $\phi'(0) >0$. In other words, replacing $x^*_i$ by $\tilde{x}_i$ with sufficiently small $\epsilon$ strictly decreases the value of $r_i(x_i)$, contradicting to the choice of $x^*$. Therefore, \[\langle x^*, x^*_i\rangle \leq B_i\] 
for all $i$. Completely analogously, we can also show that \[\langle p^*, x^*_i\rangle \geq B_i\] 
for all $i$. Therefore, for each buyer $i$, $x_i^*$ is feasible and depletes its budget $B_i$ under prices $p^*$. Hence, for any $x_i \in \RR^m_+$, $\langle p^*, x_i\rangle \leq B_i$, since $x^*_i$ maximizes $r_i(x_i, p^*)$, we have
\[ B_i \log u_i(x^*_i) -\langle p^*, x^*_i \rangle \geq B_i \log u_i(x_i) - \langle p^*, x_i\rangle. \]
Since $\langle p^*, x_i\rangle \leq B_i = \langle p^*, x^*_i\rangle$, the above implies
\[ B_i \log u_i(x^*_i) \geq B_i \log u_i(x_i). \]
Therefore, $u_i(x^*_i) \geq u_i(x_i)$. In other words, $x^*_i \in D_i(p^*)$ (buyer $i$ is optimal) for all $i$.
\paragraph{Market clearance} By the complementary slackness condition, for item $j$ such that $\sum_i x^*_{ij} < 1$, it must holds that $p^*_j = 0$, completing the proof.
\paragraph{Remark} We can also assign any leftover of item $j$ to any buyer $i$ without violating its budget constraint, in order to ``clear'' the market. Meanwhile, since $u_i$ is CCNH, it is also ``monotone'' in the following sense: for any $\alpha\geq 0$, 
\[ u_i(x^*_i + \alpha \mathbf{e}^j) \geq u_i(x^*_i) + \alpha u_i(\mathbf{e}^j) \geq 0. \]
In other words, buyer $i$'s optimality is not affected by the assignment of any leftover of any item whose price is zero.

\subsection{Characterizations of Hoffman constant}\label{app:char-hoffman-const}
We compare our definition of Hoffman constant and another common, explicit characterization. Recall that $H_\cX(A)$ is the smallest $H$ such that, for any $b$, $\cS = \{x: Ax=b\}$,
\[ \|x-\Pi_{\cX\cap \cS}(x)\|\leq H\|Ax-b\|,\ \forall\, x\in \cX. \]

For any matrix $M$, let $\cB(M)$ be the set of \textit{nonsingular} submatrices consisting of rows of $M$. Define
\begin{align}
H(M) = \max_{B\in \mathcal{B}(M)} \frac{1}{\sigma_{\min}(B)} < \infty.  \label{eq:hoffman-unconstrained}
\end{align}
The following fact is known (see, e.g., \cite[\S 11.8]{guler2010foundations} and \cite[\S 2.1]{beck2017linearly}).
\begin{lemma}
	Suppose the reference polyhedral set can be represented by inequality constraints $\cX = \{x: Cx\leq d\}$.  Then,
	\[ H_{\cX}(A) \leq H\left( \begin{bmatrix}
	A \\ C
	\end{bmatrix} \right). \]
\end{lemma}
Clearly, $H(M)$ is finite for any $M$. In fact, this is the most well-known characterization of Hoffman constant, and is tight in the following sense: let $\cS = \{x: Ax=b\}$ for some arbitrary right hand side $b$, then it is the smallest constant $H$ such that  
\[ \|x - \Pi_{\cX \cap \cS}(x)\| \leq H \left\| \begin{bmatrix}
Ax - b \\ (Cx-d)_+
\end{bmatrix} \right\|\]
for all $x$ (not necessarily $\in \cX$). However, for all of our purposes, that is, analysis of PG, $x$ is always restricted to be $\in \cX$. Therefore, we choose to define $H_\cX(A)$ as such, consistent with \cite{beck2003mirror} and \cite{pena2020new}.
Meanwhile, the following is clear. 
\begin{lemma}
	For any matrices $A\in \RR^{m\times n}$, $m\leq n$ and $C\in \RR^{\ell \times n} $, it holds that \[
	H\left(\begin{bmatrix}
	A \\ C
	\end{bmatrix} \right) \geq \max \left\{\frac{1}{\sigma_{\min}(A)}, H(A) \right\}. \label{lemma:compare-H-H(A)-sigma_min}\]
\end{lemma}
\textit{Proof.} By definition \eqref{eq:hoffman-unconstrained}, $H' := H\left(\begin{bmatrix}
A \\ C^\top
\end{bmatrix} \right) \geq H(A)$. If ${\rm rank}(A) = m$, then $H' \geq \frac{1}{\sigma_{\min}(A)}$ because $A\in \cB\left( \begin{bmatrix}
A \\ C
\end{bmatrix} \right)$. If $r = {\rm rank}(A) < m$, let the (nonzero) singular values of $A$ be $\sigma_1 \geq \dots \geq \sigma_r = \sigma_{\min}(A) > 0$. Consider any $B\in \cB(A) \subseteq \cB\left( \begin{bmatrix}
A \\ C
\end{bmatrix} \right)$ with rank $r$ (having exactly $r$ rows), let its nonzero singular values be $\sigma'_1\geq \dots \geq \sigma'_r = \sigma_{\min}(B) > 0$. Applying Cauhchy's Interlacing Theorem (see, e.g., \cite[Theorem 1]{fan1957imbedding}) on $AA^\top$ and its principal submatrix $BB^\top$, we have
\[ \sigma_1 \geq \sigma'_1 \geq \dots \geq \sigma_r \geq \sigma_r'. \]
Therefore, $H' \geq \frac{1}{\sigma_{\min}(B)} \geq \frac{1}{\sigma_{\min}(A)}$.
\qed

\subsection{Proof of Theorem \ref{thm:h(Ax)+poly_ind}} \label{app:proof-h(Ax)+poly_ind-conv}
We follow the development in \cite[\S 4 \& Appendix F]{karimi2016linear} and further articulate the constants. There, the authors show that proximal gradient achieves linear convergence under the so-called \textit{Proximal-\PL} inequality. Consider the following general nonsmooth problem
\begin{align}
F^* = \min_x F(x) = f(x) +g(x) \label{eq:min-F=f+g}
\end{align} where $f$ is smooth convex with $L_f$-Lipschitz continuous gradient, $g$ is simple closed proper convex and $\dom g \subseteq \dom f$.
One iteration of the proximal gradient method with stepsize $\gamma > 0$ is as follows:
\begin{align}
x^{t+1} = \prox_{g}\left( x^t - \gamma \nabla f(x^t) \right) = \argmin_x \left[ \langle \gamma \nabla f(x), x-x^t\rangle + \frac{1}{2}\|x - x^t\|^2 + g(x)\right]. \label{eq:pg-iter}
\end{align}

For any $\alpha>0$ and any $x\in \dom g$, define
\begin{align}
\cD(x, \alpha) =  -2 \alpha \min_{x'} \left[ \langle \nabla f(x), x' - x \rangle  + \frac{\alpha}{2}\|x'-x\|^2 + g(x') - g(x) \right].
\label{eq:def-D_g(x,a)}
\end{align}
Say that $F = f+g$ satisfies the {proximal-\PL} inequality at $x$ w.r.t. $\Lambda \geq \lambda > 0$ if
\begin{align}
\frac{1}{2} \cD(x, \Lambda) \geq \lambda(F(x) - F^*),  \label{eq:prox-PL}
\end{align}


Below is essentially \cite[Theorem 5]{karimi2016linear}, which shows that the so-called {Proximal-\PL} condition is sufficient for linear convergence. Note that, different from \cite[Theorem 5]{karimi2016linear}, we only require \eqref{eq:prox-PL} to hold for $x\in \cX$ such that $F(x)\leq F(x^0)$ instead of all $x\in \cX$. In addition, we note that in some cases \eqref{eq:prox-PL} may hold with $\Lambda > L_f$, in which case the rate needs to be slightly adjusted. Since $\cD(x, \cdot)$ is monotone \cite[Lemma 1]{karimi2016linear}, \eqref{eq:prox-PL} holds when $\Gamma$ is replaced by $\Gamma' \geq \Gamma$. The statement and proof are the same as \cite[pp. 9]{karimi2016linear} otherwise.
\begin{theorem}
	Let $x^0 \in \dom g$. If $f$ and $g$ satisfies \eqref{eq:prox-PL} for all $x\in \dom g$ such that $F(x) \leq F(x^0)$, then $x^t$ defined by \eqref{eq:pg-iter} starting from $x^0$ with constant stepsize $\gamma = 1/L_f$ converges linearly with rate $1 - \frac{\lambda}{\bar{L}}$, where $\bar{L} = \max\{\Lambda, L_f \}$. In other words,
	\[ F(x^t) - F^* \leq \left( 1- \frac{\lambda}{\bar{L}}\right)^t (F(x^0) - F^*), \ t = 1, 2, \dots \]
	\label{thm:pg-lin-conv}
\end{theorem}
\textit{Proof.} By assumption, \eqref{eq:prox-PL} holds for all $x\in \dom g$, $x\leq F(x^0)$. In particular, it holds for $x^t$, $t, 1, 2, \dots$, since proximal gradient is a \textit{descent} method, i.e., $F(x^0) \geq F(x^1) \geq \dots$ (see, e.g., \cite[Corollary 10.18]{beck2017first}). 


Therefore, by $L_f$-Lipschitz continuity of $\nabla f$, proximal gradient update \eqref{eq:pg-iter}, definition of $D(x, \cdot)$, its monotonicity, and \eqref{eq:prox-PL} for all $x^t$,
\begin{align*}
F(x^{t+1}) 
& \leq F(x^t) + \langle \nabla f(x^t), x^{t+1} - x^t\rangle  + \frac{L_f}{2} \|x^{t+1} - x^t\|^2 + g(x^{t+1}) - g(x^t) \\
&\leq F(x^t) + \left[\langle \nabla f(x^t), x^{t+1} - x^t\rangle  + \frac{\bar{L}}{2} \|x^{t+1} - x^t\|^2 + g(x^{t+1}) - g(x^t) \right] \\
& \leq F(x^t) - \frac{1}{2\bar{L}} \cD(x^t, \bar{L}) \\
& \leq F(x^t) - \frac{\lambda}{\bar{L}}\left( F(x^t) - F^* \right) \\
\Rightarrow\ \ 	& F(x^{t+1}) - F^* \leq \left( 1 - \frac{\lambda}{\bar{L}} \right)\left( F(x^t) - F^* \right).
\end{align*}
Repeatedly applying the above inequality completes the proof. \qed


Then, we prove Theorem \ref{thm:h(Ax)+poly_ind}. Clearly, problem \eqref{eq:std-form-proj-grad-desc} is \eqref{eq:min-F=f+g} with $g(x)=\delta_\cX(x)$ and PG is a special case of proximal gradient. By Theorem \ref{thm:pg-lin-conv}, in order to prove Theorem \ref{thm:h(Ax)+poly_ind}, it suffices to establish the {Proximal-\PL} condition \eqref{eq:prox-PL} (for all $x\in \cX, f(x) \leq f(x^0)$ for some initial iterate $x^0$). Let $\cX^*$ be the set of optimal solutions to \eqref{eq:std-form-proj-grad-desc} and $f^*$ be the optimal objective value.  Since $h$ is $\mu$-strongly convex and $f(x) = h(Ax)$, there exists $z^*\in \dom f$ such that $\mathcal{S} = \{x: A x = z^*\}$ and $\cX^* = \cX \cap \mathcal{S}$. Therefore, for any $x\in \cX$, $x_p := \Pi_{\cX^*}(x)$, we have
\[ f(x_p) = h(A x_p) \geq h(A x) + \langle \nabla h(A x), A(x_p - x)\rangle + \frac{\mu}{2}\|A(x_p-x)\|^2. \]
Note that 
\[ \langle \nabla h(A x), A(x_p - x)\rangle =  \langle A^\top \nabla h(A x ), x_p - x\rangle = \langle \nabla f (x), x_p - x\rangle. \]


Hence, for any $x\in \cX$, by strong convexity of $h$ and definition of $H = H_\cX(A)$, we have
\begin{align*}
f(x_p) & \geq f(x) + \langle \nabla f(x), x_p - x\rangle + \frac{\mu}{2}\|A(x-x_p)\|^2 \\
& = f(x) + \langle \nabla f(x), x_p - x\rangle + \frac{\mu}{2}\|A x- z^*\|^2 \\

& \geq f(x) + \langle \nabla f(x), x_p - x\rangle + \frac{\mu}{2 H^2} \| x - x_p \|^2, 
\end{align*}

Therefore, 
\begin{align*}
f^* &\geq f(x) + \langle \nabla f(x), x_p - x\rangle + \frac{\mu}{2 H^2}\|x-x_p\|^2 \\
&\geq f(x) + \min_{y\in \cX} \left\{ \langle \nabla f(x), y-x\rangle + \frac{\mu}{2 H^2} \|y-x\|^2 \right\} \\
& \geq f(x) - \frac{H^2}{2\mu} \cD\left(x, \frac{\mu}{H^2}\right)  \\
\Rightarrow\ \ & \frac{1}{2} \cD\left(x, \frac{\mu}{H^2}\right) \geq \frac{\mu}{H^2}(f(x) - f^*).
\end{align*}


Thus, \eqref{eq:prox-PL} holds for all $x\in \cX$, $f(x)\leq f(x^0)$ with \[\Lambda = \lambda = \frac{\mu}{H^2}.\] 
Since $\nabla f(x) = A^\top \nabla h(A x)$ and $h$ is $(\mu, L)$-s.c., its Lipschitz constant can be chosen as 
\[L_f = L \|A\|^2.\] 

By Theorem \ref{thm:pg-lin-conv}, PG with stepsize $\gamma = \frac{1}{L_f}$ converges linearly with rate

\[ 1 - \frac{\frac{\mu}{H^2}}{\max\left\{ \frac{\mu}{H^2}, L\|A\|^2 \right\}} = 1 - \frac{\mu}{\max\{ \mu, LH^2 \|A\|^2  \}}.  \]
Finally, convergence of the distance to optimality $\|x^t - \Pi_\cX(x^t)\|$ is straightforward: for any $x\in \cX$, by the strong convexity of $h$ and definition of $H$,
\[ f(x) - f^* = h(Ax) - h(Ax_p) \geq \frac{\mu}{2}\|Ax - Ax_p\|^2 = \frac{\mu}{2}\|Ax-z^*\|^2 \geq \frac{\mu}{2 H} \|x - x_p\|^2. \] \qed

\paragraph{Remark}
A special case is when $d\geq r$ (recall that $A\in \RR^{d\times r}$) and ${\rm rank}(A) =r$. In this case, $f(x) = h(Ax)$ itself is strongly convex with modulus $\mu \sigma_{\min}(A)^2$. In this case, classical analysis (e.g., \cite[\S 10.6]{beck2017first}) implies linear convergence with rate $1 - \frac{\mu \sigma_{\min}(A)^2}{L \|A\|^2}$. Meanwhile, in the above analysis, we have $\cX^* = \{x^*\} = \cS = \{x: Ax = z^*\} = \cX \cap \cS$ (since $x^*, z^*$ are unique and ${\rm rank}(A)=r$). By Lemma \ref{lemma:hoffman-const-subspace}, for any $x$, it holds that
\[\|x - \Pi_{\cX^*}(x) \|\leq \frac{1}{\sigma_{\min}(A)^2}\|Ax-z^*\|.\]
Therefore, by the definition of Hoffman constant, $H_\cX(A)\leq \frac{1}{\sigma_{\min}(A)^2}$ and the classical rate under strong convexity is recovered.

\subsection{Proof of Theorem \ref{thm:h(Ax)+<q,x>-lin-conv}} \label{app:proof-h(Ax)+<q,x>+poly_ind-conv}
Let $\cX^*$ be the set of optimal solutions to \eqref{eq:std-form-proj-grad-desc}. First, recall the following lemma \cite[Lemma 14]{wang2014iteration}, which ensures the first part of the theorem, that is, uniqueness of $Ax^*$ and $q^\top x^*$ for all $x^*\in \cX^*$.
\begin{lemma}
	There exist unique $z^*\in \RR^r$ and $w^*\in \RR$ such that for any $x^*\in \cX^*$, 
	\[ Ax^* = z^*,\ \langle q^*, x\rangle = w^*. \]
	\label{lemma:Ax^*=z^*,qTx^*=w^*}
\end{lemma}
The next lemma is essentially \cite[Lemma 2.5]{beck2017linearly}. Different from the statement of \cite[Lemma 2.5]{beck2017linearly}, we keep $\|\nabla h(z^*)\|$ instead of bounding it by $\sup_{x\in \cX} \|\nabla h(Ax)\|$. We also define $C = f(x^0) - f^*$ instead of $C = \sup_{x\in \cX} f(x) - f^*$, since subsequent application of the lemma only involves PG iterates $x^t$, which have monotone decreasing objective values $f(x^0)\geq f(x^1)\geq \dots$ The proof remains unchanged otherwise.
\begin{lemma}
	Let $z^*$ be as in Lemma \ref{lemma:Ax^*=z^*,qTx^*=w^*} and $x^0 \in \cX$. For any $x\in \cX$ such that $f(x) \leq f(x^0)$, it holds that
	\begin{align*}
	\| x - \Pi_{\cX^*}(x)\|^2 \leq \kappa \left( f(x) - f^* \right),
	
	\end{align*}
	where, same as in Theorem \ref{thm:h(Ax)+<q,x>-lin-conv},
	$\kappa = H_\cX(A)^2 \left( C + 2 G D_A + \frac{2(G^2 + 1)}{\mu} \right)$, 
	
	$C = f(x^0) - f^*$,\newline 
	$G = \|\nabla h(z^*)\|$, $D_A = \sup_{x, y\in \cX} \|A(x - y)\|$.
	\label{lemma:quad-growth-h(Ax)+<q,x>}
\end{lemma}


Finally, take $L_f = L \|A\|^2$ as a Lipschitz constant of $\nabla f$. By Lemma \ref{lemma:quad-growth-h(Ax)+<q,x>} and \cite[\S 4.1]{karimi2016linear}, it holds that \eqref{eq:std-form-proj-grad-desc} satisfies the {proximal-\PL} inequality \eqref{eq:prox-PL} with \[\Lambda = \lambda = \frac{1}{2\kappa}\] 
for all $x\in \cX$ such that $f(x) \leq f(x^0)$ (in particular, for all $x^t$, $t = 1,2,\dots$). 


By Theorem \ref{thm:pg-lin-conv}, PG converges linearly with rate $1 - \frac{\lambda}{\max\{ \Lambda, L_f \}} = 1 - \frac{1}{\max\{ 1, 2\kappa L \|A\|^2 \}}$.


\paragraph{Remark} Lemma \ref{lemma:quad-growth-h(Ax)+<q,x>} shows that QG holds. Similar convergence guarantees can also be derived from other QG-based analysis, e.g., \cite[Corollary 3.7]{drusvyatskiy2018error}.

\subsection{Linear convergence of PG with linesearch} \label{app:pg-ls-details}
First, we consider the more general proximal gradient setup \eqref{eq:min-F=f+g}. Let $L_f$ be a Lipschitz constant of $\nabla f$ and the {Proximal-\PL} inequality \ref{eq:prox-PL} holds with $\Lambda \geq \lambda \geq 0$ for all $x\in \dom g$ such that $F(x)\leq F(x^0)$. 
Let $\alpha\geq 1$, $\beta\in (0,1)$, $\Gamma >0$ (increment factor, decrement factor, upper bound on stepsize, respectively). The linesearch subroutine $\mathcal{LS}_{\alpha, \beta, \Gamma}$ is defined in Algorithm \ref{alg:ls-subroutine}. 

\begin{algorithm}
	
	\caption{$x_{t+1}, \gamma_t, k_t\leftarrow \mathcal{LS}_{\alpha,\beta, \Gamma}(x, \gamma, k_{\rm prev})$ with parameters $\alpha\geq 1$, $\beta\in (0,1)$, $\Gamma>0$.}
	\label{alg:ls-subroutine}
	\begin{algorithmic}
		\STATE If $k_{\rm prev} = 0$, set $\gamma^{(0)} = \min\{ \alpha\gamma, \Gamma\}$. Otherwise, set $\gamma^{(0)} = \gamma$. 
		\STATE For $k=0, 1, 2, \dots$
		\begin{enumerate}
			\item Compute $x^{(k)} = \prox_{\lambda^{(k)} g}(x - \gamma^{(k)}\nabla f(x) )$.
			\item Break if 
			\begin{align}
			f(x^{(k)}) \leq f(x) + \langle \nabla f(x), x^{(k)} - x\rangle + \frac{1}{2 \gamma^{(k)}} \| x^{(k)} - x\|^2.  \label{eq:ls-break-condition}
			\end{align}
			\item Set $\gamma^{(k+1)} = \beta \gamma^{(k)}$ and continue to $k+1$. 
		\end{enumerate}
		Return $x_{t+1} = x^{(k)}$, $\gamma_t = \gamma^{(k)}$, $k_t = k$.
	\end{algorithmic}
\end{algorithm}

In this way, proximal gradient with linesearch can be described formally as follows: starting from $x^0 \in \dom f$, $\gamma_{-1} = \Gamma$, $k_{-1} = 0$, perform the following iterations
\begin{align*}
(x^{t+1}, \gamma_t, k_t) \leftarrow \mathcal{LS}_{\alpha, \beta, \Gamma}(x^t, \gamma_{t-1}, k_{t-1}), \ \ t = 1, 2, \dots 

\end{align*}


Note that \eqref{eq:ls-break-condition} holds for any $\gamma^{(k)} \leq \frac{1}{L_f}$ (see, e.g., \cite[Theorem 10.16]{beck2017first}). Therefore, Algorithm $1$ terminates when $ \gamma^{(0)} \beta^k \leq \frac{1}{L_f}$. This means 
\begin{align}
\gamma_t \geq \tilde{\gamma} :=  \min\left\{ \Gamma, \frac{\beta}{L_f}\right\}. \label{eq:gamma-min-def}
\end{align}
for all $t$. Note that we explicitly include the case of $\Gamma \leq \frac{1}{L_f}$, although in practice $\Gamma$ is often set very large. Clearly,
\[ \Gamma \beta^k \leq \tilde{\gamma}\ \Leftrightarrow \ k \geq \frac{\log \frac{\Gamma}{\tilde{\gamma}}}{\log \frac{1}{\beta}}. \]
Therefore, in Algorithm \ref{alg:ls-subroutine}, the backtracking iteration index satisfies $k_t \leq \frac{\log \frac{\Gamma}{\tilde{\gamma}} }{\log \frac{1}{\beta}}$ for all $t$. Note that if the loop breaks at $k_t$, the number of $\prox$ evaluations is exactly $k_t + 1$.


Let 
\begin{align}
\bar{L} = \max\left\{\frac{1}{\tilde{\gamma}}, \Lambda\right\} = \max\left\{ \frac{1}{\Gamma}, \frac{L_f}{\beta}, \Lambda \right\}. \label{eq:def-bar-L-triple-max}
\end{align}
Then, monotonicity of $D(x, \cdot)$ implies, for all $x\in \dom g$ such that $F(x)\leq F(x^0)$, 
\[\frac{1}{2} \cD(x, \bar{L}) \geq \frac{1}{2}\cD(x, \Lambda) \geq \lambda\left( F(x) - F^* \right). \]

Following the proof of Theorem \ref{thm:pg-lin-conv} (or that of \cite[Theorem 5]{karimi2016linear}), we have
\begin{align*}
F(x^{t+1}) & \leq F(x^t) + \langle  \nabla f(x^t), x^{t+1}-x^t\rangle + \frac{L_f}{2}\|x^{t+1}-x^t\|^2 + g(x^{t+1}) - g(x^t) \\
& \leq F(x^t) + \langle  \nabla f(x^t), x^{t+1}-x^t\rangle + \frac{\bar{L}}{2}\|x^{t+1}-x^t\|^2 + g(x^{t+1}) - g(x^t) \\
& \leq F(x^t) - \frac{1}{2\bar{L}}\cD\left(x^t, \bar{L}\right) \\
& \leq F(x^t) - \frac{\lambda}{\bar{L}}(F(x^t) - F^*) \\
\Rightarrow \ \ & F(x^{t+1}) - F^* \leq \left(1 - \frac{\lambda}{\bar{L}}\right) \left( F(x^t) - F^* \right). 
\end{align*}
Summarizing the above discussion, we have the following convergence guarantee for PG with linesearch.
\begin{theorem}
	Let $\alpha\geq 1$, $\beta\in (0,1)$ and $\Gamma>0$. For problem \eqref{eq:min-F=f+g} satisfying the Proximal-\PL inequality with $\Lambda \geq \lambda > 0$ for all $x\in \dom g$ such that $F(x)\leq F(x^0)$, proximal gradient \eqref{eq:pg-iter} with linesearch subroutine $\mathcal{LS}_{\alpha, \beta, \Gamma}$ described in Algorithm \ref{alg:ls-subroutine} generates iterates $x^t$ such that
	\begin{align}
	F(x^{t+1}) - F^* \leq \left(1 - \frac{\lambda}{\bar{L}} \right)^t \left( F(x^0)  - F^*\right),\ \ t = 1, 2, \dots, \label{eq:prox-grax-f+g-with-ls-lin-conv}
	\end{align}
	where $\bar{L}$ is defined in \eqref{eq:def-bar-L-triple-max}. Furthermore, each iteration requires at most $1 + \frac{\log \frac{\Gamma}{\tilde{\gamma}}}{\log \frac{1}{\beta}}$ number of $\prox$ evaluations.
	\label{thm:prox-grad-f+g-with-ls-lin-conv}
\end{theorem}

\textit{Proof of Theorem \ref{thm:pg-ls-conv}.} In the above discussion, when $g(x) = \delta_\cX(x)$, we can replace the Lipschitz constant $L_f$ by the restricted one $L_f^\cX$ throughout, since Algorithm \ref{alg:ls-subroutine} ensures $x^t \in \cX$ for all $t$. It remains to apply Theorem \ref{thm:prox-grad-f+g-with-ls-lin-conv}. For $q = 0$, $\Lambda = \lambda = \frac{\mu}{H^2}$ and $\bar{L} = \max\left\{ \frac{1}{\Gamma}, \frac{L_f^\cX}{\beta}, \frac{\mu}{H^2} \right\}$. Therefore, the rate is 
\[ 1 - \frac{\lambda}{\bar{L}} = 1 - \frac{\mu}{\max\{\mu,  H^2/\Gamma,  H^2L_f^\cX/\beta  \}}. \]
For $q\neq 0$, $\Lambda = \lambda = \frac{1}{2 \kappa}$ and $\bar{L} = \max\left\{ \frac{1}{\Gamma}, \frac{L_f^\cX}{\beta}, \frac{1}{2\kappa} \right\}$. Therefore, the rate is
\[1 - \frac{1}{\max\{1, 2\kappa L_f^\cX /\beta, 2\kappa/\Gamma \}}. \]
\qed

\subsection{Other utility functions}\label{app:CD-and-CES}

Recall that, by Theorem \eqref{thm:eg-gives-me-for-certain-ui}, for any CCNH utilities $u_i$, optimal solutions to the EG convex program \eqref{eq:eisenberg-gale-primal} correspond to equilibrium allocation and prices.

\textbf{CES utilities} are parametrized by a nondegenerate $v$ and exponent $\rho\in (-\infty, 1]\backslash \{ 0\}$:
\[ u_i(x_i) = \left(\sum_{j=1}^m v_{ij}x_{ij}^\rho \right)^{1/\rho}.  
\]
Clearly, $\rho=1$ gives linear utilities. For $\rho < 1$, it has been shown that Proportional Response dynamics achieves linear convergence in prices and utilities \cite[Theorem 4]{zhang2011proportional} under their notion of $\epsilon$-approximate market equilibrium \cite[pp. 2693]{zhang2011proportional}. 


\textbf{Cobb-Douglas utilities} represent substitutive items and take the following form, for parameters $\lambda = (\lambda_i)$, $\lambda_i \in \Delta_m$:
\[ u_i(x_i) = \Pi_j x_{ij}^{\lambda_{ij}}. \]  
In this case, EG \eqref{eq:eisenberg-gale-primal} decomposes item-wise into simple problems with explicit solutions. Specifically, for each item $j$, the minimization problem is 
\[\min_{x_{:, j} \in \Delta_n} - \sum_i B_i \lambda_{ij} \log x_{ij}. \]
Let $p_j$ be the Lagrangian multiplier associated with constraint $\sum_i x_{ij} = 1$. The Lagrangian is
\[ \mathcal{L}(x_{:, j}, p_j) = -\sum_i B_i \lambda_{ij} \log x_{ij} + p_j\left( \sum_i x_{ij} - 1\right). \]
By first-order stationarity condition, for any $p_j\in \RR$, $\mathcal{L}(x_{:, j}, p_j)$ is minimized when 
\begin{align}
	x_{ij} =\frac{B_i \lambda_{ij}}{p_j}. \label{eq:cobb-douglas-x-in-p}
\end{align}
Substituting it into $\mathcal{L}$ and discarding the constants w.r.t. $p_j$, we have
\[ g(p_j) = \left(\sum_i B_i \lambda_{ij}\right) \log p_j - p_j, \]
which is maximized at equilibrium prices
\[ p^*_j  = \sum_i B_i \lambda_{ij}. \]
Therefore, by \ref{eq:cobb-douglas-x-in-p}, the equilibrium $x^*$ under Cobb-Douglas utilities is given by
\[ x^*_{ij} = \frac{B_i \lambda_{ij}}{\sum_i B_i \lambda_{ij}},\ \forall\, i,j. \]


\section{Linear utilities}
\subsection{Shmyrev's convex program} \label{app:shmyrev-cp}
Under linear utilities, it turns out that we can also compute market equilibrium via the following convex program due to Shmyrev \cite{shmyrev2009algorithm, birnbaum2011distributed}. In this convex program, the variables are the \textit{bids} $b_{ij}$, $i\in [n]$, $j\in [m]$ and prices $p_j$, $j\in [m]$. 
\begin{align}
\max\, \sum_{i, j} b_{ij} \log v_{ij} - \sum_j p_j \log p_j\ \ {\rm s.t.}\ \sum_i b_{ij} = p_j,\ j \in [m],\ \sum_j b_{ij} = B_i,\ i\in [n], \ b\geq 0. \label{eq:shmyrev}
\end{align}
Given an optimal solution $b^*$, equilibrium prices and allocations are then given by $p^*_j = \sum_i b^*_{ij}$ and $x^*_{ij} = \frac{b^*_{ij}}{p^*_j}$, respectively.

\subsection{Proof of Lemma \ref{lemma:eg-linear-u_min-u_max}}\label{app:eg-linear-u_min-u_max}
Any $x\in \cX$ satisfies $x\leq 1$. Therefore, $\langle v_i, x_i\rangle \leq \|v_i\|_1 \|x^*_i\|_\infty \leq \|v_i\|_1 = \bar{u}_i$. For the lower bound, recall that at an equilibrium allocation $x^*$ ensures that every buyer gets at least the utility of the proportional share, that is, 
\[ \langle v_i, x^*_i\rangle \geq \left \langle v_i, \frac{B_i}{\|B\|_1} \ones\right \rangle = \frac{B_i\|v_i\|_1}{\|B\|_1} = \ubar{u}_i. \]


\subsection{Uniqueness of equilibrium quantities and convergence of $u^t$, $p^t$} \label{app:eg-convergence-of-simplex-prices}

Convergence of $u^t$ to $u^*$ can be easily seen as follows. Let $x^t$ be the PG iterates and $\tilde{h}$, $A$, $f = \tilde{h}(Ax)$, $\mu$ be defined as in \S \ref{sec:linear-utilities} and $f^* = \min_{x\in \cX} f(x)$. Since $\tilde{h}$ is $\mu$-strongly convex, we have
\[ \frac{\mu}{2}\|u^t - u^*\|^2\leq \tilde{h}(u^t) - \tilde{h}(u^*) \leq \tilde{h}(Ax^t) - f^*, \]
which converges linearly. Next, we show uniqueness of $p^*$ via simple arguments and construct a sequence of linearly convergent prices $p^t$.

\begin{lemma}
	Assume that $v$ is nondegenerate. Then, the equilibrium prices $p^*$ under linear utilities are unique. \label{lemma:unique-p*}
\end{lemma}
\textit{Proof.} By Theorem \ref{thm:eg-gives-me-for-certain-ui} and \cite[Lemma 3]{cole2017convex}, $p^*$ is an optimal solution (together with some $\beta^*$) to the following problem (dual of \eqref{eq:eisenberg-gale-primal} with linear utilities): \qed
\begin{align}
	 \min_{p,\, \beta}\ \sum_j p_j - \sum_i B_i \log \beta_i \ \ {\rm s.t.}\  p \geq 0,\, \beta\geq 0,\, p_j \geq v_{ij}\beta_i,\, \forall\, i,j.  \label{eq:eg-dual-p-beta}
\end{align}
Here, strong duality holds since there clearly exist primal and dual strictly feasible solutions with finite objective values given nondegenerate $v$ (c.f. Theorem \ref{thm:eg-gives-me-for-certain-ui} and Appendix \ref{app:thm-eg-and-mkt-equi}). We can eliminate $p$ by letting $p_j = \max_i v_{ij}\beta_i$ for all $j$ and rewrite \eqref{eq:eg-dual-p-beta} as 
\[ \min_\beta\, \sum_j \max_i v_{ij}\beta_j - \sum_i B_i \log \beta_i\ \ {\rm s.t.}\ \beta\geq 0. \]
In the above, since the objective is strongly convex and the feasible region is $\beta\geq 0$, the optimal solution $\beta^*$ is clearly unique. Furthermore, it must hold that $\beta^* >0$ (since the optimal objective value is finite and strong duality holds). For $p^*$ optimal to \eqref{eq:eg-dual-p-beta}, it must hold that $p^*_j = \max_i v_{ij}\beta^*_i$. In fact, $p^*_j \geq \max_i v_{ij}\beta^*_i$ by feasibility and, for any strict inequality, decreasing the corresponding $p^*_j$ strictly decreases the objective. \qed

The following lemma provides simple upper and lower bounds on feasible and equilibrium prices, respectively. The lower bounds are slightly strengthened over the existing one \cite[Lemma 17]{birnbaum2011distributed}.
\begin{lemma}
	Let $p^*$ be equilibrium prices under linear utilities with nondegenerate valuations $v$. Then, $\ubar{p}_j \leq p^*_j \leq \bar{p}_j$ for all $j$, where $\ubar{p}_j = \max_i \frac{v_{ij} B_i}{\|v_i\|_1}$ and $\bar{p}_j = \|B\|_1$.  \label{lemma:linear-pmin-pmax}
\end{lemma}
\textit{Proof.} It is essentially the same as the proof of Lemma \ref{lemma:ql-shmyrev-upper-lower-bounds}, except that, at optimality, $u_i \leq \|v_i\|_1 + B_i$ can be strengthened to $u_i \leq \|v_i\|_1$ (utility of each buyer is at most that of having a unit of every item). \qed


\paragraph{A linearly convergent sequence of $p^t$.} Here, all norms are vector norms. Note that each step of PG is of the form $x^{t+1} = \Pi_\cX(\bar{x}^t)$, where $\bar{x}^t = x^t - \gamma \nabla f(x^t)$. Since $\nabla f$ is $L_f$-Lipschitz, the mapping \[\phi_1: x \mapsto x- \gamma \nabla f(x^t)\] 
is Lipschitz continuous (w.r.t. $\|\cdot\|_2$) with constant $1 + \gamma L_f = 2$ (where $\gamma = \frac{1}{L \|A\|^2}$ is the fixed stepsize). Meanwhile, we have the following. 
\begin{lemma}
	Let $y \in \RR^n$ and $y^* = \Pi_{\Delta^n}(y)$. There exists a unique multiplier $\lambda \in \RR$, which can be computed in $O(n \log n)$ time, such that
	\begin{align}
	\sum_{i=1}^n (y_i - \lambda)_+ = 1. \label{eq:simplex-unique-multi}
	\end{align}
	Moreover, the mapping $\phi_2: y \mapsto \lambda$ is 
	piecewise linear and 
	$1$-Lipschitz continuous w.r.t. $\|\cdot \|_1$. \label{simplex-projection-multiplier-Lipschitz}
\end{lemma}
\textit{Proof.} 
By the KKT conditions for simplex projection (see, e.g., \cite[\S 3]{wang2013projection}), it holds that there exists unique $\lambda$ such that
\[ y^* = (y - \lambda \ones)_+.  \]

Suppose there exists $\lambda_1 < \lambda_2$ that satisfy \eqref{eq:simplex-unique-multi}. Then, since the left-hand side of \eqref{eq:simplex-unique-multi}, denoted as $w(\lambda)$, is monotone decreasing in $\lambda$, it must hold that $w(\lambda)=1$ for all $\lambda\in [\lambda_1, \lambda_2]$. In other words, $w(\cdot)$ is \textit{constant} on $[w_1, w_2]$. This further implies $w(\lambda) = 0$ for all $w\in [w_1, w_2]$, a contradiction. Therefore, $\lambda = \phi_2(y)$ is uniquely defined. Let $I^+(y)$, $I^0(y)$, $I^-(y)$ denote the set of indices $i\in [n]$ such that $y_i > \lambda$, $y_i = \lambda$, $y_i < \lambda$, respectively (where $\lambda = \phi_2(y)$). We have 
\[ \lambda = \frac{\sum_{i\in I^+(y)} y_i - 1}{ |I^+(y)|} = \frac{\sum_{i\in I^+(y)\cup I^0(y)} y_i - 1}{|I^+(y)| + |I^0(y)|}, \]
which is piecewise linear in $y$ since there are only finitely many index possible sets of indices and $I^+(y)$ is always nonempty (otherwise $\sum_i (y_i - \lambda)_+ = 0$). To see Lipschitz continuity, let $y'$ be such that $\|y'-y\|_1 \leq \epsilon$, where $0<\epsilon < \min \{|y_i - y_j|: i,j\in [n],\, y_i \neq y_j \}$. It must hold that $I^+(y)\subseteq I^+(y')$. In other words, $\lambda' = \phi_2(y')$ does not deactivate any $i\in I^+(y)$, only bringing new $i\in I^0(y)$. Hence, it holds that $|\lambda' - \lambda| \leq \frac{\|y - y'\|_1}{|I^+(y)|} \leq \|y-y'\|_1$. In other words, $\phi_2$ is $1$-Lipschitz continuous w.r.t. $\|\cdot \|_1$. 

Finally, \cite[Algorithm 1]{wang2013projection}) computes $\lambda$ and $y^*$ in $O(n \log n)$ time. 
\qed

Abusing the notation, let $\phi_2$ also denote the mapping from $x\in \RR^{n\times m}$ to $\lambda\in \RR^m$, that is, $\lambda_j = \varphi_2(x_{1j}, \dots, x_{nj})$. 
Let \[\phi(x) = \phi_2 (\phi_1(x))/\gamma\] 
and $p^t = \phi(x^t)$. Here, $\phi_1$ is $2$-Lipschitz continuous and $\phi_2$ is $1$-Lipschitz continuous w.r.t. $\|\cdot \|_1$. For any optimal solution $x^* \in \cX^*$, by $x^* = \Pi_\cX (x^*)$ and KKT conditions for \eqref{eq:eisenberg-gale-primal} and \eqref{eq:eg-dual-p-beta}, it can be seen that 
\[p^* = \phi(x^*).\] 

Using the Lipschitz continuity properties of $\phi_1, \phi_2$ and Theorem \ref{thm:h(Ax)+poly_ind}, we have
\begin{align*}
\|p^t - p^*\|_1 &= \|\phi(x^t) - \phi( \Pi_{\cX^*}(x^t))\|_1 \leq \frac{1}{\gamma}\| \phi_1(x^t)  - \phi_1(\Pi_{\cX^*}(x^t)) \|_1
\\ & \leq \frac{n}{\gamma} \|\phi_1(x^t) - \phi_1(\Pi_{\cX^*}(x^t))\|  \leq \frac{2n}{\gamma} \cdot \|x^t - \Pi_{\cX^*}(x^t)\| \\
& \leq \frac{2n }{\gamma} \cdot \sqrt{\frac{2 H_\cX(A)}{\mu} \left( f(x^t) - f^* \right) } \\
& \leq \frac{2n}{\gamma} \sqrt{ \frac{2 H_\cX(A)}{\mu}}\cdot  \left( 1 - \frac{\mu}{2 H L \|A\|^2}\right)^{t/2} \cdot \sqrt{f(x^0) - f^*}. 

\end{align*}
Therefore, we can take $C = \frac{2n}{\gamma} \sqrt{ \frac{2 H_\cX(A)}{\mu}} \cdot \sqrt{f(x^0) - f^*}$ and $\rho = \sqrt{1 - \frac{\mu}{2 H L \|A\|^2}} \in (0,1)$. 


Since $p^* \geq \ubar{p} > 0$, we can bound the maximum relative price error $\eta^t = \max_j \frac{|p^t_j - p^*_j|}{p^*_j}$ as follows, where $\ubar{p}_{\min} = \min_j \ubar{p}_j$. 
\[ \eta^t \leq \frac{\|p^t - p^*\|_1}{\ubar{p}_{\min}} \leq \frac{C}{\ubar{p}_{\min}} \cdot \rho^t. \]
In other words, $\eta^t$ converges (R-)linearly to zero.  


\section{QL utilities}
\subsection{Derivation of the QL-Shmyrev convex program \eqref{eq:ql-shmyrev}} \label{app:ql-shmyrev-dual}
In \cite[Lemma 5]{cole2017convex}, the convex program for the equilibrium prices is as follows:
\begin{align}
\min \, \sum_j p_j - \sum_i B_i \log \beta_i \ \ 
{\rm s.t.} \ v_{ij}\beta_i \leq p_j,\ \forall\, i,j,\ \ 0 \leq \beta \leq 1. \label{eq:ql-dual-p-beta}
\end{align}
Note that it is simply the dual of EG under linear utilities \eqref{eq:eg-dual-p-beta} with additional constraints $\beta\leq 1$.
Assuming $v$ is nondegenerate, by a change of variable and Lagrange duality, we can derive the dual of \eqref{eq:ql-dual-p-beta}. First, at optimality, it must holds that $\beta_i>0$ for all $i$. Therefore, by nondegeneracy of $v$, $p_j > 0$ for all $j$ at optimality. Let $p_j = e^{q_j}$ and $\beta_i = e^{-\gamma_i}$. The above problem is equivalent to 
\begin{align}
\begin{split}
	\min \  & \sum_j e^{q_j} + \sum_i B_i \gamma_i \\
	\ {\rm s.t.}& \quad q_j + \gamma_i \geq \log v_{ij}, \ \forall\, i,j, \\
	& \quad \gamma \geq 0. 
\end{split} \label{eq:ql-dual-exponentiated}
\end{align}
Let $b_{ij} \geq 0$ be the dual variable associated with constraint $q_j + \gamma_i \geq \log v_{ij}$. The Lagrangian is
\begin{align*}
	L(q, \gamma, b) & := \sum_j e^{q_j} + \sum_i B_i \gamma_i - \sum_{i, j} b_{ij}\left( q_j + \gamma_i - \log v_{ij} \right) \\
	& =  \sum_j \left( e^{q_j} -\left(\sum_i b_{ij}\right) q_j \right) + \sum_i (B_i - \sum_j b_{ij}) \gamma_i + \sum_{i, j} (\log v_{ij}) b_{ij}. 
\end{align*}
Clearly, when $\sum_j b_{ij} \leq B_i$ for all $i$, $\gamma\geq 0$, $L(q, \gamma, b)$ is minimized at $q_j = \log \sum_i b_{ij}$ and $\gamma = 0$. When $\sum_j b_{ij} > B_i$ for some $i$, $L \rightarrow -\infty$ as $\gamma_i \rightarrow \infty$. Therefore, when $\sum_j b_{ij} \leq B_i$ for all $i$, we have
\[ g(b) = \sum_j \left[\sum_{i, j} b_{ij} - \left( \sum_i b_{ij}\right) \log \sum_i b_{ij}  \right] + \sum_{i, j} (\log v_{ij}) b_{ij}. \]
Therefore, the dual is
\[ \max\, g(b) \ \ {\rm s.t.}\ b\geq 0,\ \sum_j b_{ij} \leq B_i,\ \forall\, i. \]
Adding slack variables $\delta = (\delta_1 \dots, \delta_n)$ and writing it in minimization form yield \eqref{eq:ql-shmyrev}. 
\paragraph{Remark.} As mentioned in \S \ref{sec:ql-utilities}, when some $v_{ij} = 0$ (but $v$ is still nondegenerate), by the above derivation, the first summation in \eqref{eq:ql-shmyrev} should be replaced by $\sum_{(i,j)\in \mathcal{E}}$, where $\mathcal{E} = \{(i,j): v_{ij}>0\}$. The dual remains the same otherwise.


\subsection{Proof of Lemma \ref{lemma:ql-shmyrev-upper-lower-bounds}}\label{app:proof-ql-shmyrev-bounds}
Similar to the proof of Lemma \ref{lemma:unique-p*}, this can be seen via the uniqueness of the optimal solution $(p^*, \beta^*)$ of \eqref{eq:ql-dual-p-beta}, that is, from uniqueness of $\beta^*$ to that of $p^*_j = \max_i v_{ij}\beta^*_i$. 

Let $(b^*, \delta^*)$ be an optimal solution to \eqref{eq:ql-shmyrev}. Note that strong duality holds for \eqref{eq:ql-dual-exponentiated} and \eqref{eq:ql-shmyrev}, since there exit simple strictly feasible solutions. By the derivation in Appendix \ref{app:ql-shmyrev-dual}, it holds that $q^*_j = \log \sum_i b^*_{ij}$ gives an optimal solution to \eqref{eq:ql-dual-exponentiated} (the first-order optimality condition). Therefore, 
\[p^*_j = e^{q^*_j} = \sum_i b^*_{ij}. \]


Next we establish the upper and lower bounds on $p^*$. By the derivation in Appendix \ref{app:ql-shmyrev-dual} and Lagrange duality, for any optimal solution $b^*$ to \eqref{eq:ql-shmyrev}, it holds that $p^*_j:= \sum_i b^*_{ij}$ and $\beta^*_j = \min_{j\in J_i}\frac{p^*_j}{v_{ij}}$ give the (unique) optimal solution to \eqref{eq:ql-dual-p-beta}. Clearly, $\beta^*\leq 1$ and therefore
\[ p^*_j = \max_i v_{ij}\beta^*_i \leq \max_i v_{ij} = \bar{p}_j. \]
By \cite[Lemma 5]{cole2017convex}, the dual of \eqref{eq:ql-dual-p-beta} is (c.f. the original EG primal \ref{eq:eisenberg-gale-primal})
\begin{align}
\begin{split}
	 \max_{u,\, x,\, s} & \sum_i B_i \log u_i - s_i \\
	\st & u_i \leq v_i^\top x_i + s_i,\ \forall\, i,\ \\ 
	& \sum_i x_{ij} \leq 1,\ \forall\, j, \\
	&  x, s\geq 0.
\end{split} \label{eq:ql-eg-form}
\end{align}
Clearly, strong duality holds for \eqref{eq:ql-dual-p-beta} and \eqref{eq:ql-eg-form}. Furthermore, notice the following. 
\begin{itemize}
	\item $\beta^*_i = \frac{B_i}{u^*_i}$ at optimality, where $u^*_i$ is the amount of utility of buyer $i$. This is by the stationarity condition in the KKT optimality conditions.
	\item $u^*_i \leq \|v_i\|_1 + B_i$, where the right hand side is the amount of utility of all items and the entire budget. This can also be seen as follows. When $s_i > B_i$, decreasing $s_i$ strictly increases the objective of \eqref{eq:ql-eg-form}. Therefore, the optimal $s^*$ must satisfy $s^*_i \leq B_i$. It then follows from the constraint $u_i \leq v_i^\top x_i + s_i$. 
\end{itemize}
Therefore, 
\[ p^*_j \geq \max_i v_{ij} \beta^*_i \geq \max_i \frac{v_{ij} B_i }{\|v_i\|_1 + B_i} = \ubar{p}_j. \] 


\subsection{Proof of Theorem \ref{thm:pr-ql-shmyrev-1/T}} \label{app:proof-MD-1/T}
Similar to \cite[Lemma 7]{birnbaum2011distributed}, we first establish the following ``generalized Lipschitz condition'' for $\varphi$, which is key to the claimed last-iterate convergence.
\begin{lemma}
	For all $(b, \delta), (b', \delta')\in \cB$, it holds that 
	\begin{align}
	\varphi(b') \leq \varphi(b) + \langle \nabla \varphi(b), b' - b \rangle + D(b', \delta'\| b, \delta). \label{eq:generalized-Lipschitz}
	\end{align}
	\label{lemma:generalized-Lipschiz}
\end{lemma}
\textit{Proof.} 
Recall that $p_j(b) = \sum_i b_{ij}$, $\frac{\partial}{\partial b_{ij}} \varphi(b) = \log \frac{p_j(b)}{v_{ij}}$. For $(a, \delta^a), (b, \delta^b)\in \cB$, we have
\begin{align}
& \varphi(b) - \varphi(a) - \langle \nabla \varphi(a), b - a \rangle \nnnl
&= -\sum_{i,j} (1+\log v_{ij})(b_{ij} - a_{ij}) + \sum_j p_j(b)\log p_j(b) - \sum_j p_j(a)\log p_j(a) \nnnl
& \quad \quad - \sum_{i,j} (b_{ij} - a_{ij})\log \frac{p_j(a)}{v_{ij}} \nnnl
&= -\sum_{i,j} (b_{ij} - a_{ij}) + \sum_{j} p_j(b)\log \frac{p_j(b)}{p_j(a)} \nnnl
& = \sum_i (\delta^b_i - \delta^a_i) + \sum_j p_j(b)\log \frac{p_j(b)}{p_j(a)}. \label{eq:bounding-phi(b)-phi(a)-...}
\end{align}
Note that convexity and smoothness of $x \mapsto x\log \frac{x}{y}$ ($y>0$) implies
\begin{align}
\delta^b_i - \delta^a_i \leq \delta_i^b \log \frac{\delta^b_i}{\delta^a_i}. \label{eq:delta(b)-delta(a)<=...}
\end{align}
As in the proof of \cite[Lemma 7]{birnbaum2011distributed}, by convexity of $q(x, y) = x\log \frac{x}{y}$, it holds that
\begin{align}
\sum_j p_j(b)\log \frac{p_j(b)}{p_j(a)} \leq \sum_{i,j} b_{ij}\log \frac{b_{ij}}{a_{ij}}. \label{eq:D(p(b)||p(a))<=D(b||a)}
\end{align}
By \eqref{eq:delta(b)-delta(a)<=...} and \eqref{eq:D(p(b)||p(a))<=D(b||a)}, the right hand side of \eqref{eq:bounding-phi(b)-phi(a)-...} can be bounded by $D(b, \delta^b\| a, \delta^a)$. Therefore, \eqref{eq:generalized-Lipschitz} holds. \qed

Next, we prove the inequality on the right. Clearly, $(b^0, \delta^0)\in \cB$. By \cite[Theorem 3]{birnbaum2011distributed} (with objective $f = \varphi$, constraint set $C = \cB$ and stepsize $\gamma$), we have 
\[ \varphi(b^t) - \varphi(b^*) \leq \frac{D(b^*, \delta^*\|b^0, \delta^0)}{t}. \]
Similar to the proof of \cite[Lemma 13]{birnbaum2011distributed}, we can bound the Bregman divergence on the right hand side as follows, where $b_{ij} = \delta_i = \frac{B_i}{m+1}$. 
\begin{align*}
D(b^*, \delta^*\|b^0, \delta^0) &= \sum_{i,j} b^*_{ij} \log \frac{b^*_{ij} }{B_i} + \sum_i \delta^*_i \log \frac{\delta^*_i}{B_i} + \sum_{i,j} b^*_{ij}\log (m+1) + \sum_i \delta^*_i \log (m+1) \\ 
& \leq \sum_{i,j} b^*_{ij}\log (m+1) + \sum_i \delta^*_i \log (m+1) \\
& \leq \|B\|_1 \log (m+1),
\end{align*}
where the first inequality is because $\frac{b^*_{ij}}{B_i} \leq 1$. 
Combining the above yields the desired inequality.\footnote{In fact, the bound $\log(mn)$ in  \cite[Lemma 13]{birnbaum2011distributed} (which assumes $\|B\|_1 = 1$) can be easily strengthened to $\log m$ via the above derivation. In other words, it does not depend explicitly on the number of buyers (but implicitly through $\|B\|_1$ in general).}

Finally, we show the inequality on the left. By optimality of $(b^*, \delta^*)$, we have \[\langle \nabla \varphi(b^*), b - b^*\rangle \geq 0,\ \ \forall \, (b, \delta)\in \cB. \] 
Recall that $p_j(b) = \sum_i b_{ij}$. By \eqref{eq:bounding-phi(b)-phi(a)-...}, we have
\[  D(p^t\|p^*) = -\sum_{i,j} (b^t_{ij} - b^*_{ij}) + \sum_{j} p_j(b^t)\log \frac{p_j(b^t)}{p_j(b^*)}  \leq \varphi(b^t) - \varphi^*. \]\qed

\subsection{Details from MD \eqref{eq:md-udpate} to PR \eqref{eq:pr-ql-shmyrev-final}} \label{app:details-md-to-pr-ql}
Note that \eqref{eq:md-udpate} is buyer-wise separable: for each $i$, we have (where $\frac{\partial}{\partial b_{ij}} \varphi_b(b) = \log \frac{p_j(b)}{v_{ij}}$ and $\cB_i = B_i \cdot \Delta_{m+1}$)
\begin{align}
(b_i^{t+1}, \delta_i^{t+1}) & = \argmin_{(b_i, \delta_i)\in \cB_i}\, \sum_j \left( \log \frac{p_j(b^t)}{v_{ij}} - \log b^t_{ij}\right) b_{ij} - (\log \delta^t_i)\delta_i + \sum_j b_{ij} \log b_{ij} + \delta_i \log \delta_i \nonumber \\
& = \argmin_{(b_i, \delta_i)\in \cB_i}\, - \sum_j (\log b^t_{ij}) b_{ij} - (\log \delta^t_i)\delta_i + \sum_j b_{ij} \log b_{ij} + \delta_i \log \delta_i.
\label{eq:md-update-each-i}
\end{align}
By Lemma \ref{lemma:md-subproblem}, for all $i,j$, 
\begin{align}
b_{ij}^{t+1} = B_i\cdot \frac{\frac{v_{ij}b^t_{ij}}{p_j(b^t)}}{\sum_\ell \frac{v_{i\ell}b^t_{i\ell}}{p_\ell(b^t)} + \delta^t_i}
\label{eq:formula-b(t+1),delta(t+1)},\ \  \delta^{t+1}_j = B_i\cdot \frac{\delta^t_i}{\sum_\ell \frac{v_{i\ell}b^t_{i\ell}}{p_\ell(b^t)} + \delta^t_i}.
\end{align}
Let $p^t_j = p_j(b^t)$. Then, \eqref{eq:formula-b(t+1),delta(t+1)} can be written in terms of the allocations $x^t_{ij} = b^t_{ij} /p^t_j$ (which sum up to $1$ over buyers $i$ for any item $j$) and leftover $\delta^t_i$, thus giving \eqref{eq:pr-ql-shmyrev-final}. 


\subsection{Convergence of prices} \label{app:ql-price-conv}
Let $\eta^t = \max_j \frac{|p^t_j - p^*_j|}{p^*_j}$ be the relative price error, which can clearly be bounded by $\frac{\|p^t - p^*\|_1}{\ubar{p}_{\min}}$, where $\ubar{p}_{\min} = \min_j \ubar{p}_j >0$ is given in Lemma \ref{lemma:ql-shmyrev-upper-lower-bounds}. By Theorem \ref{thm:pr-ql-shmyrev-1/T} and strong convexity of KL divergence (w.r.t. $\|\cdot \|_1$), for $b^t$ and $p^t = p(b^t)$ generated by either PG or PR,
\begin{align}
	\frac{1}{2} \|p^t - p^*\|_1^2 \leq D(p^t\|p^*) \leq \varphi(b^t) - \varphi^*. \label{eq:price-1-norm-KL-obj-error}
\end{align}
Therefore, for PG, the quantities $\eta^t$, $\|p^t - p^*\|$ and $D(p^t\|p^*)$ all converge linearly to $0$. For PR, they converge at $O(1/T)$. 


We can further bound $\varphi(b^t) - \varphi^*$ by the duality gap. Specifically, given $b^t$, $p^t = p(b^t)$, let 
\[b^t_i = \min\left\{ \min_j \frac{p^t_j}{v_{ij}}, 1 \right\}.\] 
Then, $(p^t, \beta^t)$ is feasible to \eqref{eq:ql-dual-p-beta}. By weak duality,
\begin{align}
	\varphi(b^t) - \varphi^* \leq \varphi(b^t) + g(p^t, \beta^t), \label{eq:ql-shmyrev-dgap-expre}
\end{align}
where $g(p, \beta)$ is the (minimization) objective of \eqref{eq:ql-dual-p-beta}. Combining the above, we have
\[ \eta^t \leq \frac{\sqrt{ 2\left( \varphi(b^t) + g(p^t, \beta^t) \right) }}{\ubar{p}_{\min}}. \]
Note that the above holds for $b^t$ from either PG or PR. Although neat in theory, numerical experiments suggest that the above bound can be loose and is not suitable as a termination criteria.


\section{Leontief utilities}
\subsection{Derivation of \eqref{eq:eg-Leontief-dual-reform}}\label{app:leontief-dual}
The primal EG \eqref{eq:eisenberg-gale-primal} under Leontief utilities $u_i(x_i) = \min_{j\in J_i} \frac{x_{ij}}{a_{ij}}$ can be written in both $x$ and $u$:
\begin{align*}
\min_{u,\, x}\ &  -\sum_i B_i \log u_i \\
{\rm s.t.} \ \ &  u_i \leq \frac{x_{ij}}{a_{ij}} ,\, \forall j\in J_i,\ \forall i\in [n],  \\
& \sum_i x_{ij} \leq 1,\ \forall\, j\in [m],\\
& x\geq 0,\ u\geq 0.
\end{align*}
Clearly, it can also be written in terms of $u_i$ only as follows:
\begin{align}
\min\, -\sum_i B_i \log u_i\ \ {\rm s.t.}\ \sum_{i \in I_j} a_{ij} u_i \leq 1,\ \forall\, j,\ u\geq 0. \label{eq:eg-Leontief}
\end{align}
Let $p_j\geq 0$ be the dual variable associated with constraint $ \sum_{i\in I_j} a_{ij} u_i \leq 1 $. The Lagrangian is 
\begin{align*}
\mathcal{L}(u, p) & = -\sum_i B_i \log u_i + \sum_j p_j \left( \sum_{i\in I_j} a_{ij} u_i - 1 \right) \\
&= - \sum_j p_j + \sum_i \left[ -B_i \log u_i + \langle a_i, p \rangle u_i \right].
\end{align*}
Note that minimizing $\mathcal{L}$ w.r.t. $u$ can be performed separably for each $u_i$. For any $i$ such that $\sum_{j\in J_i} p_j>0$, by first-order stationarity condition, the term $-B_i \log u_i + \langle a_i, p \rangle u_i $ is minimized at $u^*_i(p) = \frac{B_i}{\langle a_i, p \rangle}$
with minimum value $B_i(1 - \log B_i) + B_i \log \langle a_i, p \rangle $.
If $\sum_{j\in J_i} p_j = 0$, the term approaches $-\infty$ as $u_i \rightarrow \infty$. Therefore, the dual objective is
\[ g(p) = \begin{cases}
-\sum_j p_j + 
\sum_i B_i \log \langle a_i, p \rangle + \sum_i B_i (1 - \log B_i) & {\rm if}\ p\geq 0\ {\rm and}\ \sum_{j\in J_i} a_{ij} p_j > 0 \\
-\infty & {\rm o.w.}
\end{cases} \]
Hence the (Lagrangian) dual problem is $\max_p g(p)$. Its minimization form, up to the constant $-\sum_i B_i (1 - \log B_i)$, is
\begin{align}
\min \left[\sum_j p_j - \sum_i B_i \log \langle a_i, p \rangle \right] \ \ {\rm s.t.}\ p \geq 0. \label{eq:eg-Leontief-dual}
\end{align}

By Theorem \ref{thm:eg-gives-me-for-certain-ui}, we have the following.
\begin{itemize}
	\item  An optimal solution to \eqref{eq:eg-Leontief-dual} gives equilibrium prices.
	\item A market equilibrium $(x^*, p^*)$ satisfies $\langle p^*, x_i\rangle = B_i$ for all $i$ and $\sum_i x^*_{ij} = 1$ for all $j$. Therefore, we have $\sum_j p^*_j = \|B\|_1$.
\end{itemize}
Therefore, we can add the constraint $\sum_j p_j = \|B\|_1$ to \eqref{eq:eg-Leontief-dual} without affecting any optimal (equilibrium) solution. This leads to \eqref{eq:eg-Leontief-dual-reform}.

\subsection{Proof of Lemma \ref{lemma:r-ubar-r-bar}} \label{app:leontief-bounds-<a(i),p>}
let $p$ be any feasible solution to \eqref{eq:eg-Leontief-dual-reform}. Since $\sum_j p_j = \|B\|_1$, we have $\langle a_i, p \rangle \leq \|a_i\|_\infty \|p\|_1 = \|a_i\|_\infty \|B\|_1 $ for all $i$. By Appendix \ref{app:leontief-dual}, at equilibrium, $p^*$ and primal variables $u^*_i$ satisfy $u^*_i = \frac{B_i}{\langle a_i, p^* \rangle}$ (the stationarity condition) and 
\[u^*_i \leq \textnormal{utility of getting one unit of every item} = \min_{j\in J_i} \frac{1}{a_{ij}} = \frac{1}{\|a_i\|_\infty}\] 
for all $i$. Therefore $\langle a_i, p^* \rangle = \frac{B_i}{u^*_i} \leq \|a_i\|_\infty \|B\|_1$.

\subsection{Linear convergence of utilities} \label{app:leontief-utilities-lin-conv}
Note that the equilibrium utilities $u^*$ are clearly unique by \eqref{eq:eg-Leontief}. By the KKT stationary condition,  
\[ u^*_i = \frac{B_i}{\langle a_i, p^*\rangle},\ \ \forall\, i \]
for equilibrium prices $p^*$. Therefore, an intuitive construction of $u^t$ is as follows. Let $p^t$ be the current iterate, $r_i^t = \langle a_i, p\rangle$. First compute $\tilde{u}^t_i = \frac{B_i}{r^t_i}$. Then, to satisfy the primal constraints $\sum_i u_i a_{ij} \leq 1$, take 
\[ u^t = \frac{\tilde{u}^t}{\max_j \sum_i u_i a_{ij}  } = \frac{\tilde{u}^t}{\|a^\top \tilde{u}\|_\infty}. \]
Let $r^* = \langle a_i, p^*\rangle = \frac{B_i}{u^*_i}$ and $f^* = \argmin_{p\in \cP} \tilde{h}(ap) = \tilde{h}(r^*) = h(r^*)$. Strong convexity of $\tilde{h}$ implies $\frac{\mu}{2}\|r^t - r^*\|^2\leq h(r^t) - f^*$.
Furthermore, the mapping $r^t\mapsto \tilde{u}^t \mapsto u^t$ is Lipschitz continuous on $r^t \in [\ubar{r}, \bar{r}]$. Therefore, $\|u^t - u^*\|$ converges to $0$ linearly as well. 

\section{Additional details on numerical experiments} \label{app:details-experiment}

For linear utilities, we generate market data $v = (v_{ij})$ where $v_{ij}$ are i.i.d. from standard Gaussian, uniform, exponential, or lognormal distribution. For each of the sizes $n=50, 100, 150, 200$ (on the horizontal axis) and $m=2n$, we generate $30$ instances with unit budgets $B_i=1$ and random budgets $B_i = 0.5 + \tilde{B}_i$ (where $\tilde{B}_i$ follows the same distribution as $v_{ij}$). See \S \ref{sec:experiments} for plots under random budgets and below for those under uniform budgets.

The termination conditions (on the vertical axis) are
\[\epsilon(p^t, p^*) \leq \eta,\ \eta = 10^{-2}, 10^{-3},\] 
where $p^*$ is the optimal Lagrange multipliers of \eqref{eq:eisenberg-gale-primal} computed by CVXPY+Mosek. Then, for $n = 100, 200, 300, 400$ and $n=2m$, we repeat the above with termination conditions
\[{\rm dgap}_t/n \leq \eta,\ \eta = 10^{-3}, 10^{-4}, 10^{-5}, 5\times 10^{-6}.\] 

For QL utilities, we repeat the above (same random $v$, same sizes and termination conditions) using budgets $B_i = 5( 1+\tilde{B}_i)$. This is to make buyers have nonzero bids and leftovers (i.e., $0< \delta^*_i < B_i$) at equilibrium in most scenarios. In this case, $p^* = p(b^*)$, where $b^*$ is the optimal solution to \eqref{eq:ql-shmyrev} computed by CVXPY+Mosek. For QL, FW does not perform well in initial trials and is excluded in subsequent experiments.

For the linesearch subroutine $\mathcal{LS}_{\alpha,\beta,\Gamma}$ in PG (see Appendix \ref{app:pg-ls-details}), we use parameters $\alpha=1.02$, $\beta = 0.8$ and $\Gamma = 100 L\|A\|^2$ throughout.

For Leontief utilities, in addition to ${\rm dgap}_t/n\leq \eta$, we also use the termination condition $\epsilon(u^t, u^*) = \max_j \frac{|u^t_j - u^*_j|}{u^*_j} \leq \eta$, where $u^*$ is the optimal solution to EG under Leontief utilities \eqref{eq:eg-Leontief} computed by CVXPY+Mosek.

\paragraph{Computing the duality gap.} For linear utilities, the objective of the original Shmyrev's convex program \eqref{eq:shmyrev} is
\[ \varphi(b) = -\sum_{i,j} (\log v_{ij})b_{ij} +  \sum_j p_j(b)\log p_j(b) \]
where $p_j(b) = \sum_i b_{ij}$. Recall the objective of the (EG) dual \eqref{eq:eg-dual-p-beta}, equivalent to the dual of Shmyrev's \eqref{eq:shmyrev},
\[g(p, \beta) = \sum_j p_j  - \sum_i B_i \log \beta_i. \]
Given iterate $b^t$, let $p^t_j = p_j(b^t)$ and $\beta^t_i = \min_j \frac{p_j}{v_{ij}}$, which is finite since $v$ is nondegenerate and $p^t>0$. The duality gap is computed via
\[ {\rm dgap}_t = \varphi(b^t) + g(p^t, \beta^t). \] 
For QL utilities, it is computed similarly, that is, through \eqref{eq:ql-shmyrev-dgap-expre}. For Leontief utilities, it is computed using the construction in Appendix \ref{app:leontief-utilities-lin-conv}.


\paragraph{Additional plots.} In \S \ref{sec:experiments}, the plots for linear utilities are generated under random $B_i$. Here we present an augmented set of plots under different random market data, different utilities, unit v.s. varying, random budgets $B_i$ and different termination conditions (${\rm dgap}_t/n \leq \eta$ or $\epsilon(p^t, p^*)\leq \eta$). All horizontal axes are market sizes $n$ (with $m=2n$) and all vertical axes are number of iterations. The error vertical bars are the standard errors of the number of iterations across $30$ repeats. The legends are in the subplot ``Linear utilities, dgap$/n\leq$1e-5''.
\begin{center}
	
	\includegraphics[width=0.8\linewidth]{plots/iters_to_p_diff_vs_mosek/linear-budgets-all-1/iters-to-p_diff-1e-02.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_p_diff_vs_mosek/linear-budgets-all-1/iters-to-p_diff-1e-03.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/linear-budgets-all-1/iters-to-normalized-dgap-1e-03.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/linear-budgets-all-1/iters-to-normalized-dgap-1e-04.pdf} 
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/linear-budgets-all-1/iters-to-normalized-dgap-1e-05.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/linear-budgets-all-1/iters-to-normalized-dgap-5e-06.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/linear-varying-budgets/iters-to-normalized-dgap-1e-03.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/linear-varying-budgets/iters-to-normalized-dgap-1e-04.pdf} 
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/linear-varying-budgets/iters-to-normalized-dgap-1e-05.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/linear-varying-budgets/iters-to-normalized-dgap-5e-06.pdf}
	
	\includegraphics[width=0.8\linewidth]{plots/iters_to_p_diff_vs_mosek/ql-varying-budgets-pg-pr/iters-to-p_diff-1e-02.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_p_diff_vs_mosek/ql-varying-budgets-pg-pr/iters-to-p_diff-1e-03.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/ql-varying-budgets-pg-pr/iters-to-normalized-dgap-1e-03.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/ql-varying-budgets-pg-pr/iters-to-normalized-dgap-1e-04.pdf} 
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/ql-varying-budgets-pg-pr/iters-to-normalized-dgap-1e-05.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/ql-varying-budgets-pg-pr/iters-to-normalized-dgap-5e-06.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/ql-varying-budgets-pg-pr/iters-to-normalized-dgap-1e-03.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/ql-varying-budgets-pg-pr/iters-to-normalized-dgap-1e-04.pdf} 
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/ql-varying-budgets-pg-pr/iters-to-normalized-dgap-1e-05.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/ql-varying-budgets-pg-pr/iters-to-normalized-dgap-5e-06.pdf}
	


	\includegraphics[width=0.8\linewidth]{plots/iters_to_p_diff_vs_mosek/leontief-varying-budgets/iters-to-normalized-dgap-1e-02.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_p_diff_vs_mosek/leontief-varying-budgets/iters-to-normalized-dgap-1e-03.pdf}


	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/leontief-varying-budgets/iters-to-normalized-dgap-1e-03.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/leontief-varying-budgets/iters-to-normalized-dgap-1e-04.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/leontief-varying-budgets/iters-to-normalized-dgap-1e-05.pdf}
	\includegraphics[width=0.8\linewidth]{plots/iters_to_ave_dgap/leontief-varying-budgets/iters-to-normalized-dgap-5e-06.pdf}
\end{center}
{"config": "arxiv", "file": "2006.06747/text/appendix.tex"}