$\displaystyle\stackrel{{\scriptstyle\eqref{eqn:gap<max+gap}}}{{\leq}}\max\\{\mathrm{Gap}(x_{t},\mathcal{C}_{t}),\gamma\mathrm{Gap}(x_{t},\mathcal{C}^{*})\\}+\frac{2L_{\Theta}}{q^{2}}\mathrm{Gap}^{2}(x_{t},\mathcal{C}^{*})+\frac{\gamma{L_{\Theta}}}{q^{2}}\mathrm{Gap}^{2}(x_{t},\mathcal{C}^{*}).$

The above recurrence inequality implies that

$\displaystyle\max\\{\mathrm{Gap}(x_{m+1},\mathcal{C}_{m+1}),\gamma\mathrm{Gap}(x_{m+1},\mathcal{C}^{*})\\}$
(3.51)
$\displaystyle\leq\max\\{\mathrm{Gap}(x_{1},\mathcal{C}_{1}),\gamma\mathrm{Gap}(x_{1},\mathcal{C}^{*})\\}+\sum^{m}_{t=1}\frac{(2+\gamma)L_{\Theta}}{q^{2}}\mathrm{Gap}^{2}(x_{t},\mathcal{C}^{*})$
$\displaystyle\leq\gamma\mathrm{Gap}(x_{1},\mathcal{C}^{*})+\sum^{m}_{t=1}\frac{(2+\gamma)L_{\Theta}}{q^{2}}\mathrm{Gap}^{2}(x_{t},\mathcal{C}^{*}),$

where we use $\mathrm{Gap}(x_{1},\mathcal{C}_{1})=0$ due to that
$\mathcal{C}_{t}$ has only one element. Note that by the definition of index
$m$, we have

$\displaystyle\max\\{\mathrm{Gap}(x_{m+1},\mathcal{C}_{m+1}),\gamma\mathrm{Gap}(x_{m+1},\mathcal{C}^{*})\\}$
(3.52)
$\displaystyle=\max\\{\mathrm{Gap}(x_{m+1},\mathcal{C}^{*}),\gamma\mathrm{Gap}(x_{m+1},\mathcal{C}^{*})\\}$
$\displaystyle=\mathrm{Gap}(x_{m+1},\mathcal{C}^{*}),$

where in the first equality we use the fact
$\mathrm{Gap}(x_{m+1},\mathcal{C}^{*})=\mathrm{Gap}(x_{m+1},\\{\pi_{m+1},\theta_{m+1}\\})$
and
$\\{\pi_{m+1},\theta_{m+1}\\}\subseteq\\{\pi_{m+1}\\}\cup\mathcal{A}_{m}=\mathcal{C}_{m+1}$
due to the definition of $m$. Substituting (3.52) into (3.51) leads to the
following inequality

$\displaystyle\mathrm{Gap}(x_{m+1},\mathcal{C}^{*})$
$\displaystyle\leq\gamma\mathrm{Gap}(x_{1},\mathcal{C}^{*})+\sum^{m}_{t=1}\frac{(2+\gamma)L_{\Theta}}{q^{2}}\mathrm{Gap}^{2}(x_{t},\mathcal{C}^{*})$
$\displaystyle\overset{\eqref{eqn:gap<(1+L)^kgap}}{\leq}\gamma\mathrm{Gap}(x_{1},\mathcal{C}^{*})+\frac{m(2+\gamma)L_{\Theta}}{q^{2}}\left(1+\frac{L_{\Theta}\delta}{q^{2}}\right)^{2m}\mathrm{Gap}^{2}(x_{1},\mathcal{C}^{*})$
$\displaystyle\leq\gamma\mathrm{Gap}(x_{1},\mathcal{C}^{*})+\frac{(2+\gamma)L_{\Theta}|\mathcal{C}^{*}|}{q^{2}}\left(1+\frac{L_{\Theta}\delta}{q^{2}}\right)^{2|\mathcal{C}^{*}|}\mathrm{Gap}^{2}(x_{1},\mathcal{C}^{*})$
$\displaystyle\leq\gamma\mathrm{Gap}(x_{1},\mathcal{C}^{*})+\frac{(2+\gamma)L_{\Theta}|\mathcal{C}^{*}|\delta}{q^{2}}\left(1+\frac{L_{\Theta}\delta}{q^{2}}\right)^{2|\mathcal{C}^{*}|}\mathrm{Gap}(x_{1},\mathcal{C}^{*})$
$\displaystyle\leq 2\gamma\mathrm{Gap}(x_{1},\mathcal{C}^{*}),$

where the last inequality follows if

$\frac{L_{\Theta}|\mathcal{C}^{*}|\delta}{q^{2}}\left(1+\frac{L_{\Theta}\delta}{q^{2}}\right)^{2|\mathcal{C}^{*}|}\leq\frac{\gamma}{2+\gamma},$

which certainly holds if the neighborhood $\mathbb{B}(x^{*},r)$ (i.e., the
$r$) is small enough.

Part (b) Let $\alpha_{t}$ and $d_{t}$ be the step size and joint gradient
generated inside Validation when it is called at the iteration $t$ of
GapReduction, and hence $x_{t+1}=x_{t}-\alpha_{t}d_{t}$, where

$\alpha_{t}=\frac{\left(1-\frac{\gamma}{2}\right)\mathrm{Gap}(x_{t},\mathcal{C}^{*})}{\left\lVert
d_{t}\right\rVert^{2}-\left\langle\nabla
f_{\theta_{t}}(x_{t}),d_{t}\right\rangle}\geq\frac{\left(1-\frac{\gamma}{2}\right)\mathrm{Gap}(x_{t},\mathcal{C}^{*})}{\left\lVert
d_{t}\right\rVert^{2}+G_{f}\left\lVert d_{t}\right\rVert},$ (3.53)

where $G_{f}$ is defined in Table 1. Therefore, moving from $x_{t}$ to
$x_{t+1}$, the function value has decreased by

$\displaystyle
f_{\pi_{t}}(x_{t})-f_{\pi_{t+1}}(x_{t+1})=f_{\pi_{t}}(x_{t})-f_{\pi_{t}}(x_{t+1})+f_{\pi_{t}}(x_{t+1})-f_{\pi_{t+1}}(x_{t+1})$
(3.54) $\displaystyle\geq\alpha_{t}\left\lVert
d_{t}\right\rVert^{2}-\frac{1}{2}\alpha^{2}_{t}L_{\Theta}\left\lVert
d_{t}\right\rVert^{2}+f_{\pi_{t}}(x_{t+1})-f_{\pi_{t+1}}(x_{t+1})$
$\displaystyle\geq\alpha_{t}\left\lVert
d_{t}\right\rVert^{2}-\frac{1}{2}\alpha^{2}_{t}L_{\Theta}\left\lVert
d_{t}\right\rVert^{2}+f_{\pi_{t}}(x_{t}-\alpha_{t}{d_{t}})-f_{\pi_{t+1}}(x_{t}-\alpha_{t}{d_{t}})$
$\displaystyle\geq\alpha_{t}\left\lVert
d_{t}\right\rVert^{2}-\frac{1}{2}\alpha^{2}_{t}L_{\Theta}\left\lVert
d_{t}\right\rVert^{2}+f_{\pi_{t}}(x_{t})-f_{\pi_{t+1}}(x_{t})-\alpha^{2}_{t}L_{\Theta}\left\lVert
d_{t}\right\rVert^{2}$ $\displaystyle\geq\alpha_{t}\left\lVert
d_{t}\right\rVert^{2}-\frac{3}{2}\alpha^{2}_{t}L_{\Theta}\left\lVert
d_{t}\right\rVert^{2},$

where we use the fact
$\pi_{t},\pi_{t+1}\in\mathcal{C}^{*}\setminus\\{\theta_{t}\\}$ and Proposition
1 in deriving the above inequalities. Substituting (3.53) into (3.54) gives

$\displaystyle f_{\pi_{t}}(x_{t})-f_{\pi_{t+1}}(x_{t+1})$
$\displaystyle\geq\frac{\left(1-\frac{\gamma}{2}\right)\mathrm{Gap}(x_{t},\mathcal{C}^{*})}{1+\frac{G_{f}}{\left\lVert
d_{t}\right\rVert}}-\frac{3L_{\Theta}\left(1-\frac{\gamma}{2}\right)^{2}\mathrm{Gap}^{2}(x_{t},\mathcal{C}^{*})}{2\left\lVert
d_{t}\right\rVert^{2}}$
$\displaystyle\geq\frac{\left(1-\gamma\right)\mathrm{Gap}(x_{t},\mathcal{C}^{*})}{1+\frac{G_{f}}{q}}-\frac{3L_{\Theta}\left(1-\frac{\gamma}{2}\right)^{2}\mathrm{Gap}^{2}(x_{t},\mathcal{C}^{*})}{2q^{2}}$
$\displaystyle\geq\frac{\left(1-\gamma\right)\mathrm{Gap}(x_{t},\mathcal{C}^{*})}{2\left(1+\frac{G_{f}}{q}\right)},$

where we use the fact that $\mathrm{Gap}(x_{t},\mathcal{C}^{*})\ll{1}$ in the
last inequality. In particular, for $t=1$ we can obtain (3.37). ∎

######  Theorem 3 (linear rate of convergence).

Suppose Assumption 2 hold. Suppose $x^{*}$ is a non-degenerate local minimizer
of $f$. Suppose $f$ satisfies the multi-component PL inequality in
$\mathbb{B}(x^{*},r)$. Let $\\{x_{n}\\}$ be the points generated by the main
part of Algorithm 2. If $x_{m}$ is sufficiently close to $x^{*}$ for some
index $m$ and if $\mathcal{S}_{n}=\mathcal{C}^{*}$ for every $n\geq{m}$, where
$\mathcal{C}^{*}$ is the complete basis at $x^{*}$, then the following
inequality holds for some constants $G_{f}$ and $q$ that only depend on the
behavior of function $f$ in the neighborhood of $x^{*}$:

$f_{m+k}-f^{*}\leq\left(1-\min\left\\{\frac{\mu_{\mathrm{dec}}\rho_{0}(1-\rho_{0})}{3\mu_{\mathrm{PL}}{L_{\Theta}}},\frac{1-\gamma}{4\left(1+\frac{G_{f}}{q}\right)\mu_{\mathrm{PL}}}\right\\}\right)^{k}(f_{m}-f^{*}),$
(3.55)

where $\mu_{\mathrm{dec}}$, $\rho_{0}$, $L_{\Theta}$, $\mu_{\mathrm{PL}}$ and
$\gamma$ are parameters defined in Table 1.

###### Proof.

Claim 1. For every $n\geq{m}$, if the line search is executed (i.e., the
condition at Line 6 holds) the succeed step-size $\alpha^{*}$ has the
following lower bound

$\displaystyle\alpha^{*}\geq\frac{2\mu_{\mathrm{dec}}(1-\rho_{0})}{3L}.$
(3.56)

Proof of Claim 1. Using Proposition 3, for every $\theta\in\mathcal{C}^{*}$
($\mathcal{S}_{n}=\mathcal{C}^{*}$) the following inequality holds for any
$\alpha>0$:

$\alpha\left\lVert
d_{n}\right\rVert^{2}-\frac{1}{2}L_{\Theta}\alpha^{2}\left\lVert
d_{n}\right\rVert^{2}\leq{f}_{\theta}(x_{n})-f_{\theta}(x_{n}-\alpha
d_{n})\leq\alpha\left\lVert
d_{n}\right\rVert^{2}+\frac{1}{2}L_{\Theta}\alpha^{2}\left\lVert
d_{n}\right\rVert^{2},$ (3.57)

where
$d_{n}=\mathrm{vec(}\mathrm{QP}\textrm{-}{x_{n}}\textrm{-}{\mathcal{S}_{n}}\mathrm{)}=\mathrm{vec(}\mathrm{QP}\textrm{-}{x_{n}}\textrm{-}{\mathcal{C}^{*}}\mathrm{)}$.
Suppose $\alpha^{*}$ is the succeeded step size and
$x_{n+1}=x_{n}-\alpha^{*}d_{n}$ by the updating rule. Suppose
$f_{\theta}(x_{n})=f(x_{n})$ and
$f_{\theta^{\prime}}(x_{n}-\alpha^{*}{d}_{n})=f(x_{n}-\alpha^{*}{d}_{n})$ for
some component indices $\theta$ and $\theta^{\prime}$. Note that we can assume
$\theta^{\prime}\in\mathcal{C}^{*}$ as $x_{n}$ is sufficiently close to
$x^{*}$. Then

$\displaystyle
f(x_{n})-f(x_{n}-\alpha^{*}{d_{n}})=f_{\theta}(x_{n})-f_{\theta^{\prime}}(x_{n}-\alpha^{*}{d_{n}})$
(3.58)
$\displaystyle={f}_{\theta}(x_{n})-f_{\theta}(x_{n}-\alpha^{*}{d_{n}})+f_{\theta}(x_{n}-\alpha^{*}{d_{n}})-f_{\theta^{\prime}}(x_{n}-\alpha^{*}{d_{n}}).$

Note that

$\displaystyle{f}_{\theta}(x_{n})-f_{\theta}(x_{n}-\alpha^{*}{d_{n}})\geq\alpha^{*}\left\lVert
d_{n}\right\rVert^{2}-\frac{1}{2}L_{\Theta}\alpha^{*2}\left\lVert
d_{n}\right\rVert^{2},$ (3.59) $\displaystyle
f_{\theta}(x_{n}-\alpha^{*}{d_{n}})-f_{\theta^{\prime}}(x_{n}-\alpha^{*}{d_{n}})\geq{f}_{\theta}(x_{n})-f_{\theta^{\prime}}(x_{n})-L_{\Theta}\alpha^{*2}\left\lVert
d_{n}\right\rVert^{2}\geq-L_{\Theta}\alpha^{*2}\left\lVert
d_{n}\right\rVert^{2}.$

Substituting (3.59) into (3.58) yields

$f(x_{n})-f(x_{n}-\alpha^{*}{d_{n}})\geq\alpha^{*}\left\lVert
d_{n}\right\rVert^{2}-\frac{3}{2}L_{\Theta}\alpha^{*2}\left\lVert
d_{n}\right\rVert^{2}=\left(1-\frac{3}{2}L_{\Theta}\alpha^{*}\right)\alpha^{*}\left\lVert
d_{n}\right\rVert^{2}.$ (3.60)

The reduction ratio has the lower bound

$\rho=\frac{f(x_{n})-f(x_{n}-\alpha^{*}{d_{n}})}{\alpha^{*}\left\lVert
d_{n}\right\rVert^{2}}\geq 1-\frac{3}{2}L_{\Theta}\alpha^{*}.$

It implies that a sufficient condition for $\rho\geq\rho_{0}$ is
$\alpha^{*}\leq\frac{2(1-\rho_{0})}{3L_{\Theta}}$. Suppose $\alpha^{\prime}$
is the last trial step size that gives $\rho<\rho_{0}$, which implies that
$\alpha^{\prime}>\frac{2(1-\rho_{0})}{3L_{\Theta}}$. Then

$\alpha^{*}=\mu_{\mathrm{dec}}\alpha^{\prime}\geq\frac{2\mu_{\mathrm{dec}}(1-\rho_{0})}{3L_{\Theta}}.$

This concludes the proof of the claim.

Note that since $m$ is sufficiently large, $\mathcal{S}_{n}=\mathcal{C}^{*}$
for every $n\geq{m}$. We analyze separately for two cases.
Case 1. $\left\lVert
d_{n}\right\rVert^{2}\geq\mathrm{Gap}(x_{n},\mathcal{S}_{n})$. In this case we
have

$\displaystyle\frac{f_{n}-f_{n+1}}{f_{n}-f^{*}}\geq\frac{\rho_{0}\alpha^{*}\left\lVert
d_{n}\right\rVert^{2}}{\mu_{\mathrm{PL}}(\left\lVert
d_{n}\right\rVert^{2}+\mathrm{Gap}_{n})}\geq\frac{\rho_{0}\alpha^{*}\left\lVert
d_{n}\right\rVert^{2}}{2\mu_{\mathrm{PL}}\left\lVert
d_{n}\right\rVert^{2}}\geq\frac{\rho_{0}\alpha^{*}}{2\mu_{\mathrm{PL}}}\overset{\eqref{eqn:alpha_succ}}{\geq}\frac{\mu_{\mathrm{dec}}\rho_{0}(1-\rho_{0})}{3\mu_{\mathrm{PL}}{L_{\Theta}}},$
(3.61)

where $f_{n}:=f(x_{n})$ and
$\mathrm{Gap}_{n}:=\mathrm{Gap}(x_{n},\mathcal{S}_{n})$ in the above
inequalities. We have applied the multi-component PL inequality and the step-
size selection rule to derive the first inequality of (3.61), and we have used
the claim to derive the last inequality of (3.61).
Case 2. $\left\lVert
d_{n}\right\rVert^{2}<\mathrm{Gap}(x_{n},\mathcal{S}_{n})$. Using (3.37) from
Lemma 8, in this case we have

$\displaystyle\frac{f_{n}-f_{n+1}}{f_{n}-f^{*}}$
$\displaystyle\geq\frac{(1-\gamma)\mathrm{Gap}_{n}}{2\left(1+\frac{G_{f}}{q}\right)}\cdot\frac{1}{\mu_{\mathrm{PL}}(\left\lVert
d_{n}\right\rVert^{2}+\mathrm{Gap}_{n})}$ (3.62)
$\displaystyle\geq\frac{(1-\gamma)\mathrm{Gap}_{n}}{2\left(1+\frac{G_{f}}{q}\right)}\cdot\frac{1}{2\mu_{\mathrm{PL}}\mathrm{Gap}_{n}}$
$\displaystyle\geq\frac{1-\gamma}{4\left(1+\frac{G_{f}}{q}\right)\mu_{\mathrm{PL}}}.$

Combining the two cases leads to the linear rate of convergence (3.55). ∎

###  3.6 Global rate of convergence

Figure 3: Neighborhood regions of the unique local minimizer $x^{*}$ used in
the proof of Theorem 4.

Put all things together, we have the following result on the global rate of
convergence for Algorithm 2.

######  Theorem 4.

Suppose Assumption 2 hold, and $f$ is a locally-max representable. Let $x^{*}$
be the unique local and global minimizer of $f$ and $x^{*}$ is non-degenerate.
Suppose the multi-component Polyak-Łojasiewicz inequality hold at a
neighborhood of of $x^{*}$. Apply Algorithm 2 to minimize $f$, and let
$\\{x_{n}\\}$ be the sequence of points generated from the main iterations of
the algorithm. Then there exists an integer $m>0$ and parameters
$c_{0},\epsilon_{0},\epsilon_{1},q_{0},q_{1}>0$ such that

$\displaystyle
m\lesssim\frac{f^{0}-f^{*}}{c_{0}}\left(\log\frac{w_{\textrm{init}}}{\min\\{\epsilon_{1},q^{2}_{0}\\}}\right)^{2}+\log\frac{w_{\textrm{init}}}{\epsilon_{0}},$
(3.63) $\displaystyle
f_{m+k}-f^{*}\leq\left(1-\min\left\\{\frac{\mu_{\mathrm{dec}}\rho_{0}(1-\rho_{0})}{3\mu_{\mathrm{PL}}{L_{\Theta}}},\frac{1-\gamma}{4\left(1+\frac{G_{f}}{q_{1}}\right)\mu_{\mathrm{PL}}}\right\\}\right)^{k-1}(f_{m}-f^{*}),$
(3.64)

where $f_{n}=f(x_{n})$, $f^{*}=f(x^{*})$ and $f^{0}=f(x_{0})$. The constants
$c_{0},\epsilon_{0},\epsilon_{1},q_{0},q_{1}$ only depend on the properties of
the function $f$ but independent of the process of the algorithm. The
$w_{\textrm{init}}$ is in the input of the algorithm, and the other parameters
$\mu_{\mathrm{dec}}$, $\rho_{0}$, $L_{\Theta}$, $\mu_{\mathrm{PL}}$, $\gamma$
and $G_{f}$ are defined in Table 1.

###### Proof.

Without loss of generality, we focus on the general case that
$\widetilde{\mathcal{M}}(x^{*})\setminus\mathcal{C}^{*}$ is non-empty. The
argument for the case of $\widetilde{\mathcal{M}}(x^{*})=\mathcal{C}^{*}$ is
simpler. Let $U(v):=\\{x:f(x)\leq{v}\\}$ be the level sets of $f$. There exist
$r_{1}$, $r_{2}$, $\epsilon_{0}$, $\epsilon_{1}$, $q_{0}$ and $v$ which only
depend on $f$ making the following properties hold:

$\displaystyle r_{2}>r_{1}>0,\;q_{0},\epsilon_{0},\epsilon_{1}>0,$ (3.65)
$\displaystyle r_{1}\textrm{ satisfies the condition in
Lemma~{}\ref{lem:lb_dec}},$ (3.66) $\displaystyle r_{2}\textrm{ satisfies the
conditions in Lemma~{}\ref{lem:validation}, Lemma~{}\ref{lem:gap_reduction}
and Theorem~{}\ref{thm:linear_rate}},$ (3.67)
$\displaystyle\mathbb{B}(x^{*},r_{1})\subseteq{U}(v)\subseteq\mathbb{B}(x^{*},r_{2}),$
(3.68)
$\displaystyle\widetilde{\mathcal{M}}(x)=\widetilde{\mathcal{M}}(x^{*})\qquad\forall
x\in\mathbb{B}(x^{*},r_{2}),$ (3.69)
$\displaystyle\left\lVert\mathrm{vec(}\mathrm{QP}\textrm{-}{x}\textrm{-}{\mathcal{C}^{*}}\mathrm{)}\right\rVert\leq\epsilon_{0}\qquad\forall
x\in\mathbb{B}(x^{*},r_{2}),$ (3.70)
$\displaystyle\mathrm{Gap}(x,\mathcal{C}^{*})\leq\epsilon_{0}\qquad\forall
x\in\mathbb{B}(x^{*},r_{2}),$ (3.71) $\displaystyle
4\epsilon_{0}\leq{f}(x)-\max_{\theta\in\mathcal{C}^{\prime}}f_{\theta}(x)\leq{6\epsilon_{0}},\;\mathcal{C}^{\prime}=\widetilde{\mathcal{M}}(x^{*})\setminus\mathcal{C}^{*}\qquad\forall
x\in\mathbb{B}(x^{*},r_{2}),$ (3.72)
$\displaystyle\left\lVert\mathrm{vec(}\mathrm{QP}\textrm{-}{x}\textrm{-}{\mathcal{C}}\mathrm{)}\right\rVert\geq{q_{0}}\qquad\forall
x\in{X}\setminus\mathbb{B}(x^{*},r_{1})\;\forall\mathcal{C}\subseteq\mathcal{S}(x,\epsilon_{1}),$
(3.73) $\displaystyle
w_{k}\geq\frac{1}{2}\min\\{\epsilon_{1},q^{2}_{0}\\}\qquad\forall
x_{k}\in{X}\setminus\mathbb{B}(x^{*},r_{1}),$ (3.74)

where the achievability of (3.68) and (3.69) is ensured by Observation 2 and
Assumption 2, respectively. The achievability of (3.70), (3.71) and (3.72) is
ensured by the continuity of the component functions and the fact that
$f(x^{*})-\max_{\theta\in\mathcal{C}^{\prime}}f_{\theta}(x^{*})>0$. The
achievability of (3.73) and (3.74) is ensured by Proposition 5 and Proposition
7, respectively.

Claim 1. $x_{n}\in{U(v)}\subseteq\mathbb{B}(x^{*},r_{2})$ for any
$n\geq{n_{0}}$, where $n_{0}$ is an integer having the following upper bound:

$n_{0}\lesssim\frac{f^{0}-f^{*}}{c_{0}}\left(\log\frac{w_{\textrm{init}}}{\min\\{\epsilon_{1},q^{2}_{0}\\}}\right)^{2}.$
(3.75)

Proof. Let $k_{1}$ and $k_{2}$ be two consecutive iteration indices satisfying
$x_{k_{1}},x_{k_{2}}\in{X}\setminus{U(v)}$, $\left\lVert
d_{k_{1}}\right\rVert^{2}\geq\mathrm{Gap}(x_{k_{1}},\mathcal{S}_{k_{1}})$, and
$\left\lVert
d_{k_{2}}\right\rVert^{2}\geq\mathrm{Gap}(x_{k_{2}},\mathcal{S}_{k_{2}})$. By
hypothesis, we have $\left\lVert
d_{k}\right\rVert^{2}<\mathrm{Gap}(x_{k},\mathcal{S}_{k})$ for every $k$
between $k_{1}$ and $k_{2}$. For every $k_{1}<k<k_{2}$, there are only two
possible cases listed as follows

  1. (1)

$w_{k+1}=w_{k}/2$ by executing Line 16;

  2. (2)

$w_{k+1}=w_{k}$ by executing Line 18 due to that the condition at Line 8 of
GapReduction is violated (returning $status=1$). The later implies that
$\mathcal{S}_{k+1}=\mathcal{S}_{k}$ and
$\mathrm{Gap}(x_{k+1},\mathcal{S}_{k+1})\leq\mathrm{Gap}(x_{k},\mathcal{S}_{k})/2$.

Suppose $m_{1}<m_{2}$ are two consecutive indices between $k_{1}$ and $k_{2}$
such that the case (1) is satisfied at $x_{m_{1}}$ and $x_{m_{2}}$. By
hypothesis, $x_{m}$ satisfies the case (2) at every index $m_{1}<m<m_{2}$.
Furthermore, for every $m_{1}<m<m_{2}$, we should have

$\left\lVert
d_{m}\right\rVert^{2}=\left\lVert\mathrm{vec(}\mathrm{QP}\textrm{-}{x_{m}}\textrm{-}{\mathcal{S}_{m}}\mathrm{)}\right\rVert^{2}<\mathrm{Gap}(x_{m},\mathcal{S}_{m})\leq{w_{m}},$
(3.76)

by hypothesis due to $k_{1}<m_{1}<m<m_{2}<k_{2}$. We argue that the following
bound holds:

$m_{2}-m_{1}\leq{O}\left(\log\frac{w_{\textrm{init}}}{\min\\{\epsilon_{1},q^{2}_{0}\\}}\right).$
(3.77)

To see this, we notice that the GapReduction procedure executed at Line 14
succeeds (i.e., $status=1$) for any $m_{1}<m<m_{2}$ by the definition of
$m_{1}$ and $m_{2}$. Furthermore, violating the condition at Line 8 of
GapReduction indicates

$\mathrm{Gap}(x_{m+1},\mathcal{S}_{m+1})\leq\mathrm{Gap}(x_{m},\mathcal{S}_{m})/2\quad\forall
m_{1}<m<m_{2}-2.$

It follows that

$\mathrm{Gap}(x_{m_{2}-1},\mathcal{S}_{m_{2}-1})\leq\left(\frac{1}{2}\right)^{m_{2}-m_{1}-2}\mathrm{Gap}(x_{m_{1}+1},\mathcal{S}_{m_{1}+1})\leq\left(\frac{1}{2}\right)^{m_{2}-m_{1}-2}w_{\textrm{init}}.$

We now prove by contradiction. Suppose (3.77) does not hold. The above
inequality then yields

$\mathrm{Gap}(x_{m_{2}-1},\mathcal{S}_{m_{2}-1})\leq\frac{1}{2}\min\\{\epsilon_{1},q^{2}_{0}\\},$
(3.78)

where the factor $1/2$ in the above inequality is valid due to the constant
factor in (3.77) is adjustable. Then Observation 1 implies that
$\mathcal{S}_{m^{\prime}}\subseteq\mathcal{S}(x_{m^{\prime}},\epsilon_{1})$.
By (3.73), one has

$\left\lVert d_{m_{2}-1}\right\rVert\geq{q_{0}}.$ (3.79)

However, substituting (3.78) and (3.79) into (3.76) (with $m=m_{2}-1$) leads
to a contradicting inequality $q^{2}_{0}<\min\\{\epsilon_{1},q^{2}_{0}\\}/2$,
which proves the bound (3.77).
Let $I=\\{m_{i}:k_{1}<m_{i}<k_{2},\;\textrm{the case (1) holds at iteration
$m_{i}$}\\}$. The condition (3.74) implies that

$|I|\leq{O}\left(\log\frac{w_{\textrm{init}}}{\min\\{\epsilon_{1},q^{2}_{0}\\}}\right).$
(3.80)

Applying (3.77) and (3.80) gives

$k_{2}-k_{1}\lesssim\left(\log\frac{w_{\textrm{init}}}{\min\\{\epsilon_{1},q^{2}_{0}\\}}\right)^{2}.$
(3.81)

Let $K=\\{k_{i}:x_{k_{i}}\in{X}\setminus{U(v)},\;\left\lVert
d_{k_{i}}\right\rVert^{2}\geq\mathrm{Gap}(x_{k_{i}},\mathcal{S}_{k_{i}})\\}$.
By (3.66), (3.68) and Lemma 2, we have

$f(x_{k_{i}})-f(x_{k_{i}+1})\geq{c_{0}}\qquad\forall{k_{i}}\in{K},$ (3.82)

for some constant $c_{0}>0$ that only depends on $r_{1},\epsilon_{1},q_{0}$,
and intrinsic parameters of $f$. It follows that

$|K|\lesssim\frac{f^{0}-f^{*}}{c_{0}}.$ (3.83)

Combining (3.81) and (3.83) gives (3.75), concluding the proof of Claim 1.

Claim 2. $x_{n}\in{U(v)}$ and $w_{n}<6\epsilon_{0}$ for any
$n\geq{n_{0}+n_{1}}$, where $n_{0}$ is the number given in Claim 1, and
$n_{1}$ is a number satisfies

$n_{1}\lesssim\log\frac{w_{\textrm{init}}}{\epsilon_{0}}.$ (3.84)

Proof. By Claim 1, the sequence of points will hit the $v$-level set $U(v)$
after at most $n_{0}$ iterations, and the sequence will stay in $U(v)$
afterwards due to Observation 2. We argue that $w_{n+1}=w_{n}/2$ for any
iteration $n\geq{n_{0}}$ satisfying $w_{n},w_{n+1}>6\epsilon_{0}$. To see it.
First, notice that for such an iteration $n$, (3.72) implies that
$\mathcal{S}_{n}=\mathcal{S}(x_{n},w_{n})\supseteq\mathcal{S}(x_{n},6\epsilon_{0})\supseteq\mathcal{C}^{*}$.
Therefore,

$\displaystyle\left\lVert
d_{n}\right\rVert=\left\lVert\mathrm{vec(}\mathrm{QP}\textrm{-}{x_{n}}\textrm{-}{\mathcal{S}_{n}}\mathrm{)}\right\rVert\leq\left\lVert\mathrm{vec(}\mathrm{QP}\textrm{-}{x_{n}}\textrm{-}{\mathcal{C}^{*}}\mathrm{)}\right\rVert\overset{\eqref{eqn:global_5}}{\leq}\epsilon_{0},$
$\displaystyle\mathrm{Gap}(x_{n},\mathcal{S}_{n})\overset{\eqref{eqn:global_7}}{\geq}4\epsilon_{0},\;\mathrm{Gap}(x_{n+1},\mathcal{S}_{n+1})\overset{\eqref{eqn:global_7}}{\geq}4\epsilon_{0}$

where the hypothesis $w_{n},w_{n+1}>6\epsilon_{0}$ has been used to get the
second row of inequalities. Therefore, the condition at Line 6 of the
algorithm cannot be satisfied at the iteration $n$. Let
$I^{\prime}=\\{n:n\geq{n_{0}},\;w_{n},w_{n+1}\geq 6\epsilon_{0}\\}$. Similarly
as the argument in the proof of Claim 1, for any $n\in{I}^{\prime}$ we either
have $w_{n+1}=w_{n}/2$ or
$\mathrm{Gap}(x_{n+1},w_{n+1})=\mathrm{Gap}(x_{n},w_{n})/2$. On the other
side, since $6\epsilon_{0}\leq{w}_{n},w_{n+1}\leq{w}_{\textrm{init}}$ and
$4\epsilon_{0}\leq\mathrm{Gap}(x_{n},\mathcal{S}_{n}),\mathrm{Gap}(x_{n+1},\mathcal{S}_{n+1})\leq{w}_{\textrm{init}}$.
It follows that

$|I^{\prime}|\lesssim\log\frac{w_{\textrm{init}}}{\epsilon_{0}},$

Therefore, after at most $O(\log\frac{w_{\textrm{init}}}{\epsilon_{0}})$
iterations starting from $n_{0}$, the gap parameter $w_{n}$ will be smaller
than $6\epsilon_{0}$, which concludes the proof of Claim 2.

Claim 3. Let $J=\\{n:x_{n}\in{U(v)},\;n>{n_{0}+n_{1}}\\}$, where $n_{0}$ and
$n_{1}$ are given in Claim 1 and Claim 2, respectively. Then (a) there is at
most one index $k_{0}\in{J}$ such that
$\mathcal{S}_{k_{0}}\setminus\mathcal{C}^{*}$ is non-empty, and (b)
$\mathcal{S}_{n}\subseteq\mathcal{C}^{*}$ for any
$n\in{J}\setminus\\{k_{0}\\}$. Proof. First, we note that
$w_{n}<6\epsilon_{0}$ and hence $\mathrm{Gap}(x_{n},w_{n})<6\epsilon_{0}$ for
any $n\in{J}$. Let $k_{0}$ be the first index in $J$ making
$\mathcal{S}_{k_{0}}\setminus\mathcal{C}^{*}$ non-empty, if it exists. With
the same argument as in the Claim 2, it can be shown that the condition at
Line 6 of the algorithm cannot be satisfied. We argue that the $status$
returned by GapReduction must be 0. Assuming that $status=1$, then the
condition at Line 8 of GapReduction should not hold, and it implies that
$\mathcal{S}_{k_{0}+1}=\mathcal{S}_{k_{0}}$, $x_{k_{0}+1}\in{U(v)}$ due to
$f(x_{k_{0}+1})\leq{f}(x_{k_{0}})$ and Observation 2, and also

$\frac{\mathrm{Gap}(x_{k_{0}+1},\mathcal{S}_{k_{0}+1})}{\mathrm{Gap}(x_{k_{0}},\mathcal{S}_{k_{0}})}\leq\frac{1}{2}.$
(3.85)

On the other hand, since $\mathcal{S}_{k_{0}+1}=\mathcal{S}_{k_{0}}$ which
contains at least one component from
$\widetilde{\mathcal{M}}(x^{*})\setminus\mathcal{C}^{*}$ by the definition of
$k_{0}$, (3.72) points out that both
$\mathrm{Gap}(x_{k_{0}},\mathcal{S}_{k_{0}})$ and
$\mathrm{Gap}(x_{k_{0}+1},\mathcal{S}_{k_{0}+1})$ are between $4\epsilon_{0}$
and $6\epsilon_{0}$. It contradicts to (3.85). Therefore, we have $status=0$
and hence Line 16 of the algorithm will be executed at the iteration $k_{0}$,
which gives $w_{k_{0}+1}=w_{k_{0}}/2<3\epsilon_{0}$, and hence
$\mathcal{S}_{n}\subseteq\mathcal{C}^{*}$ for any $n\geq{k_{0}+1}$ because of
(3.71) and (3.72). This concludes the proof of Claim 3.

Theorem 3 and Claim 3 show that the algorithm achieves a linear reduction on
$f$ at every iteration $n\in{J}\setminus\\{k_{0}\\}$ leading to (3.64), where
the power index in (3.64) is $k-1$ instead of $k$ due to Claim 3. This
concludes the proof of the theorem. ∎

######  Remark.

The above theorem points out that there exists a iteration index $m$ such that
the average decrement per iteration is lower bounded by a constant $c_{0}$ in
the first $m$ iterations, and it achieves a linear rate of convergence in the
following-up iterations.

######  Corollary 2 (complexity).

Let all the assumptions in Theorem 4 hold. For any sufficiently small
$\epsilon>0$, the computational complexity of Algorithm 2 is
$O(\log\left(\frac{1}{\epsilon}\right))$ for reaching a point $x$ satisfying
$f(x)-f(x^{*})\leq\epsilon$.

###### Proof.

By Claim 2 in the proof of Theorem 4, it takes a bounded number $n_{0}+n_{1}$
of iterations for Algorithm 2 to enter the level set $U(v)$ and have
$w_{n}<6\epsilon_{0}$ for every $n\geq{n_{0}+n_{1}}$. By Claim 3, for
$n\in{J}$ with $J=\\{n:x_{n}\in{U(v)},\;n>{n_{0}+n_{1}}\\}$, there is at most
one iteration in which the algorithm does not achieve a linear rate of
decrement in the function value. Therefore, the number of outer iterations of
while-loop at Line 3 of Algorithm 2 is of
$O(\log\left(\frac{1}{\epsilon}\right))$. It amounts to show that the number
of operations within each outer iteration is a bounded number. Indeed, it
takes a bounded number of inner iterations to identify a qualified step size
$\alpha$ in the while-loop at Line 8 which has been proved in (3.11) and
(3.56). Furthermore, the while-loop at Line 3 of GapReduction takes a bounded
number of iterations to execute since $|\mathcal{C}_{\textrm{in}}|$ is a
bounded number. This concludes the proof. ∎

##  4 Numerical investigation

Practical implementation. The Algorithm 1 and Algorithm 2 are designed to
satisfy the requirement on convergence and the rate of convergence. For
implementation, some procedures in these algorithms can be omitted, while some
empirical rules are needed to improve the speed of convergence in practice. A
template for implementing the joint gradient descent method with component
information is given in Algorithm 3, which provides a guideline for practical
implementation. The algorithm consists of four major steps: component
selection (selection of previously encountered component functions to compute
the joint gradient), joint-gradient computation, line search and component
information updating. The algorithm has been implemented in Julia, and the
solver for quadratic program is Ipopt [38]. The code can be accessed via the
following link:

https://github.com/luofqfrank/CID-JGD.git.

Our numerical experience reveals that component selection is crucial in later
iterations to maintain an effective reduction of the objective value per
iteration, especially when the solution is close to a local minimizer. The
rule of component selection also depends on whether all the component
functions have a uniform domain or not. For example, consider the function
$f(x):=\max\\{f_{1}(x),f_{2}(x)\\}$. It has two component functions $f_{1}$
and $f_{2}$ which have an uniform domain, i.e., at any point $x\in\mathbb{R}$,
the evaluation of $f_{1}$ and $f_{2}$ at $x$ is achievable. On the other hand,
for the function $f(x)=|x|^{x^{2}}$, it has two component functions
$f_{1}(x)=x^{x^{2}}$ defined on $x\in[0,\infty)$ and $f_{2}(x)=(-x)^{x^{2}}$
defined on $x\in(-\infty,0]$, and hence they do not have a uniform domain.

We have tried and implemented the following rules for component selection: 1\.
Select components whose representative points $x_{\theta}$ are within a
distance $r$ from the current point $x$ and $r$ can be an adjustable parameter
over iterations. This rule is applicable to both uniform-domain and non-
uniform-domain problem instances. 2\. Select components whose function values
are within a gap $\delta$ from the active component evaluating at the current
point $x$, i.e., $\mathcal{C}\leftarrow\\{\theta:f_{\theta_{0}}(x)-\delta\leq
f_{\theta}(x)\leq{f}_{\theta_{0}}(x)\\}$, and the gap $\delta$ is adjustable
over iterations. This rule is applicable to uniform-domain problem instances.
3\. For the situation that the function-value reduction becomes very small, we
proposed a more accurate way of managing neighborhood component functions. We
maintain a dynamic set $\mathcal{C}$ in the following manner: When the
algorithm encounters a new blocking component $\theta$, we include it in
$\mathcal{C}$ for the next iteration. Once $|\mathcal{C}|$ exceeds a certain
capacity or the proximity measure (either by the distance or the gap) of
components in $\mathcal{C}$ exceeds a certain threshold, we reset
$\mathcal{C}$ to be a singleton set containing only the current active
component. Setting the capacity for $\mathcal{C}$ intends to limit the CPU
time for solving each QP for the joint-gradient computation especially for
large scale problem instances. The capacity of $\mathcal{C}$ is set to 50 in
our experiments, and our numerical experience shows that for large instances
($n\geq 1000$), solving a single joint-gradient QP with about $100$ components
will take $20\textrm{s}\sim 30\textrm{s}$. This rule is applicable to both
uniform-domain and non-uniform-domain problem instances. Note that the above
rules can be combined to apply at different stages of the iteration course,
and their effectiveness can vary. More advanced rules can be investigated in
following-up studies.

Test problems. We consider 10 standard test problems from [18] with 6
different dimensions: $n=50,100,500,1000,2000\textrm{ and }5000$ in our
experiments. The notation and properties of the test problems are summarized
in Table 2.

Non-smooth optimization methods in the comparison. We compare the performance
of our component-information-driven joint gradient descent method with other
three advanced methods for non-smooth optimization developed in literature.
The methods and their notations are listed as follows:

  * -

JGD: the joint gradient descent method with component information (Algorithm
3);

  * -

TRB: the trust-region bundle method [33];

  * -

QNS: the quasi-Newton method with sampling [7];

  * -

VM: the variable metric method [37].

Termination criteria. Determining termination criteria for non-smooth
optimization problems is in general more difficulty than for smooth
optimization problems, as because it is almost impossible to judge if the
current solution is in a small neighborhood of a local minimizer simply based
on the sub-gradient at the point. Furthermore, since different methods work
very differently, it adds another layer of difficulty to defining a uniform
set of termination criteria for all methods. But the overall guideline should
covers the following three aspects:

  * -

Locality and stationarity: The algorithm gets trapped in a a sufficiently
small region and has collected sufficient information about the sub-gradient
information within that region. An evaluation of the “mean” sub-gradient
within the small region shows that it is lower than a tolerance threshold;

  * -

Decrement of objective in consecutive iterations: When the algorithm has a
consistent small decrement of objective that is lower than a threshold per
iteration, but the stationarity condition is not satisfied, we terminate the
algorithm;

  * -

Time out: the algorithm terminates if the CPU time exceeds a threshold.

The above termination types are referred as Type I, II and III terminations,
respectively. In our experiments, for Type I termination, we set the locality
threshold to be $1\times{10}^{-5}$ and the stationarity tolerance to be
$1\times{10}^{-3}$. For Type II termination, we monitor 10 consecutive
iterations (or consecutive successful iterations for the VM method) and check
if the decrement per iteration is less than $1\times{10}^{-8}$. For Type III
termination, we set the maximum CPU time to be 20 minutes. For the JGD method,
we collect the component functions with representative points within the
locality threshold from the current point and compute the joint gradient to
test if its magnitude is within the stationarity tolerance. This is actually
another advantage of logging the component information: It can be utilized to
evaluate the locality and stationarity condition more efficiently than other
methods. Our numerical experience points out that the quality of the joint-
gradient computation is sensitive to the tolerance error of the QP solver. A
low-quality joint gradient can fail the line search (i.e., the inner loop for
line search does not terminate) in the later stage of the solution process.
Therefore, we set the tolerance of the QP solver to be $1\times{10}^{-16}$.

Performance on median-size problem instances All the methods and numerical
experiments are implemented in Julia and the automatic differentiation
technology has been applied to compute the gradient of a given component
function. The computation is performed on a 2.3 GHz 4-core computer with 32 GB
RAM. Each job of solving an particular problem instance using a method is
managed by a single processor (i.e., no parallel processing). For median-size
problem instances ($n=50,100$), we compare the performance of JGD with other
three methods TRB, QNS and VM. Figure 4 shows the performance curves of the
four methods on solving the 10 test problems for $n=50$ and $n=100$,
respectively. The performance curve of a solver is defined as $\rho_{s}(\tau)$
versus $\tau$, where $\tau$ is a changing parameter and $\rho_{s}(\tau)$ is
defined as

$\rho_{s}(\tau):=\frac{\textrm{no. of problems where
}\log_{2}(r_{p,s})\leq\tau}{\textrm{total no. of problems}},$

where $s$ is the index of the solver, $p$ is the index of problem instance,
$r_{p,s}$ is the ratio between the time to solve problem $p$ by solver $s$
over the lowest time required by any of the solvers. The ratio $r_{p,s}$ is
set to infinity if solver $s$ cannot terminate with the computational time
limit (20 minutes in this case). Table 3 gives the termination times and
objective values achieved by every method on each problem instance, and Figure
5 shows the reduction of relative optimality gap (ROG) over time for all
objective functions with $n=100$. The plot of performance curves in Figure 4
points out that the JGD is the best of all methods considered which can solve
9 out of 10 instances with Type I and Type II termination criteria for both
$n=50$ and $n=100$. From Table 3, the JGD method achieves shorter termination
time relative to the other three methods on all instances with only few
exceptions, and the quality of solution given by the JGD method is in general
magnitudes better than that of the other three methods, which can also be seen
from Figure 5.

Performance on large problem instances We have noticed that the TRB, QNS and
VM method cannot make effectively progress on solving large scale problem
instances due to that they need to solve a $n$-dimensional quadratic program
with matrix size being $n\times{n}$ and update the matrix in each iteration,
which is very computationally expensive for $n\geq 1000$. On the other side,
for the JGD method, the dimension of the joint-gradient QP is equal to the
number of selected components (i.e., the size of $\mathcal{C}$) which is no
greater than 50 as we have set a capacity for $\mathcal{C}$. Furthermore, the
joint-gradient QP has a much lower complexity relative to QP’s in the other
three methods as it does not require a explicit construction of the matrix.
These features make the scale of the joint-gradient QP more robust towards the
growth of dimension, and hence more efficient to be solved. In addition to
these properties, the joint gradient descent is a more informative search
direction by intuition. Therefore, it is not surprising that the JGD has
essential advantages in solving large scale problem instances compared to
other methods. The performance of using the JGD method to solve large problem
instances ($n=500,1000,2000,5000$) is summarized in Table 4 and the relative
optimality gap over iterations is visualized in Figure 6. In particular, Table
4 shows that for 15 out of 40 instances, the JGD method can terminate in less
than 1 minute. For all instances of gen_MXHILB, Chained_CB3_II and
Chained_Crescent_I, the JGD method can solve them within $1\sim 60$ seconds
with the relative optimality gap (ROG) from $1\times 10^{-12}$ to $1\times
10^{-4}$. For all instances of Chained_LQ, Chained_CB3_I, brown_func2,
Chained_Mifflin2 and Chained_Crescent_II, the JGD method can effectively
reduce the ROG to the range of $1\times 10^{-3}\sim 1\times 10^{-6}$. For
gen_MAXQ instances, the JGD can reduce ROG to $1\times 10^{-4}\sim 1\times
10^{-3}$ for $n=500,1000$, but it does not perform well on $n=2000,5000$,
leaving ROG up to $0.4$. For all instances of num_active_faces, the ROG is up
to $0.1\sim 0.4$. Figure 6 shows the efficiency of reducing ROG of 6 selected
objective functions for different dimensions. We notice that for Chained_LQ
and Chained_CB3_I, the efficiency of ROG reduction is relatively consistent
among dimensions. For Chained_CB3_II and Chained_Crescent_I, the ROG can be
reduced down to roughly a same level ($<1\times 10^{-}9$) among dimensions,
but it takes about 5 time amount of computational time for $n=5000$ compared
to $n=500,1000,2000$. For brown_func2, the efficiency very typically depends
on the dimension, as the gaps among these curves are very clear. For
Chained_Crescent_II, the efficiency for $n=1000,2000,5000$ is very close,
while the efficiency for $n=500$ is about 100 times better than the other
three problem sizes. The violin plots in Figure 7 show the number
$|\mathcal{C}|$ of selected components for joint-gradient calculation and its
distribution over iterations. The number $|\mathcal{C}|$ ranges from 1 to 80
as because we have set a capacity for $\mathcal{C}$ and reset $\mathcal{C}$
once its size exceeds the capacity.

##  5 Concluding Remarks

Potential research directions. The paper demonstrates the power of using
component information in designing algorithms for solving practical non-smooth
optimization problems from the theoretical and numerical perspectives.
Different ways of using the component information in algorithm design and how
to integrate this information with empirical rules or machine learning based
rules for component selection to further improve the algorithm performance
becomes an emerging research direction. In addition, following-up methodology
development is needed for various problem settings. For example, suppose some
component functions are in a stochastic form, e.g.,
$f_{\theta}(x)=\mathbb{E}_{\xi}[\psi(x,\xi)]$. Then the evaluation of the
function value and gradient of $f_{\theta}$ becomes stochastic. One can
imagine that the stochastic ingredient in this problem setting can enrich the
algorithm design and the convergence analysis. For another example, suppose
the objective function is in a sample-average form:
$f(x)=\frac{1}{N}\sum^{N}_{i=1}g(x,\xi_{i})$, and $g(x,\xi)$ is an encodable
piecewise-smooth function of $x$. In this case, it is not a good idea to log
the component of each $g(x,\xi_{i})$ if the sample size $N$ is a big number,
as because the complexity of codes grows linearly with $N$, and the total
number of codes grows exponentially with $N$. As a consequence, the importance
of a specific component function is diluted. Therefore, how to extract high-
level component information and justify it in theoretical and numerical
analysis becomes an interesting problem to investigate.

An open problem. In the theorems for linear rate of convergence, i.e., Theorem
3 and Theorem 4, it is assumed that the minimizer is non-degenerate. An open
problem is: Is it possible to design an algorithm that can achieve a linear
rate of convergence for a general case that the minimizer is possibly
degenerate. As an intuition, the analysis for this case can be more intricate
and additional validation procedures are needed in the design of such an
algorithm.

Algorithm 3 A practical implementation for joint gradient descent method with
component information

Parameters: $0<\alpha_{0},\mu_{0},\mu_{\mathrm{dec}}<1$

Data structure: A dictionary dict that maps the code $\theta$ to
$[x_{\theta},f_{\theta}(x_{\theta}),\nabla{f}_{\theta}(x_{\theta})]$

Step 1. Initialization: $x\leftarrow{x_{0}}$. Get the code $\theta_{0}$ of
active component function at $x_{0}$. Set
$\textsc{dict}[\theta_{0}]=[x_{0},f_{\theta_{0}}(x_{0}),\nabla{f}_{\theta_{0}}(x_{0})]$.

Step 2. Component selection: Select a set $\mathcal{C}$ (active set) of codes
from dict according to a certain rule based on the proximity of $x_{\theta}$
with respect to the current point $x$.

Step 3. Computation of the joint gradient: Solve the quadratic program

$\min_{\lambda}\left\lVert\sum_{\theta\in\mathcal{C}}\lambda_{\theta}\nabla{f}_{\theta}(x_{\theta})\right\rVert^{2}\;\textrm{s.t.
}\sum_{\theta\in\mathcal{C}}\lambda_{\theta}=1$ to determine the joint
gradient $d$ and let $\hat{d}=d/\left\lVert d\right\rVert$ be the normalized
vector of $d$.

Step 4. Line search: Starting from $\alpha\leftarrow\alpha_{0}$, reduce
$\alpha\leftarrow\mu_{\mathrm{dec}}\alpha$ until one finds a step size
$\alpha$ such that $[f(x)-f(x-\alpha\hat{d})]/(\alpha\left\lVert
d\right\rVert)\geq\mu_{0}$.

Step 5. Component information update: During Step 4, for each test step size
$\alpha^{\prime}$, determine the active code $\theta$ at
$x-\alpha^{\prime}\hat{d}$. If $\theta\notin\textsc{dict}$ or
$\theta\in\textsc{dict}$ but $\|\alpha^{\prime}\hat{d}\|<\left\lVert
x_{\theta}-x\right\rVert$, set
$x_{\theta}\leftarrow{x}-\alpha^{\prime}\hat{d}$ and set
$\textsc{dict}[\theta]=[x_{\theta},f_{\theta}(x_{\theta}),\nabla{f}_{\theta}(x_{\theta})]$.

Step 6. Moving forward and termination check: Set
$x\leftarrow{x}-\alpha\hat{d}$. If certain termination criteria are satisfied
at $x$, then the algorithm terminates; otherwise, go to Step 2.

Table 2: Test problems [18] and their properties

function name | analytical form
---|---
gen_MAXQ | $\max_{1\leq{i}\leq{n}}\;x^{2}_{i}$
gen_MXHILB | $\max_{1\leq{i}\leq{n}}\left|\frac{x_{i}}{i+j-1}\right|$
Chained_LQ | $\sum^{n-1}_{i=1}\max\\{-x_{i}-x_{i+1},-x_{i}-x_{i+1}+(x^{2}_{i}+x^{2}_{i+1}-1)\\}$
Chained_CB3_I | $\sum^{n-1}_{i=1}\\{x^{4}_{i}+x^{2}_{i+1},(2-x_{i})^{2}+(2-x_{i+1})^{2},2e^{-x_{i}+x_{i+1}}\\}$
Chained_CB3_II | $\sum^{n-1}_{i=1}\max\\{x^{2}_{i}+(x_{i+1}-1)^{2}+x_{i+1}-1,-x^{2}_{i}-(x_{i+1}-1)^{2}+x_{i+1}+1\\}$
num_active_faces | $\max_{1\leq{i}\leq{n}}\left\\{\log\left|\sum^{n}_{i=1}x_{i}\right|,\log|x_{i}|\right\\}$
brown_func2 | $\sum^{n-1}_{i=1}\big{(}|x_{i}|^{x^{2}_{i+1}+1}+|x_{i+1}|^{x^{2}_{i}+1}\big{)}$
Chained_Mifflin2 | $\sum^{n-1}_{i=1}(-x_{i}+2(x^{2}_{i}+x^{2}_{i+1}-1)+1.75|x^{2}_{i}+x^{2}_{i+1}-1|$
Chained_Crescent_I | $\max\left\\{\sum^{n-1}_{i=1}(x^{2}_{i}+(x_{i+1}-1)^{2}+x_{i+1}-1),\sum^{n-1}_{i=1}(-x^{2}_{i}-(x_{i+1}-1)^{2}+x_{i+1}+1)\right\\}$
Chained_Crescent_II | $\sum^{n-1}_{i=1}\max\\{x^{2}_{i}+(x_{i+1}-1)^{2}+x_{i+1}-1,-x^{2}_{i}-(x_{i+1}-1)^{2}+x_{i+1}+1\\}$

function name | optimal objective | Convexity | uniform domain
---|---|---|---
gen_MAXQ | 0.0 | $+$ | $+$
gen_MXHILB | 0.0 | $+$ | $+$
Chained_LQ | $-\sqrt{2}(n-1)$ | $+$ | $+$
Chained_CB3_I | $2(n-1)$ | $+$ | $+$
Chained_CB3_II | $2(n-1)$ | $+$ | $+$
num_active_faces | 0.0 | $-$ | $-$
brown_func2 | 0.0 | $-$ | $-$
Chained_Mifflin2 |  unclear∗ | $-$ | $+$
Chained_Crescent_I | 0.0 | $-$ | $+$
Chained_Crescent_II | 0.0 | $-$ | $+$

∗$f(x^{*})\approx-6.51$ for $n=10$, $f(x^{*})\approx-70.15$ for $n=100$, and
$f(x^{*})\approx-706.55$ for $n=1000$.

Table 3: Performance of non-smooth optimization methods on median-size problem instances function name | dim | time(JGD) | obj(JGD) | time(TRB) | obj(TRB) | time(QN) | obj(QN) | time(VM) | obj(VM)
---|---|---|---|---|---|---|---|---|---
gen_MAXQ | 50 | 3.47 | 3.011E-03 | 3.49 | 2.305E-01 | - | 9.492E-01 | - | 9.357E-01
gen_MAXQ | 100 | 11.39 | 1.224E-04 | 13.70 | 2.384E-01 | 2.31 | 9.553E-01 | - | 9.606E-01
gen_MXHILB | 50 | 1.43 | 1.112E-04 | 15.91 | 2.100E-02 | 1.45 | 1.120E-01 | 7.45 | 1.552E-02
gen_MXHILB | 100 | 0.67 | 9.929E-02 | 3.41 | 1.091E-01 | 1.66 | 8.559E-02 | 8.57 | 5.592E-02
Chained_LQ | 50 | - | 2.483E-03 | 27.07 | 2.837E+00 | 10.08 | 3.712E+01 | - | 5.902E-01
Chained_LQ | 100 | - | 2.017E-02 | 243.54 | 2.810E+00 | 11.90 | 7.118E+01 | - | 2.373E+00
Chained_CB3_I | 50 | 48.17 | 1.784E-05 | 613.86 | 7.054E+00 | 18.08 | 6.998E+01 | - | 2.908E+01
Chained_CB3_I | 100 | 150.65 | 4.408E-01 | - | 2.096E+00 | - | 1.122E+02 | - | 1.499E+02
Chained_CB3_II | 50 | 1.77 | 1.018E-06 | 250.54 | 4.028E+00 | 4.39 | 2.444E+01 | - | 7.854E-01
Chained_CB3_II | 100 | 0.61 | 9.966E-08 | - | 1.493E+00 | 5.20 | 3.682E+01 | - | 1.930E+00
num_active_faces | 50 | 15.57 | 4.413E-05 | 43.36 | 4.546E-01 | 3.84 | 5.558E-01 | - | 4.166E-01
num_active_faces | 100 | 6.59 | 9.375E-05 | 1.59 | 1.094E+00 | 2.64 | 6.860E-01 | - | 5.082E-01
brown_func2 | 50 | 9.91 | 3.279E-04 | - | 1.438E+00 | 832.63 | 5.375E-03 | 82.58 | 4.326E+00
brown_func2 | 100 | 5.90 | 2.516E-04 | - | 1.278E+00 | - | 3.345E+00 | 24.01 | 3.028E+01
Chained_Mifflin2 | 50 | 69.44 | -3.476E+01 | 34.39 | -3.348E+01 | 13.39 | -2.240E+01 | - | -2.957E+01
Chained_Mifflin2 | 100 | 44.67 | -7.000E+01 | 89.64 | -6.917E+01 | - | -4.892E+01 | - | -5.408E+01
Chained_Crescent_I | 50 | 0.44 | 2.869E-08 | 34.62 | 2.891E+00 | 3.01 | 3.592E+01 | - | 3.857E+00
Chained_Crescent_I | 100 | 0.46 | 4.954E-07 | 789.75 | 7.449E-01 | 4.24 | 6.288E+01 | - | 1.408E+00
Chained_Crescent_II | 50 | 21.10 | 1.882E-05 | 66.26 | 1.785E+00 | - | 4.817E+00 | - | 1.725E+00
Chained_Crescent_II | 100 | 19.67 | 2.881E-02 | 235.80 | 2.485E+00 | - | 3.706E+01 | - | 6.338E+00
Table 4: Performance of the joint gradient descent method on solving large-size problem instances. For columns of objective value checkpoints, the symbol ‘-’ indicates that the termination time is shorter than that checkpoint time.  function name | dim | obj(init) | term. time (sec) | obj(final) | obj(1min) | obj(5min) | obj(10min) | components
---|---|---|---|---|---|---|---|---
gen_MAXQ | 500 | 9.950E-01 | 158.69 | 2.228E-04 | 9.525E-03 | - | - | 497
gen_MAXQ | 1000 | 9.995E-01 | 570.74 | 3.444E-04 | 2.460E-01 | 4.802E-03 | - | 990
gen_MAXQ | 2000 | 9.995E-01 | 1200.02 | 1.876E-02 | 5.565E-01 | 1.645E-01 | 6.840E-02 | 1736
gen_MAXQ | 5000 | 9.995E-01 | 1201.13 | 4.829E-01 | 8.733E-01 | 7.137E-01 | 6.128E-01 | 1503
gen_MXHILB | 500 | 1.328E+00 | 0.77 | 4.655E-04 | - | - | - | 44
gen_MXHILB | 1000 | 1.343E+00 | 3.00 | 3.334E-04 | - | - | - | 54
gen_MXHILB | 2000 | 1.318E+00 | 9.16 | 4.143E-04 | - | - | - | 110
gen_MXHILB | 5000 | 1.316E+00 | 57.78 | 2.234E-04 | - | - | - | 113
Chained_LQ | 500 | 7.141E+02 | 4729.99 | 1.534E-01 | 1.907E-01 | 1.570E-01 | 1.534E-01 | 2769
Chained_LQ | 1000 | 1.468E+03 | 12646.27 | 2.095E-01 | 3.157E-01 | 2.096E-01 | 2.096E-01 | 1222
Chained_LQ | 2000 | 2.868E+03 | 1200.26 | 5.463E-01 | 1.755E+00 | 9.489E-01 | 6.637E-01 | 3861
Chained_LQ | 5000 | 7.260E+03 | 1201.88 | 1.732E+00 | 4.344E+00 | 2.632E+00 | 2.195E+00 | 1933
Chained_CB3_I | 500 | 3.286E+03 | 128.51 | 1.390E+00 | 1.392E+00 | - | - | 2082
Chained_CB3_I | 1000 | 6.723E+03 | 348.09 | 1.083E+00 | 1.162E+00 | 1.083E+00 | - | 3103
Chained_CB3_I | 2000 | 1.317E+04 | 727.93 | 1.446E+00 | 1.833E+00 | 1.459E+00 | 1.447E+00 | 2925
Chained_CB3_I | 5000 | 3.323E+04 | 1200.61 | 1.557E+00 | 2.616E+00 | 1.928E+00 | 1.675E+00 | 1535
Chained_CB3_II | 500 | 3.259E+03 | 1.05 | 7.715E-09 | - | - | - | 3
Chained_CB3_II | 1000 | 6.667E+03 | 2.19 | 4.829E-09 | - | - | - | 3
Chained_CB3_II | 2000 | 1.306E+04 | 5.52 | 1.121E-08 | - | - | - | 3
Chained_CB3_II | 5000 | 3.298E+04 | 27.75 | 1.334E-06 | - | - | - | 3
num_active_faces | 500 | 2.286E+00 | 39.98 | 3.499E-01 | - | - | - | 1215
num_active_faces | 1000 | 1.024E+00 | 58.49 | 4.976E-01 | - | - | - | 818
num_active_faces | 2000 | 3.603E+00 | 105.31 | 5.999E-01 | 6.003E-01 | - | - | 604
num_active_faces | 5000 | 3.636E+00 | 888.69 | 6.511E-01 | 6.946E-01 | 6.915E-01 | 6.865E-01 | 729
brown_func2 | 500 | 4.415E+02 | 20.89 | 7.558E-04 | - | - | - | 2055
brown_func2 | 1000 | 8.693E+02 | 112.39 | 5.314E-04 | 1.284E-03 | - | - | 2930
brown_func2 | 2000 | 1.753E+03 | 173.31 | 6.742E-04 | 1.628E+00 | - | - | 2741
brown_func2 | 5000 | 4.314E+03 | 909.55 | 2.513E-03 | 6.112E+02 | 5.994E+00 | 9.356E-02 | 2847
Chained_Mifflin2 | 500 | 4.378E+01 | 436.39 | -3.524E+02 | -3.519E+02 | -3.524E+02 | - | 8833
Chained_Mifflin2 | 1000 | 1.024E+02 | 1200.21 | -7.052E+02 | -7.025E+02 | -7.049E+02 | -7.052E+02 | 13750
Chained_Mifflin2 | 2000 | 1.977E+02 | 1200.07 | -1.411E+03 | -1.408E+03 | -1.410E+03 | -1.411E+03 | 6668
Chained_Mifflin2 | 5000 | 4.594E+02 | 1201.23 | -3.527E+03 | -3.514E+03 | -3.522E+03 | -3.525E+03 | 4152
Chained_Crescent_I | 500 | 3.271E+02 | 0.78 | 8.950E-08 | - | - | - | 2
Chained_Crescent_I | 1000 | 6.641E+02 | 1.70 | 2.128E-09 | - | - | - | 2
Chained_Crescent_I | 2000 | 1.311E+03 | 2.54 | 1.880E-09 | - | - | - | 2
Chained_Crescent_I | 5000 | 3.244E+03 | 10.83 | 1.826E-06 | - | - | - | 2
Chained_Crescent_II | 500 | 5.312E+02 | 122.62 | 3.236E-02 | 3.302E-02 | - | - | 2479
Chained_Crescent_II | 1000 | 1.052E+03 | 649.33 | 1.654E-01 | 2.741E+00 | 6.928E-01 | 1.655E-01 | 7748
Chained_Crescent_II | 2000 | 2.117E+03 | 1200.09 | 8.867E-01 | 5.277E+00 | 1.812E+00 | 1.027E+00 | 7318
Chained_Crescent_II | 5000 | 5.238E+03 | 1201.34 | 4.886E+00 | 9.379E+00 | 7.649E+00 | 6.221E+00 | 2968

Figure 4: Performance curves of 4 non-smooth optimization methods on solving
median size problem instances. Figure 5: The efficiency of reducing objective
values compared for 4 non-smooth optimization methods. The dimension of all
problem instances is 100. Figure 6: Plot of objective over time generated by
the JGD method for selected objective functions and $n=500,1000,2000,5000$.
Figure 7: Violin plots of the number of component functions used in computing
the joint gradient (referred as the number of active component functions in
the figure) in each iteration.

##  Appendix A Omitted proofs

###### proof of Proposition 1.

(a) Consider the following quadratic program:

$\min_{\lambda}\;\left\lVert\sum^{k}_{i=1}\lambda_{i}a_{i}\right\rVert^{2}\quad\textrm{s.t.
}\;\sum^{k}_{i=1}\lambda_{i}=1,$ (A.1)

which can be written equivalently as: $\min_{\lambda}\;\left\lVert
A\lambda\right\rVert^{2}\quad\textrm{s.t. }\;e^{\top}\lambda=1$. Using the
Lagrangian multiplier method, one can verify that the optimal solution
$\lambda^{*}$ of the above quadratic program is

$\lambda^{*}=\frac{1}{e^{\top}(A^{\top}A)^{-1}e}\cdot(A^{\top}A)^{-1}e,$ (A.2)

and we have $d^{*}=A\lambda^{*}$. Furthermore, we have

$[\left\langle a_{1},d^{*}\right\rangle,\left\langle
a_{2},d^{*}\right\rangle,\ldots,\left\langle
a_{k},d^{*}\right\rangle]^{\top}=A^{\top}d^{*}=\frac{e}{e^{\top}(A^{\top}A)^{-1}e},$

which implies that

$\left\langle
a_{i},d^{*}\right\rangle=\frac{1}{e^{\top}(A^{\top}A)^{-1}e}\quad\forall{i\in[k]}.$

It means that moving along the direction $d_{i}$, each $f_{i}$ has the same
amount of increment.

(b) We compute the decrement value of $f$ along any normalized direction
$-\hat{d}$ at $x_{0}$, where $\|\hat{d}\|=1$. First, $\hat{d}$ can be uniquely
decomposed as $\hat{d}=d^{\prime}+\alpha\hat{d}^{*}$ for some $\alpha\leq 1$
such that $\langle d^{\prime},\hat{d}^{*}\rangle=0$, where
$\hat{d}^{*}=d^{*}/\left\lVert d^{*}\right\rVert$. The decrement of $f$ along
$-\hat{d}$ with step-size $\delta$ can be computed as

$\displaystyle f(x_{0})-f(x_{0}-\delta\hat{d})$
$\displaystyle=\max\\{a^{\top}_{i}x_{0}+b_{i}:\forall
i\in[k]\\}-\max\\{a^{\top}_{i}(x_{0}-\delta\hat{d})+b_{i}:\forall i\in[k]\\}$
$\displaystyle=f(x_{0})-\left(f(x_{0})-\min\\{\delta
a^{\top}_{i}\hat{d}:\forall i\in[k]\\}\right)$ $\displaystyle=\min\\{\delta
a^{\top}_{i}\hat{d}:\forall i\in[k]\\}.$

Note that we have

$\displaystyle\langle a_{i},\hat{d}\rangle$ $\displaystyle=\alpha\langle
a_{i},\hat{d}^{*}\rangle+\langle a_{i},d^{\prime}\rangle$
$\displaystyle=\alpha\langle a_{i},\hat{d}^{*}\rangle+\langle
a_{i}-d^{*},d^{\prime}\rangle+\langle d^{*},d^{\prime}\rangle$
$\displaystyle=\alpha\left\lVert d^{*}\right\rVert+\left\langle
a_{i}-d^{*},d^{\prime}\right\rangle.$

Since we have $\sum^{k}_{i=1}\lambda^{*}_{i}\left\langle
a_{i}-d^{*},d^{\prime}\right\rangle=0$ by the property of $d^{*}$ and
$\lambda^{*}$, there exists an index $j$ such that $\left\langle
a_{j}-d^{*},d^{\prime}\right\rangle\leq 0$ given that $\lambda^{*}\geq 0$.
Using this property, we further compute

$\displaystyle f(x_{0})-f(x_{0}-\delta\hat{d})$
$\displaystyle=\delta\alpha\left\lVert
d^{*}\right\rVert+\min\\{\delta\left\langle
a_{i}-d^{*},d^{\prime}\right\rangle:\forall i\in[k]\\}$
$\displaystyle\leq\delta\alpha\left\lVert
d^{*}\right\rVert\leq\delta\left\lVert d^{*}\right\rVert,$

where the equality holds if $\alpha=1$, i.e., $\hat{d}=\hat{d}^{*}$. This
indicates that indeed, $-d^{*}$ is the steepest descent direction of $f$ at
$x_{0}$. ∎

###### proof of Proposition 2.

Let $u=\sum_{i\in\mathcal{I}}\lambda^{*}_{i}v_{i}$ and consider the following
two-variate problem for any $i\in\mathcal{I}^{\prime}$

$\min_{\alpha_{1},\alpha_{2}}\left\lVert\alpha_{1}u+\alpha_{2}v_{i}\right\rVert^{2}\quad\mathrm{s.t.}\;\alpha_{1}+\alpha_{2}=1,\;\alpha_{1},\alpha_{2}\geq
0.$ (A.3)

By the optimality condition of (2.5), the above problem must have an optimal
solution $[\alpha^{*}_{1},\alpha^{*}_{2}]=[1,0]$, otherwise $\lambda^{*}$
cannot be the optimal solution of (2.5). The Lagrangian of (A.3) is given by

$L(\alpha,\beta)=\left\lVert\alpha_{1}u+\alpha_{2}v_{i}\right\rVert^{2}+\beta(\alpha_{1}+\alpha_{2}-1)-\tau_{1}\alpha_{1}-\tau_{2}\alpha_{2}.$

The KKT conditions of (A.3) read as

$\displaystyle 2\alpha_{1}\left\lVert
u\right\rVert^{2}+2\alpha_{2}\left\langle
u,v_{i}\right\rangle+\beta-\tau_{1}=0,$ $\displaystyle 2\alpha_{2}\left\lVert
v_{i}\right\rVert^{2}+2\alpha_{1}\left\langle
u,v_{i}\right\rangle+\beta-\tau_{2}=0,$
$\displaystyle\alpha_{1}+\alpha_{2}=1,\;\alpha_{1}\tau_{1}=0,\;\alpha_{2}\tau_{2}=0,\;\alpha_{1},\alpha_{2},\tau_{1},\tau_{2}\geq
0,$

and they lead to the following condition after eliminating $\alpha_{1}$ and
$\beta$

$(\left\lVert u\right\rVert-\left\lVert
v_{i}\right\rVert)^{2}\alpha_{2}+\frac{\tau_{1}-\tau_{2}}{2}=\left\lVert
u\right\rVert^{2}-\left\langle u,v_{i}\right\rangle.$

Since $[\alpha^{*}_{1},\alpha^{*}_{2}]=[1,0]$ is a solution of the KKT system,
the complementary slackness condition implies that $\tau_{1}=0$. Substituting
$\alpha^{*}_{2}=0$ and $\tau_{1}=0$ into the above equation gives

$\left\lVert u\right\rVert^{2}-\left\langle
u,v_{i}\right\rangle=-\frac{\tau_{2}}{2}\leq 0,$

which concludes the proof. ∎

######  Observation 2.

For any continuous function $g$, if $x^{*}$ is the unique local minimizer of
$g$ the set $U(v):=\\{x:g(x)\leq{v}\\}$ for $v\geq{g(x^{*})}$ is a connected
closed set containing $x^{*}$

###### Proof.

Suppose there exists a value $v\geq{g(x^{*})}$ such that $U(v)$ is a closed
set consisting of multiple disjointed (closed) subsets $V_{1},\ldots,V_{k}$
($k\geq 2$). Suppose $x^{*}\in{V}_{1}$. By definition of the level sets, there
must exist a local minimizer in $V_{2}$ which is different from $x^{*}$. This
contradicts to that $x^{*}$ is the unique local minimizer of $g$. ∎

######  Proposition 10.

For any $x_{0}\in X$, the function
$v(x)=\left\lVert\mathrm{vec(}\mathrm{QP}\textrm{-}{x}\textrm{-}{\mathcal{C}}\mathrm{)}\right\rVert$
is continuous at $x_{0}$ in a neighborhood $B$ of $x_{0}$ satisfying
$\mathcal{C}\subseteq\widetilde{\mathcal{M}}(x)$ for every $x\in{B}$.

###### Proof.

By definition
$\mathrm{vec(}\mathrm{QP}\textrm{-}{x}\textrm{-}{\mathcal{C}}\mathrm{)}$ is an
optimal vector of the following QP:

$\min_{\lambda}\;\left\lVert\sum_{\theta\in\mathcal{C}}\lambda_{\theta}\nabla{f_{\theta}(x)}\right\rVert\;\textrm{s.t.
}\sum_{\theta\in\mathcal{C}}\lambda_{\theta}=1,\;\lambda\geq 0.$ (A.4)

Let $\lambda^{*}$ be an optimal solution at $x_{0}$ of the above QP. The
continuity of $v(x)$ at $x_{0}$ follows from the fact that
$\nabla{f_{\theta}}$ is continuous at $x_{0}$ for $\theta\in\mathcal{C}$. ∎

## References

  * [1] Felipe Atenas, Claudia Sagastizábal, Paulo J. S. Silva, and Mikhail Solodov.  A unified analysis of descent sequences in weakly convex optimization, including convergence rates for bundle methods.  SIAM Journal on Optimization, 33(1):89–115, 2023.
  * [2] Adil Bagirov, Napsu Karmitsa, and Marko M. Mäkelä.  Bundle methods.  In Introduction to Nonsmooth Optimization, chapter 12, pages 305–310. Springer International Publishing, Switzerland, 2014.
  * [3] James V. Burke, Adrian S. Lewis, and Michael L. Overton.  A robust gradient sampling algorithm for nonsmooth, nonconvex optimization.  SIAM Journal on Optimization, 15(3):751–779, 2005.
  * [4] Vasileios Charisopoulos and Damek Davis.  A superlinearly convergent subgradient method for sharp semismooth problems.  Mathematics of Operations Research, in press.
  * [5] Ying Cui and Jong-Shi Pang.  Modern nonconvex nondifferentiable optimization.  MOS-SIAM Series on Optimization. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2021.
  * [6] Frank E. Curtis and Minhan Li.  Gradient sampling methods with inexact subproblem solves and gradient aggregation.  arXiv, 2021.
  * [7] Frank E. Curtis and Xiaocun Que.  A quasi-newton algorithm for nonconvex, nonsmooth optimization with global convergence guarantees.  Mathematical Programming Computation, 7:399–428, 2015.
  * [8] Minh N. Dao.  Bundle method for nonconvex nonsmooth constrained optimization.  Journal of Convex Analysis, 22:1061–1090, 2015.
  * [9] Minh N. Dao, Joachim Gwinner, Dominikus Noll, and Nina Ovcharova.  Nonconvex bundle method with application to a delamination problem.  Computational Optimization and Applications, 65:173–203, 2016.
  * [10] Mateo Díaz and Benjamin Grimmer.  Optimal convergence rates for the proximal bundle method.  SIAM Journal on Optimization, 33:424–454, 2023.
  * [11] Yu Du and Andrzej Ruszczyński.  Rate of convergence of the bundle method.  Journal of Optimization Theory and Applications, 173:908–922, 2017\.
  * [12] Sabrina Fiege, Andrea Walther, and Andreas Griewank.  An algorithm for nonsmooth optimization by successive piecewise linearization.  Mathematical Programming, 177:343–370, 2019.
  * [13] Sabrina Fiege, Andrea Walther, Kshitij Kulshreshtha, and Andreas Griewank.  Algorithmic differentiation for piecewise smooth functions: a case study for robust optimization.  Optimization Methods and Software, 33(4-6):1073–1088, 2018.
  * [14] Andreas Griewank.  On stable piecewise linearization and generalized algorithmic differentiation.  Optimization Methods and Software, 28:1139–1178, 2013.
  * [15] Andreas Griewank and Andrea Walther.  First- and second-order optimality conditions for piecewise smooth objective functions.  Optimization Methods and Software, 31(5):904–930, 2016.
  * [16] Andreas Griewank and Andrea Walther.  Relaxing kink qualifications and proving convergence rates in piecewise smooth optimization.  SIAM Journal on Optimization, 29(1):262–289, 2019.
  * [17] Andreas Griewank, Andrea Walther, Sabrina Fiege, and Torsten Bosse.  On Lipschitz optimization based on gray-box piecewise linearization.  Mathematical Programming, 158:383–415, 2016.
  * [18] M. Haarala, K. Miettinen, and M. M. Makela.  New limited memory bundle method for large-scale nonsmooth optimization.  Optimization Methods and Software, 19(6):673–692, 2004.
  * [19] X. Y. Han and Adrian S. Lewis.  Survey descent: a multipoint generalization of gradient descent for nonsmooth optimization.  SIAM Journal on Optimization, 33(1):36–62, 2023.
  * [20] W. Hare, C. Sagastizábal, and M. Solodov.  A proximal bundle method for nonsmooth nonconvex functions with inexact information.  Computational Optimization and Applications, 63:1–28, 2015.
  * [21] Warren Hare and Claudia Sagastizábal.  A redistributed proximal bundle method for nonconvex optimization.  SIAM Journal on Optimization, 20(5):2442–2473, 2010.
  * [22] Michael Hintermüller.  A proximal bundle method based on approximate subgradients.  Computational Optimization and Applications, 20:245–266, 2001.
  * [23] Krzysztof C. Kiwiel.  Restricted step and Levenberg-Marquardt techniques in proximal bundle methods for nonconvex nondifferentiable optimization.  SIAM Journal on Optimization, 6(1):227–249, 1996.
  * [24] Krzysztof C. Kiwiel.  A proximal bundle method with approximate subgradient linearization.  SIAM Journal on Optimization, 16(4):1007–1023, 2006.
  * [25] Krzysztof C. Kiwiel.  Convergence of the gradient sampling algorithm for nonsmooth nonconvex optimization.  SIAM Journal on Optimization, 18(2):379–388, 2007.
  * [26] C. Lemaréchal.  Bundle methods in nonsmooth optimization.  In C. Lemaréchal and R. Mifflin, editors, Nonsmooth Optimization, pages 79–102, Laxenburg, 1977. Pergamon Press, Oxford.
  * [27] Adrian S. Lewis and Michael L. Overton.  Nonsmooth optimization via quasi-Newton methods.  Mathematical Programming, 141:135–163, 2013.
  * [28] Ladislav Lukšan and Jan Vlček.  A bundle-Newton method for nonsmooth unconstrained minimization.  Mathematical Programming, 83:373–391, 1998.
  * [29] Jian Lv, Li-Ping Pang, and Fan-Yun Meng.  A proximal bundle method for constrained nonsmooth nonconvex optimization with inexact information.  Journal of Global Optimization, 70:517–549, 2018.
  * [30] Morteza Maleknia and Mostafa Shamsi.  A quasi-Newton proximal bundle method using gradient sampling technique for minimizing nonsmooth convex functions.  Optimization Methods and Software, 37(4):1415–1446, 2022.
  * [31] Robert Mifflin, Defeng Sun, and Liqun Qi.  Quasi-Newton bundle-type methods for nondifferentiable convex optimization.  SIAM Journal on Optimization, 8(2):583–603, 1998.
  * [32] N. Hoseini Monjezi and S. Nobakhtian.  Convergence of the proximal bundle algorithm for nonsmooth nonconvex optimization problems.  Optimization Letter, 16:1495–1511, 2022.
  * [33] Najmeh H. Monjezi.  A bundle trust region algorithm for minimizing locally Lipschitz functions.  SIAM Journal on Optimization, 33(1):319–337, 2023.
  * [34] Najmeh H. Monjezi and S. Nobakhtian.  A filter proximal bundle method for nonsmooth nonconvex constrained optimization.  Journal of Global Optimization, 79:1–37, 2021.
  * [35] Najmeh Hoseini Monjezi and S. Nobakhtian.  A new infeasible proximal bundle algorithm for nonsmooth nonconvex constrained optimization.  Computational Optimization and Applications, 74:443–480, 2019.
  * [36] Claudia Sagastizábal.  Composite proximal bundle method.  Mathematical Programming, 140:189–233, 2013.
  * [37] J. Vlcek and L. Luksan.  Globally convergent variable metric method for nonconvex nondifferentiable unconstrained minimization.  Journal of Optimization Theory and Applications, 111(2):407–430, 2001.
  * [38] Andreas Wächter and Lorenz T. Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Mathematical Programming, 106(6):25–57, 2006.


remarkRemark hypothesisHypothesis claimClaim Bounding the approach to
oligarchyDavid W. Cohen and Bruce M. Boghosian

# Bounding the approach to oligarchy in a variant of the yard-sale model

David W. Cohen Department of Mathematics, Tufts University, Medford, MA (,
https://sites.tufts.edu/davidcohen/).The first author acknowledges the support
of the NDSEG fellowship<EMAIL_ADDRESS>Bruce M. Boghosian Department
of Mathematics, Tufts University, Medford, MA<EMAIL_ADDRESS>
###### Abstract

We present analytical results for the Gini coefficient of economic inequality
under the dynamics of a modified Yard-Sale Model of kinetic asset exchange. A
variant of the Yard-Sale Model is introduced by modifying the underlying
binary transaction of the classical system. It is shown that the Gini
coefficient is monotone under the resulting dynamics but the approach to
oligarchy, as measured by the Gini index, can be bounded by a first-order
differential inequality used in conjunction with the differential Grönwall
inequality. The asymptotics of the modified system, with a redistributive tax,
are derived and shown to agree with the original, taxed Yard-Sale Model, which
implies the modified system is as suitable for matching real wealth
distributions.

###### keywords:

yard-sale model, econophysics, non-linear Fokker-Planck equation, mean-field
theory, McKean-Vlasov equations

91B80, 82C22, 82C31

##  1 Introduction

The Yard-Sale Model is a well-studied model of kinetic asset exchange
introduced by A. Chakraborti and named by B. Hayes [13, 19]. At its core, the
Yard-Sale Model is a specification of how a wealth transaction is to be
conducted between two agents randomly selected from a population. From this
point, the system may be altered and then studied as a stochastic finite-agent
system or, through the use of thermodynamic limits and other techniques from
mathematical physics, as a deterministic, continuum equation of motion for the
distribution of wealth.

The classical Yard-Sale Model transaction is stochastic: At each integer time
$t$, two agents from an $N$ agent population are selected at random without
replacement and their wealths are updated according to the rule

(1) $\begin{pmatrix}w_{t+1}^{i}\\\
w_{t+1}^{j}\end{pmatrix}=\begin{pmatrix}w_{t}^{i}\\\
w_{t}^{j}\end{pmatrix}+\sqrt{\gamma}\left(w_{t}^{i}\wedge
w_{t}^{j}\right)\begin{pmatrix}1\\\ -1\end{pmatrix}\eta,$

where $\gamma\in(0,1)$ is a transaction intensity parameter, $\eta$ is a
random variable with outcomes $-1$ and $+1$ with equal probability, and
$\wedge$ is the $\min$ operator.

Under this rule the expected change of wealth for each agent is zero yet it is
well known that wealth condenses and an oligarchy forms as time progresses.
The definition of the classical transaction rule is premised on the belief
that the poorer agent’s wealth should determine the magnitude of a binary
transaction and that the changes in wealth are associated to a lack of perfect
information in a given transaction.111Informally, in a transaction with a very
wealthy agent, a poorer agent is willing to stake a small fraction of their
wealth but cannot possibly stake the same fraction of the other’s wealth as
losing would horrifically bankrupt the poorer agent. Adhering to this
principle, the Yard-Sale Model is the simplest, non-trivial stochastic
transaction that cannot send an agent to negative wealth when all agents are
initialized with positive wealth.

In [3] a time step is introduced into the transaction Eq. 1 to give an
infinitesimal characterization of the process by a limiting procedure
involving the number of agents and time step combined with an analogue of the
molecular chaos assumption (Stosszahlansatz) – termed the random-agent
approximation. This procedure produces a deterministic continuum equation. The
equation of motion for the probability distribution of agents in wealth-space
under the classical Yard-Sale Model dynamics of Eq. 1 is

(2) $\frac{\partial\rho(w,t)}{\partial t}=\frac{\partial^{2}}{\partial
w^{2}}\left[\left(\frac{\gamma}{2}\int_{0}^{\infty}\,dx\,\left(w\wedge
x\right)^{2}\rho(x,t)\right)\rho(w,t)\right].$

A frequently studied summary statistic of distributions of wealth is the Gini
coefficient of economic inequality, which was introduced by Corrado Gini in
1912 [18]. The Gini coefficient maps a wealth distribution to a value in
$[0,1]$. Values near zero correspond to more egalitarian distributions,
whereas values of the Gini approaching unity indicate inequality and
oligarchy. If $\mu$ is the mean wealth of a distribution of wealth with
density $\rho$ and $X,Y$ are i.i.d. random variables with density $\rho$ then

$G[\rho]:=\frac{\mathbb{E}\left[|X-Y|\right]}{2\mu}.$

Boghosian et al. in [6] established that the Gini coefficient is monotone
under the dynamics of the continuum model Eq. 2 of the classical Yard-Sale
Model. This result was shown both for the master equation and the resulting
non-linear partial integro-differential Fokker-Planck equation first derived
in [3]. Thus the Gini coefficient is a Lyapunov functional [17, 20] for the
Yard-Sale Model. Chorro showed via a martingale convergence theorem argument
in [14] that the finite-agent system likewise approaches oligarchy (in a
probabilistic sense). Börgers and Greengard in [7] produced yet simpler proofs
of wealth condensation for the classical finite-agent system. Börgers and
Greengard [7] and Boghosian [5] have similar results for transactions that are
biased in favor of the wealthier agent. Beyond the results for the Yard-Sale
Model and its variants, Cardoso et al. in [9, 10] have shown that wealth
condensation is more likely the rule than the exception for more general
unbiased binary exchanges.

Despite the many methods of showing that wealth condensation occurs, we are
not aware of any explicit bounds on the rate of increase of the Gini
coefficient under models for which the Gini coefficient increases
monotonically.

Cao in [8] investigated the rate of change of the Gini coefficient under the
dynamics of the repeated averaging model (sometimes called the divorce model)
where the binary transaction sends each agent to the mean wealth of the two
transacting agents. Cao mentions that in “econophysics literature, analytical
results on [the] Gini index are comparatively rare.”

Here we produce an analytic result on the Gini coefficient for a modified, yet
reasonable, Yard-Sale Model that steps beyond the standard propositions of
monotonicity. In particular, the result is a bound on the rate of production
of inequality, as measured by the Gini coefficient, under the modified
dynamics. From a physical perspective, this is akin to finding a bound on the
rate of entropy production while still having a thermodynamic second law.

The paper is organized as follows. The variant of the Yard-Sale Model on which
the present paper focuses is motivated and defined in Section 2. The Gini
coefficient is briefly reviewed in Section 3 with particular focus on its
invariance under a normalization of the equations of motion. In Section 4 it
is proven both that the Gini coefficient increases monotonically in time under
the induced dynamics and that its rate of increase may be bounded. Finally the
asymptotics of the modified system when a redistributive tax is incorporated
are derived in Section 5 and shown to match the classical Yard-Sale Model with
taxation.

##  2 The modified Yard-Sale Model

Henceforth we sometimes use the abbreviation YSM for the yard-sale model.

Let there be $N>1$ agents each with a dimensionally-meaningful wealth
$\theta^{i}_{s}$ indexed by $i=1,\ldots,N$ and a time $s\geq 0$. Time
subscripts $s$ are occasionally omitted for clarity when the time is
unimportant for the expression. Let $W_{\theta}:=\sum\theta^{i}$ be the total
wealth and $\mu_{\theta}:=W_{\theta}/N$ be the average wealth per agent. To
each of $N$ agents, associate a dimensionless quantity
$w^{i}:=\theta^{i}/\mu_{\theta}$ obtained by dividing wealth by mean
wealth.222In what follows, we still call $w$ wealth despite its dimensionless
nature.

Let $k=0,1,2,\ldots$ and $\Delta t\in(0,1)$. The modified version of the Yard-
Sale Model has a binary transaction between agents indexed by $i$ and $j$
given by

(3) $\begin{pmatrix}w_{(k+1)\Delta t}^{i}\\\ w_{(k+1)\Delta
t}^{j}\end{pmatrix}=\begin{pmatrix}w_{k\Delta t}^{i}\\\ w_{k\Delta
t}^{j}\end{pmatrix}+\sqrt{\gamma\Delta t}\phi\left(w^{i}_{k\Delta
t},w^{j}_{k\Delta t}\right)\begin{pmatrix}1\\\ -1\end{pmatrix}\eta,$

with

(4) $\phi\left(w^{i},w^{j}\right)=\begin{cases}(w^{i}\wedge w^{j})&\text{\quad
if $w^{i}\wedge w^{j}<1$;}\\\ \sqrt{(w^{i}\wedge w^{j})}&\text{\quad if
$w^{i}\wedge w^{j}\geq 1$,}\end{cases}$

where $\gamma\in(0,1)$, $\wedge$ is the $\min$ operator, and $\eta$ is a
random variable that takes values $-1$ and $+1$ with equal probability. The
standard YSM has $w^{i}\wedge w^{j}$ outside of the square root function
regardless of the size of $w^{i}\wedge w^{j}$. The importance of the square
root in Eq. 4 will become clearer below in light of the more general statement
of Corollary 4.7 and the relation of the diffusion coefficient kernel to that
of the Gini coefficient.

This modification maintains the important features that the poorer agent’s
wealth is the determining quantity in the exchange and no exchange can send an
agent to negative wealth. Note that the comparison of dimensionless wealth in
the piecewise micro-transaction is equivalent to checking if the minimum of
the two dimension-full wealths is above or below mean wealth $\mu_{\theta}$.

Under the same assumptions common in mathematical physics that are laid out in
[2] and [3] – namely assuming independence of the laws of each agent as the
number of agents goes to infinity – the limit as $\Delta t\rightarrow 0$ and
$N\rightarrow\infty$ leads to a continuum equation for the law of a single
prototypical agent in the population. The nonlinear evolution equation is
found via the Kramers-Moyal expansion of the Chapman-Kolmogorov equation after
making the random-agent approximation [3] and applying Pawula’s Theorem. A
more formal derivation from the stochastic process would use methods from the
propagation of chaos literature on interacting particle systems of the
Boltzmann type, see [11, 12, 21] for more details. The yard-sale model and its
variants result in McKean-Vlasov stochastic differential equations [1] for
which the drift and/or diffusion coefficients are functionals of the dependent
variable.

Let $\rho(w,t)$ be the probability density of agents in the dimensionless
wealth variable. Since the dimensionless quantity $w$ was obtained by dividing
wealth $\theta$ by mean wealth $\mu_{\theta}$, the first moment of $\rho(w,t)$
is also $1$, that is $\int_{\mathbb{R}^{+}}w\rho=1$. Therefore, by
construction, the mean (dimensionless) wealth is $1$.

This leads to a Fokker-Planck equation for the normalized agent distribution

(5) $\frac{\partial\rho(w,t)}{\partial t}=\frac{\partial^{2}}{\partial
w^{2}}\left[\underbrace{\frac{\gamma}{2}\left(\int_{0}^{\infty}\,dx\,\kappa(w,x)\rho(x,t)\right)}_{:=D[w,\rho(\cdot,t)]}\rho(w,t)\right],$

where

(6) $\kappa(w,x):=\begin{cases}(w\wedge x)^{2}&\text{\quad if $w\wedge
x<1$;}\\\ (w\wedge x)&\text{\quad if $w\wedge x\geq 1$.}\end{cases}$

Both the zeroth and first moments of $\rho$ are conserved quantities; the
former being the probability mass (implying conservation of number of agents)
and the latter corresponding to conservation of total wealth, which we take
canonically to be 1.

Throughout this paper we assume that each agent has positive wealth and the
support of all distributions are $\mathbb{R}^{+}$.

##  3 The Gini coefficient under normalization

Here for the reader’s convenience we define the Gini coefficient and show its
invariance under a particular transformation of a distribution of wealth.

Let $Q:[0,\infty)\rightarrow[0,\infty)$ have finite zeroth and first moment
denoted $N_{Q}$ and $W_{Q}$, respectively. Define $\mu_{Q}=W_{Q}/N_{Q}$. The
transformation to a distribution with unit zeroth and first moment is given by

$q(w):=\frac{\mu_{Q}}{N_{Q}}Q\left(\mu_{Q}w\right).$

$Q$ can be viewed as a wealth distribution of a population of $N_{Q}$ agents
with total population wealth $W_{Q}$. The Gini coefficient of economic
inequality can be expressed as

$\displaystyle G[Q]$
$\displaystyle=1-\frac{1}{N_{Q}W_{Q}}\int_{0}^{\infty}\,dw\,\int_{0}^{\infty}\,dx\,\left(w\wedge
x\right)Q(w)Q(x)$ (7)
$\displaystyle=1-\int_{0}^{\infty}\,dw\,\int_{0}^{\infty}\,dx\,\left(w\wedge
x\right)q(w)q(x),$

which we can equivalently call $G[q]$ under the normalizing transformation
$Q\mapsto q$. Thus $G$ is invariant under the transformation between $q$ and
$Q$.

The Gini coefficient is $0$ for a density concentrated on mean wealth (that
is, for a wealth-egalitarian society) whereas it approaches its upper limit of
$1$ as the wealth is concentrated into an ever-vanishing proportion of the
population. See [4, 15] for a discussion of the nonstandard properties of
wealth distributions that extremize the Gini coefficient under the dynamics of
Eq. 2.

We will work in the space of normalized wealth distributions in which both the
zeroth and first moments are unity.

######  Lemma 3.1.

Two derivatives of the Fréchet derivative of the Gini coefficient under
normalization is twice the density, that is

(8) $\frac{d^{2}}{d{w}^{2}}\frac{\delta G}{\delta\rho}=2\rho(w).$

This is a routine calculation and follows from the symmetry of the kernel in
the double integral that defines $G[\rho]$. We state the proof, despite its
simplicity, since the result is repeatedly used.

######  Proof 3.2.

Since the integral kernel of $G$ is symmetric in its arguments, the Fréchet
derivative of $G$ is

$\frac{\delta G}{\delta\rho}=-2\int_{\mathbb{R}^{+}}\,dx\,\left(w\wedge
x\right)\rho(x).$

The first $w$-derivative of the Fréchet derivative is

$\frac{d}{dw}\frac{\delta G}{\delta\rho}=-2\int_{w}^{\infty}\,dx\,\rho(x)$

and the next $w$-derivative is

$\frac{d^{2}}{d{w}^{2}}\frac{\delta G}{\delta\rho}=2\rho(w).$

##  4 Bounding the rate of increase of the Gini coefficient

We now turn to the main results of the paper: That the Gini coefficient,
despite being monotonically increasing under the modified Yard-Sale Model
dynamics, can have a non-trivial bound on its rate of change in time. This
bound may be carried over into a bound on the value of the Gini coefficient at
a future time.

By $G(t)$, we mean $G[\rho(\cdot,t)]$ where $\rho$ is a solution to Eq. 5.

######  Theorem 4.1 (Increasing inequality).

The Gini coefficient Eq. 7 is monotone increasing under the dynamics of Eq. 5.

######  Proof 4.2.

$\displaystyle\frac{dG}{dt}$
$\displaystyle=\int_{0}^{\infty}\,dw\,\frac{\delta
G}{\delta\rho}\frac{\partial\rho}{\partial t}$
$\displaystyle=\int_{0}^{\infty}\,dw\,\frac{\delta
G}{\delta\rho}\frac{\partial^{2}}{\partial
w^{2}}\left[D\left[w,\rho(\cdot,t)\right]\rho(w,t)\right]$ (by Eq. 5)
$\displaystyle=\int_{0}^{\infty}\,dw\,\left(\frac{\partial^{2}}{\partial
w^{2}}\frac{\delta
G}{\delta\rho}\right)D\left[w,\rho(\cdot,t)\right]\rho(w,t)$ (via two IbP)
$\displaystyle=2\int_{0}^{\infty}\,dw\,D\left[w,\rho(\cdot,t)\right]\rho(w,t)^{2}$
(by Lemma 3.1) $\displaystyle\geq 0$ (since $D\left[w,\rho\right]\geq 0$)

The above result can be strengthened when there is mass that is not
concentrated at the origin.

######  Corollary 4.3 (Strictly increasing inequality).

Let $s>0$ and $\epsilon\in(0,1).$ If there exists $a>0$ such that

$\int_{a}^{\infty}\,dw\,\rho(w,s)>\epsilon$

then there exists $\delta>0$ such that

$\frac{dG}{dt}\bigg{|}_{t=s}>\delta.$

######  Proof 4.4.

Let $C=(a^{2}\wedge a)$, which is positive. Making use of the calculation in
Theorem 4.1, we have

$\displaystyle\frac{dG}{dt}\bigg{|}_{t=s}$
$\displaystyle=\int_{\mathbb{R}^{+}}\,dw\,D[w,\rho(s,\cdot)]\left(\rho(s,w)\right)^{2}$
$\displaystyle\geq\int_{a}^{\infty}\,dw\,\int_{a}^{\infty}\,dx\,\kappa(w,x)\rho(s,x)\left(\rho(s,w)\right)^{2}$
$\displaystyle\geq C\epsilon\int_{a}^{\infty}\,dw\,\left(\rho(s,w)\right)^{2}$
$\displaystyle\geq C\epsilon\frac{\epsilon^{2}}{4(b-a)}$
$\displaystyle(\text{Jensen's inequality})$ $\displaystyle>0.$

The step prior to the application of Jensen’s inequality uses that there must
exist $b$ such that $0<a<b<\infty$ for which

$\int_{a}^{b}\,dx\,\rho(s,x)>\frac{\epsilon}{2}.$

Letting

$\delta=\frac{C\epsilon^{3}}{4(b-a)}>0$

completes the proof.

Corollary 4.3 also holds for the classical Yard-Sale Model.

######  Theorem 4.5 (Bounding inequality).

If the initial datum $\rho_{0}(w)=\rho(w,0)$ is in
$L^{\infty}(\mathbb{R}^{+})$ and the dynamics keep $\rho(w,t)$ in
$L^{\infty}(\mathbb{R}^{+})$ up to time $s>0$ then the rate of change of the
Gini coefficient Eq. 7 is bounded by

(9) $\frac{dG}{dt}\leq\gamma||\rho(\cdot,t)||_{\infty}\left(1-G(t)\right)$

for $0<t<s.$

######  Proof 4.6.

Let $t\in(0,s).$ Starting from the penultimate line in the main computation of
the proof of Theorem 4.1, we have

$\displaystyle\frac{dG}{dt}$
$\displaystyle=2\int_{0}^{\infty}\,dw\,D[w,\rho]\left(\rho(w)\right)^{2}$
$\displaystyle\leq 2||\rho||_{\infty}\int_{0}^{\infty}\,dw\,D[w,\rho]\rho(w)$
(Hölder’s inequality)
$\displaystyle=2||\rho||_{\infty}\int_{0}^{\infty}\,dw\,\rho(w)\frac{\gamma}{2}\int_{0}^{\infty}\,dx\,\kappa(w,x)\rho(x)$
(definition of $D[w,\rho]$)
$\displaystyle\leq\gamma||\rho||_{\infty}\int_{0}^{\infty}\,dw\,\rho(w)\int_{0}^{\infty}\,dx\,(w\wedge
x)\rho(x)$ ($\kappa(w,x)\leq(w\wedge x)$ )
$\displaystyle=\gamma||\rho||_{\infty}(1-G)$ $\displaystyle\text{\quad(by
\lx@cref{creftype~refnum}{eq:gini})}.$

Note that both Theorem 4.1 and Theorem 4.5 can be extended to other binary
transactions that modify the original Yard-Sale Model if the transaction
kernel meets two conditions.

######  Corollary 4.7.

For a system with transaction kernel
$\widetilde{\phi}:\mathbb{R}^{+}\times\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$
such that

$\widetilde{\phi}\left(w^{i},w^{j}\right)\leq\left(w^{i}\wedge w^{j}\right)$

and

$\widetilde{\phi}\left(w^{i},w^{j}\right)^{2}\leq\left(w^{i}\wedge
w^{j}\right)$

on $\mathbb{R}^{+}\times\mathbb{R}^{+}$, we have that

$0\leq\frac{dG}{dt}\leq\gamma||\rho||_{\infty}(1-G).$

The two conditions on the kernel correspond to not sending agents to negative
wealth in the finite model for $\Delta t\in(0,1]$ and the fourth step in the
proof of Theorem 4.5, respectively. However, the class of kernels satisfying
these conditions need not have the stake of a transaction be independent of
the wealthier agent.

The lower bound implies that the Gini coefficient increases under the dynamics
of the modified YSM, whereas the upper bound on $\dot{G}$ bounds $G(t)$.

######  Remark 4.8.

If the upper bound of Theorem 4.5 were in fact attained as
$\dot{G}=\gamma||\rho||_{\infty}(1-G)$, it could be viewed as asserting that
at a lower inequality the accessible rearrangements of wealth lead quickly to
increasing inequality. This speculative remark seems to agree with numerical
simulations of the dynamics in which the random initial condition quickly
changes in conjunction with a rapid increase of the Gini coefficient followed
by a slowing of the change in the distribution and the Gini coefficient.

To produce the bound on $G(t)$ we need a slightly modified version of the
differential Grönwall’s inequality as presented in [16].

######  Lemma 4.9 (A version of Grönwall’s differential inequality).

Let $\eta,\phi,$ and $\psi$ be functions of $t\in[0,T]$. Further let $\eta$ be
absolutely continuous and $\phi$ and $\psi$ be integrable on $[0,T]$. If for
almost every $t$

$\dot{\eta}\leq\phi\eta+\psi,$

then

$\eta(t)\leq\exp\left(\int_{0}^{t}\,ds\,\phi(s)\right)\left[\eta(0)+\int_{0}^{t}\,ds\,\psi(s)\exp\left(-\int_{0}^{s}\,dr\,\phi(r)\right)\right].$

######  Proof 4.10.

By direct calculation first and then applying the assumed inequality, we have
that

$\displaystyle\frac{d}{ds}\left[\eta(s)\exp\left(-\int_{0}^{s}\,dr\,\phi(r)\right)\right]$
$\displaystyle=\exp\left(-\int_{0}^{s}\,dr\,\phi(r)\right)\left[\dot{\eta}(s)-\phi(s)\eta(s)\right]$
$\displaystyle\leq\exp\left(-\int_{0}^{s}\,dr\,\phi(r)\right)\psi(s).$

Integrating from $0$ to $t$ in $s$ yields

$\eta(t)\exp\left(-\int_{0}^{t}\,dr\,\phi(r)\right)-\eta(0)\leq\int_{0}^{t}\,ds\,\psi(s)\exp\left(-\int_{0}^{s}\,dr\,\phi(r)\right),$

which upon re-arranging proves the lemma.

If $\phi$ were assumed to be nonnegative then
$\exp\left(-\int_{0}^{s}\,dr\,\phi(r)\right)$ is at most unity for all $s$, in
which case the upper bound on that term by unity gives the standard
inequality.

######  Claim 1.

Let

$M_{T}:=\sup_{t\in[0,T]}||\rho(\cdot,t)||_{\infty}.$

For $t\in(0,T]$,

$G(t)\leq G(0)\exp\left(-\gamma M_{T}t\right)+\frac{\exp\left(\gamma
M_{T}t\right)-1}{\exp\left(\gamma M_{T}t\right)}.$

######  Proof 4.11.

Apply Lemma 4.9 with $\eta=G$, $\phi=-\gamma M_{T}$, and $\psi=\gamma M_{T}$.

Thus the magnitude of the exponential rate of convergence to oligarchy under
the modified YSM is at most $\gamma M_{T}$ in a finite time horizon $[0,T]$.

######  Claim 2.

Let $M(t):=||\rho(\cdot,t)||_{\infty}$. For $t>0$,

$G(t)\leq
G(0)\exp\left(-\gamma\int_{0}^{t}\,ds\,M(s)\right)+\frac{\exp\left(\gamma\int_{0}^{t}\,ds\,M(s)\right)-1}{\exp\left(\gamma\int_{0}^{t}\,ds\,M(s)\right)}.$

######  Proof 4.12.

Apply Lemma 4.9 with $\eta=G$, $\phi=-\gamma M(s)$, and $\psi=\gamma M(s)$.
Note that the term

$\int_{0}^{t}\,ds\,\psi(s)\exp\left(-\int_{0}^{s}\,dr\,\phi(r)\right)$

simplifies to

$\exp\left(\gamma\int_{0}^{t}\,ds\,M(s)\right)-1.$

##  5 Asymptotic analysis

Theorem 4.1 asserts that inequality monotonically increases under the dynamics
of Eq. 5, which leads to total wealth condensation. The present state of
affairs in the world is not one of total wealth condensation so in order for
this model to have explanatory or predictive power, it is necessary to
introduce a regularization term and ask if the resulting wealth distributions
are close to those of the world.

In [5], real distributions of wealth were compared to those arising
asymptotically from the YSM with differing combinations of: (1) taxation and
redistribution, (2) a bias in the stochastic transaction in favor of the
wealthier agent, and (3) negative wealth. It was found that the model’s
asymptotics match very well with real wealth distributions. Thus, rather than
carrying out such comparisons from scratch, we instead show that the
asymptotics of the modified YSM match those of the classical model.

To regularize the modified YSM, we introduce a deterministic taxation and
redistribution term that reverts the population to the mean wealth. This
scheme can be viewed in two equivalent ways:

  * •

As an exogenous force that collects the same fraction of each agent’s wealth,
pools the total, and redistributes the total equally amongst the agents; or

  * •

As a partial divorce model so that in addition to each stochastic transaction
between agents there is also a deterministic exchange towards their pairwise
mean wealth.

These two perspectives yield the same continuum equation of motion.

The binary transaction for the modified YSM with redistribution is

(10) $\begin{pmatrix}w_{(k+1)\Delta t}^{i}\\\ w_{(k+1)\Delta
t}^{j}\end{pmatrix}=\begin{pmatrix}w_{k\Delta t}^{i}\\\ w_{k\Delta
t}^{j}\end{pmatrix}+\left[\chi\Delta t\left(w^{j}_{k\Delta t}-w^{i}_{k\Delta
t}\right)+\sqrt{\gamma\Delta t}\phi\left(w^{i}_{k\Delta t},w^{j}_{k\Delta
t}\right)\eta\right]\begin{pmatrix}1\\\ -1\end{pmatrix},$

where $\chi\in[0,\frac{1}{2})$ and the other definitions and parameters from
Eq. 3 and Eq. 4 carry over.

The associated equation of motion for the agent density in wealth-space is

(11) $\frac{\partial\rho(w,t)}{\partial t}=-\frac{\partial}{\partial
w}\left[\chi(1-w)\rho(w,t)\right]+\frac{\partial^{2}}{\partial
w^{2}}\left[\frac{\gamma}{2}\left(\int_{0}^{\infty}\,dx\,\kappa(w,x)\rho(x,t)\right)\rho(w,t)\right],$

where $\kappa$ is as before.333Starting from Eq. 11, both $\gamma$ and $\chi$
lose their initial restrictions imposed from the binary transactions and we
demand only their non-negativitity.

In this section, we closely follow the procedure of asymptotic analysis laid
out in [5]. In what follows, $\rho_{\infty}(w)$ is the asymptotic state of Eq.
11 for which the time derivative is set to zero and the ODE is studied.

###  5.1 Analysis for equilibrium at $w\ll 1$

For $w\ll 1$, $\kappa(w,x)=\left(w\wedge x\right)^{2}$ thus

$D[w,\rho]=\frac{\gamma}{2}\left(\int_{0}^{w}\,dx\,x^{2}\rho(x,t)+w^{2}\int_{w}^{\infty}\,dx\,\rho(x,t)\right).$

We approximate this by $D[w,\rho]\approx\frac{\gamma}{2}w^{2}$. This
assumption demands the same a posteriori justification as noted in [5]. At
equilibrium

$\chi(1-w)\rho_{\infty}(w)=\frac{d}{dw}\left[\frac{\gamma}{2}w^{2}\rho_{\infty}(w)\right].$

This is solved by

$\rho_{\infty}(w)=\frac{c_{0}}{w^{2+2\chi/\gamma}}\exp\left(-\frac{2\chi}{\gamma
w}\right),$

where $c_{0}$ is a positive constant. This agrees with the results in [5]
where this analysis was carried out for the classic Yard-Sale Model with
redistribution. Thus the justification for the earlier assumption about the
behavior of $D[w,\rho]$ at $w\ll 1$ is equally valid as in the original paper.

###  5.2 Analysis for equilibrium at $w\gg 1$

Upon integrating once and re-arranging, the large-$w$ equilibrium condition is

(12)
$\frac{d\log\rho_{\infty}(w)}{dw}=\frac{\chi(1-w)-\frac{dD[w,\rho_{\infty}]}{dw}}{D[w,\rho_{\infty}]}.$

We first investigate the behavior of $D[w,\rho]$ for $w\gg 1$. The piecewise
conditions of $\kappa(w,x)$ partially lose their dependence on $w$ so that

$\displaystyle D[w,\rho]$
$\displaystyle=\frac{\gamma}{2}\int_{0}^{\infty}\,dx\,\rho(x)\begin{cases}x^{2}&\text{\quad
if $x<1$;}\\\ (w\wedge x)&\text{\quad if $x\geq 1$,}\end{cases}$
$\displaystyle=\frac{\gamma}{2}\left(\int_{0}^{1}\,dx\,x^{2}\rho(x)+\int_{1}^{w}\,dx\,x\rho(x)+w\int_{w}^{\infty}\,dx\,\rho(x)\right).$

Let

$c_{1}:=\int_{0}^{1}\,dx\,x^{2}\rho(x)$

and

$c_{2}:=\int_{0}^{\infty}\,dx\,x\rho(x)-\int_{0}^{1}\,dx\,x\rho(x),$

noting that $c_{2}>0$. Thus

(13)
$D[w,\rho]=\frac{\gamma}{2}\left(c_{1}+c_{2}-\int_{w}^{\infty}\,dx\,x\rho(x)+w\int_{w}^{\infty}\,dx\,\rho(x)\right).$

Therefore we also have

(14) $\frac{dD[w,\rho]}{dw}=\frac{\gamma}{2}\int_{w}^{\infty}\,dx\,\rho(x).$

We make the ansatz that

(15) $\rho_{\infty}(w)\approx c_{\infty}\exp(-aw^{2}-bw)$

for $w\gg 1$, $a>0$, $b\in\mathbb{R}$, and $c_{\infty}>0$. In this case,

(16) $\frac{d\log\rho_{\infty}(w)}{dw}=-2aw-b.$

We simplify the right hand side of Eq. 12 term-by-term under the ansatz Eq. 15
to check whether it is an affine function of $w$, and if so identify the
appropriate constants.

In doing so, we make repeated use of the asymptotic approximation

$\text{erfc}(z)\approx\frac{1}{z\sqrt{\pi}}\exp(-z^{2}).$

Starting first with Eq. 14

$\displaystyle\int_{w}^{\infty}\,dx\,\rho_{\infty}(x)$
$\displaystyle=\frac{c_{\infty}}{2}\sqrt{\frac{\pi}{a}}\exp\left(\frac{b^{2}}{4a}\right)\text{erfc}\left(\frac{2aw+b}{2\sqrt{a}}\right)$
$\displaystyle\approx\frac{c_{\infty}}{2}\sqrt{\frac{\pi}{a}}\exp\left(\frac{b^{2}}{4a}\right)\frac{2\sqrt{a}}{2aw+b}\frac{1}{\sqrt{\pi}}\exp\left(-\left(\frac{2aw+b}{2\sqrt{a}}\right)^{2}\right)$
(17) $\displaystyle=\frac{c_{\infty}}{2aw+b}\exp(-aw^{2}-bw).$

Next we turn to the terms with $w$ dependence in Eq. 13. The first of which is

$\displaystyle\int_{w}^{\infty}\,dx\,x\rho_{\infty}(x)$
$\displaystyle=c_{\infty}\frac{\exp\left(-aw^{2}-bw\right)}{2a}-\frac{b}{2a}\left[\frac{c_{\infty}}{2}\sqrt{\frac{\pi}{a}}\exp\left(\frac{b^{2}}{4a}\right)\text{erfc}\left(\frac{2aw+b}{2\sqrt{a}}\right)\right]$
$\displaystyle\approx
c_{\infty}\frac{\exp\left(-aw^{2}-bw\right)}{2a}-\frac{b}{2a}\left[\frac{c_{\infty}}{2aw+b}\exp(-aw^{2}-bw)\right]$
$\displaystyle=\frac{c_{\infty}}{2a}\exp(-aw^{2}-bw)\left(1-\frac{b}{2aw+b}\right)$
(18) $\displaystyle=w\frac{c_{\infty}}{2aw+b}\exp(-aw^{2}-bw).$

Combining Eqs. 17 and 18, we have that

$-\int_{w}^{\infty}\,dx\,x\rho_{\infty}(x)+w\int_{w}^{\infty}\,dx\,\rho_{\infty}(x)\approx
0,$

which implies $D[w,\rho]\approx\frac{\gamma}{2}(c_{1}+c_{2}).$

Taking these results together, we see that upon using Eq. 15 and approximating
asymptotically, Eq. 12 simplifies to

$\frac{d\log\rho_{\infty}(w)}{dw}=\frac{2\chi}{\gamma(c_{1}+c_{2})}(1-w)-\underbrace{\frac{c_{\infty}}{c_{1}+c_{2}}\frac{\exp\left(-aw^{2}-bw\right)}{2aw+b}}_{\text{subdominant
corrections}}.$

In particular, we note from Eq. 16 that
$a=\frac{\chi}{\gamma(c_{1}+c_{2})}>0$.

##  6 Conclusions

We introduced a variant of the Yard-Sale Model for which the Gini coefficient
of economic inequality monotonically increases under the resulting continuum
dynamics yet the rate of change in time of the Gini coefficient permits an
upper bound.

The method of bounding $\dot{G}$ used techniques from the analysis of
deterministic equation of motions. It is not clear if or how this bound could
be applied to the stochastic, finite-agent models that precede the diffusion
approximation. Nor is it clear how a similar bound, phrased in a probabilistic
way, could be derived directly from the stochastic, finite-agent system.

Since the bound in Theorem 4.5 includes $||\rho||_{\infty}$, the $L^{\infty}$
norm of the density, it is natural to ask how the bound could be tested with
empirical data since the transition from empirical data to a density function
involves user choices (e.g., width in kernel density estimation) that affect
the essential supremum of the resulting density.

By showing that the asymptotics of this modified model with redistribution
match that of the original Yard-Sale Model with redistribution, we have put
forward an argument that the macroscopic asymptotics are robust to changes of
the microscopic transactions. Thus many microscopic transactions (and not just
a practitioner’s favorite) may give rise to simple yet accurate descriptions
of the evolution of wealth that depend on dramatically fewer parameters than
most economic theories.

Finally, the upper bound on $\dot{G}$ was shown to hold not only for the
particular modification studied but also for a broader class of models
described in Corollary 4.7. In doing so, we can consider classes of kinetic
asset exchange models with “entropy production bounds” where the phrase should
instead be considered as “inequality production bounds.”

## Acknowledgments

We thank M. Johnson and D. Gentile of Tufts University for productive
conversations.

## References

  * [1] V. Bogachev, N. Krylov, M. Röckner, and S. Shaposhnikov, Fokker–Planck–Kolmogorov Equations, Mathematical Surveys and Monographs, American Mathematical Society, 2015.
  * [2] B. M. Boghosian, Fokker-Planck Description of Wealth Dynamics and the Origin of Pareto’s Law, International Journal of Modern Physics C, 25 (2014).
  * [3] B. M. Boghosian, Kinetics of Wealth and the Pareto Law, Physical Review E, 89 (2014).
  * [4] B. M. Boghosian and C. Börgers, The Mathematics of Poverty, Inequality, and Oligarchy, SIAM News, (2023).
  * [5] B. M. Boghosian, A. Devitt-Lee, M. Johnson, J. Li, J. A. Marcq, and H. Wang, Oligarchy as a Phase Transition: The Effect of Wealth-Attained Advantage in a Fokker-Planck Description of Asset Exchange, Physica A, 476 (2017).
  * [6] B. M. Boghosian, M. Johnson, and J. A. Marcq, An H-Theorem for Boltzmann’s Equation for the Yard-Sale Model of Asset Exchange, Journal of Statistical Physics, 161 (2014).
  * [7] C. Börgers and C. Greengard, A new probabilistic analysis of the Yard-Sale Model, 2023, https://arxiv.org/abs/2308.01485.
  * [8] F. Cao, Explicit Decay Rate for the Gini Index in the Repeated Averaging Model, Mathematical Methods in the Applied Sciences, 46 (2023).
  * [9] B. H. F. Cardoso, S. Gonçalves, and J. R. Iglesias, Why Equal Opportunities Lead to Maximum Inequality? The Wealth Condensation Paradox Generally Solved, Chaos, Solitions, and Fractals, 168 (2023).
  * [10] B. H. F. Cardoso, J. R. Iglesias, and S. Gonçalves, Wealth Concentration in Systems with Unbiased Binary Exchanges, Physica A, 579 (2021).
  * [11] L. P. Chaintron and A. Diez, Propagation of chaos: A review of models, methods and applications. I. Models and methods, Kinetic and Related Models, 15 (2022), p. 895, https://doi.org/10.3934/krm.2022017, https://doi.org/10.3934%2Fkrm.2022017.
  * [12] L. P. Chaintron and A. Diez, Propagation of chaos: A review of models, methods and applications. II. Applications, Kinetic and Related Models, 15 (2022), p. 1017, https://doi.org/10.3934/krm.2022018, https://doi.org/10.3934%2Fkrm.2022018.
  * [13] A. Chakroborti, Distributions of Money in Model Markets of Economy, International Journal of Modern Physics C, 13 (2002).
  * [14] C. Chorro, A Simple Probabilistic Approach of the Yard-Sale Model, Statistics & Probability Letters, 112 (2016).
  * [15] A. Devitt-Lee, H. Wang, J. Li, and B. M. Boghosian, A Nonstandard Description of Wealth Concentration in Large-Scale Economies, SIAM Journal of Applied Mathematics, 78 (2018).
  * [16] L. C. Evans, Partial Differential Equations, American Mathematical Society, 2 ed., 2010.
  * [17] T. D. Frank, Nonlinear Fokker-Planck Equations, Springer-Verlag, 2010.
  * [18] C. Gini, Variabilità e mutabilità: contributo allo studio delle distribuzioni e delle relazioni statistiche, Tipogr. di P. Cuppini, Bologna, Italy, 1912.
  * [19] B. Hayes, Follow the Money, American Scientist, 90 (2002).
  * [20] A. N. Michel, L. Hou, and D. Liu, Stability of Dynamical Systems: On the Role of Monotonic and Non-Monotonic Lyapunov Functions, Birkhäuser, 2 ed., 2015.
  * [21] A. S. Sznitman, Topics in Propagation of Chaos, Springer-Verlag, 1991.