
\subsection{Convergence}
We begin by showing that the algorithm's termination states are \aclp{NE}:

\begin{theorem}
\label{thm:stationary}
Assume that Algorithm \ref{alg:MXL} is run with
a decreasing step-size sequence $\gamma_{n}$ such that $\sum_{n=1}^{\infty} \gamma_{n}^{2} < \sum_{n=1}^{\infty} \gamma_{n} = \infty$
and
gradient observations satisfying \eqref{eq:zeromean} and \eqref{eq:MSE}.
If $\mathbf{X}(n)$ converges, it does so to a \acl{NE} of the game \textup(a.s.\textup)\xspace.
\end{theorem}


The proof of Theorem \ref{thm:stationary} is presented in detail in Appendix \ref{app:stationary} and is essentially by contradiction.
To provide some intuition, if the limit of \eqref{eq:MXL} is not a \acl{NE}, at least one player of the game must be dissatisfied,  thus experiencing a repelling drift \textendash\ in the limit and on average.
Owing to the algorithm's exponentiation step, this ``repulsion'' can be quantified via the so-called \emph{von Neumann} (or \emph{quantum}) \emph{entropy} \cite{Ved02};
then, by using the law of large numbers, it is possible to control the impact of the noise and show that this entropy diverges to infinity, thus obtaining the desired contradiction.

\smallskip
Stated concisely, Theorem \ref{thm:stationary} shows that if \eqref{eq:MXL} converges, it converges to a \acl{NE}.
However, the theorem does not provide any guarantee that the algorithm converges in the first place. A sufficient condition for convergence based on equilibrium stability is given below:

\begin{theorem}
\label{thm:global}
Assume that Algorithm \ref{alg:MXL} is run with
a step-size sequence $\gamma_{n}$ such that $\sum_{n=1}^{\infty} \gamma_{n}^{2} < \sum_{n=1}^{\infty} \gamma_{n} = \infty$
and
gradient observations satisfying \eqref{eq:zeromean} and \eqref{eq:MSE}.
If $\act^{\ast}$ is globally stable, then $\mathbf{X}(n)$ converges to $\act^{\ast}$ \textup(a.s.\textup)\xspace.
\end{theorem}

The proof of Theorem \ref{thm:global} is given in Appendix \ref{app:global}.
In short, it comprises the following steps:
First, we consider a deterministic, continuous-time variant of \eqref{eq:MXL} and we show that globally stable equilibria are global attractors of said system.
To do this, we introduce a matrix version of the so-called \emph{Fenchel coupling} \cite{MS16} and we show that it plays the role of a Lyapunov function in continuous time.
Subsequently, we derive the evolution of the discrete-time stochastic system \eqref{eq:MXL} by using the method of stochastic approximation \cite{Ben99} and the theory of concentration inequalities \cite{dlP99} to control the aggregate error between continuous and discrete time.

\smallskip

From a practical point of view, an immediate corollary of Theorem \ref{thm:global} is the following second-derivative convergence test:

\begin{corollary}
\label{cor:global}
If the game's Hessian matrix $\bD(\mathbf{X})$ is negative-definite for all $\mathbf{X}\in\boldsymbol{\mathcal{X}}$, Algorithm \ref{alg:MXL} converges to the game's \textup(necessarily\textup) unique \acl{NE} \textup(a.s.\textup)\xspace.
\end{corollary}

The above results show that \eqref{eq:MXL} converges to stable equilibria under very mild assumptions on the underlying stochasticity (zero-mean errors and finite conditional variance).
The following theorem shows that, under the slightly sharper assumption \eqref{eq:subexp}, local stability implies local convergence with arbitrarily high probability:

\begin{theorem}
\label{thm:local}
Assume that Algorithm \ref{alg:MXL} is run with
gradient observations satisfying \eqref{eq:zeromean} and \eqref{eq:subexp},
and small enough step-sizes $\gamma_{n}$ such that $\sum_{n=1}^{\infty} \gamma_{n}^{1+q} < \sum_{n=1}^{\infty} \gamma_{n} = \infty$ for some $q\in(0,1)$.
If $\act^{\ast}$ is \textup(locally\textup) stable, then it is locally attracting with arbitrarily high probability;
specifically, for every $\varepsilon>0$, there exists a neighborhood $U_{\varepsilon}$ of $\act^{\ast}$ such that
\begin{equation}
\label{eq:prob-local}
\probof{\mathbf{X}(n) \to \act^{\ast} \given \mathbf{X}(0) \in U_{\varepsilon}}
	\geq 1-\varepsilon.
\end{equation}
\end{theorem}

Theorem \ref{thm:local} is proven in Appendix \ref{app:local}, building on the (global) proof of Theorem \ref{thm:global}.
The key difference with Theorem \ref{thm:global} is that, since we only assume the existence of a \emph{locally stable} equilibrium $\act^{\ast}$, the drift of \eqref{eq:MXL} has no reason to point towards $\act^{\ast}$ globally.
As a result, $\mathbf{X}(n)$ may exit the basin of $\act^{\ast}$ in the presence of high uncertainty.
However, by invoking \eqref{eq:subexp} and the Borel\textendash Cantelli lemma, it is possible to show that this happens with arbitrarily small probability, leading to the probabilistic convergence result \eqref{eq:prob-local}.


As in the case of globally stable equilibria, Theorem \ref{thm:local} leads to the following easy-to-check condition for local convergence:
\begin{corollary}
\label{cor:local}
Let $\act^{\ast}$ be a \acl{NE} of $\mathcal{G}$ such that $\bD(\act^{\ast}) \prec 0$.
Then \eqref{eq:MXL} converges locally to $\act^{\ast}$ with arbitrarily high probability.
\end{corollary}



\subsection{Rate of convergence}

Combining Theorems \ref{thm:stationary}\textendash \ref{thm:local} gives a fairly complete picture of the \emph{qualitative} convergence properties of the exponential learning algorithm \eqref{eq:MXL}:
the only possible end-states of \eqref{eq:MXL} are \aclp{NE}, and the existence of a globally (resp.~locally) stable \acl{NE} implies the algorithm's global (resp.~local) convergence to said equilibrium. On the other hand, these theorems do not address the \emph{quantitative} aspects of the algorithm's long-term behavior, such as its rate of convergence.
In what follows, we study precisely this question.

We begin by introducing the so-called \emph{quantum Kullback\textendash Leibler divergence} (or \emph{von Neumann relative entropy}) to measure distances on $\boldsymbol{\mathcal{X}}$ \cite{Ved02}.
Specifically, the quantum divergence between $\mathbf{X}$ and $\act^{\ast}$ is defined as
\begin{equation}
\label{eq:Bregman}
D_{\mathrm{KL}}(\act^{\ast},\mathbf{X})
	= \trof{\act^{\ast} (\log\act^{\ast} - \log\mathbf{X})}.
\end{equation}
By Klein's inequality \cite{Ved02}, $D_{\mathrm{KL}}(\act^{\ast},\mathbf{X}) \geq 0$ with equality if and only if $\mathbf{X}=\act^{\ast}$, so $D_{\textup{KL}}(\act^{\ast},\mathbf{X})$ represents a (convex) measure of the distance between $\mathbf{X}$ and $\act^{\ast}$.
With this in mind, we introduce below the following quantitative measure of stability:
\begin{definition}
\label{def:stable-strong}
Given $B>0$, $\act^{\ast}$ is called \emph{$B$-strongly stable} if
\begin{equation}
\label{eq:stable-strong}
\trof{(\mathbf{X} - \act^{\ast})\,\bV(\mathbf{X})}
	\leq -B\,D_{\mathrm{KL}}(\act^{\ast},\mathbf{X})
	\quad
	\text{for all $\mathbf{X}\in\boldsymbol{\mathcal{X}}$}.
\end{equation}
\end{definition}

Under this refinement of equilibrium stability,%
\footnote{Since $D_{\mathrm{KL}}(\act^{\ast},\mathbf{X}) \geq 0$ with equality if and only if $\mathbf{X} = \act^{\ast}$, it follows that strongly stable states are automatically stable.}
we obtain the following quantitative result:

\begin{theorem}
\label{thm:convrate}
Assume that Algorithm \ref{alg:MXL} is run with the step-size sequence $\gamma_{n} = \gamma/n$
and
gradient observations satisfying \eqref{eq:zeromean} and \eqref{eq:MSE}.
If $\act^{\ast}$ is $B$-strongly stable and $\gamma > B^{-1}$, we have
\begin{equation}
\label{eq:convrate}
\exof{D_{\mathrm{KL}}(\act^{\ast},\mathbf{X}(n))}
	\leq \frac{\gamma^{2} V^{2}}{B\gamma - 1} \frac{1}{n}.
\end{equation}
\begin{subequations}
In particular, if $\act^{\ast}$ lies in the interior of $\boldsymbol{\mathcal{X}}$, then
\begin{equation}
\label{eq:convrate-int}
\exof{\norm{\mathbf{X}(n) - \act^{\ast}}}
	= \bigoh(n^{-1/2}).
\end{equation}
Otherwise, if $\act^{\ast}$ is an extreme point of $\boldsymbol{\mathcal{X}}$, we have
\begin{equation}
\label{eq:convrate-ext}
\exof{\norm{\mathbf{X}(n) - \act^{\ast}}}
	= \bigoh(n^{-1}).
\end{equation}
\end{subequations}
\end{theorem}

The explicit bounds of Theorem \ref{thm:convrate} (proven in Appendix \ref{app:rates}) have several interesting consequences:

\begin{enumerate}
[\itshape i\upshape)]
\item
The rate of convergence to interior versus extreme states of $\boldsymbol{\mathcal{X}}$ is substantially different because the von Neumann divergence grows quadratically when $\act^{\ast}$ is interior and linearly when $\act^{\ast}$ is extreme.
Intuitively, the reason for this is that, in the case of interior states, the algorithm has to slow down considerably near $\act^{\ast}$ in order for random oscillations around $\act^{\ast}$ to dissipate.
On the other hand, in the case of extreme points, the algorithm has a constant drift towards the boundary of $\boldsymbol{\mathcal{X}}$ and there are no possibilities of overshooting (precisely because the equilibrium state in question is an extreme point of $\boldsymbol{\mathcal{X}}$).
As a result, random oscillations become irrelevant near extreme points, thus allowing for faster convergence in that case. 

\smallskip

\item
Thanks to the explicit expression \eqref{eq:convrate}, we see that if the stability constant $B$ can be estimated ahead of time, we can attain the convergence rate
\begin{equation}
\label{eq:convrate-opt}
\exof{D_{\mathrm{KL}}(\act^{\ast},\mathbf{X}(n))}
	\leq \frac{4V^{2}}{B^{2}n},
\end{equation}
achieved for the optimized step-size sequence $\gamma_{n} = 2/(Bn)$.
It is also possible to estimate the algorithm's convergence rate for equilibrium states that are only locally strongly stable.
However, given that the algorithm's convergence is probabilistic in that case, the resulting bounds are also probabilistic and, hence, more complicated to present.
\end{enumerate}

\subsection{Practical implementation}
\label{par:implementation}

We close this section with a discussion of certain issues pertaining to the practical implementation of Algorithm \ref{alg:MXL}:

\paragraph{Updates and feasibility}
\label{par:feasible}

As we discussed earlier, the exponentiation/normalization step of Algorithm \ref{alg:MXL} ensures that the players' action variables $\mathbf{X}_{\play}(n)$ satisfy the game's feasibility constraints \eqref{eq:feasible} at each stage $n$.
In the noiseless case ($\bZ=0$), feasibility is automatic because, by construction, $\bV_{\play}$ (and, hence, $\mathbf{Y}_{\play}$) is Hermitian.
In the presence of noise however, there is no reason to assume that the estimates $\hat\bV_{\play}$ are Hermitian,
so the updates $\mathbf{X}_{\play}$ may also fail to be feasible.
To rectify this, we tacitly assume that
each player corrects such errors by replacing $\hat\bV_{\play}$ with $(\hat\bV_{\play} + \hat\bV_{\play}^{\dag})/2$.
Since this error-correcting operation is linear in its input, Hypotheses \eqref{eq:zeromean} and \eqref{eq:MSE} continue to apply, so our analysis and proofs hold as stated.


\paragraph{On the step-size sequence $\gamma_{n}$}
\label{par:stepsize}

Using a decreasing step-size $\gamma_{n}$ in \eqref{eq:MXL} may appear counter-intuitive because it implies that new information enters the algorithm with decreasing weights.
As evidenced by Theorem \ref{thm:convrate}, the reason for this is that
a constant step-size might cause the process to overshoot and lead to oscillatory behavior around the algorithm's end-state.
In the deterministic regime, these oscillations can be dampened by using forward-backward splitting or accelerated descent methods \cite{Nes04}.
However, in the presence of noise, the use of a decreasing step-size is essential in order to dissipate measurement noise and other stochastic effects, explaining why it is not possible to improve on \eqref{eq:convrate} by using a constant step.

\paragraph{Computational complexity and numerical stability}
\label{par:complexity}

From the point of view of computational complexity, the bottleneck of each iteration of Algorithm \ref{alg:MXL} is the matrix exponentiation step $\mathbf{Y}_{\play} \mapsto \exp(\mathbf{Y}_{\play})$.
Since matrix exponentiation has the same complexity as matrix multiplication, this step has polynomial complexity with a low degree on the input dimension of $\mathbf{X}_{\play}$ \cite{MvL03}.
In particular, each exponentiation requires $\bigoh(M_{\play}^{\omega})$ floating point operations, where the complexity exponent can be as low as $\omega = 2.373$ if the players  employ fast Coppersmith\textendash Winograd matrix multiplication methods \cite{DS13}. 

Finally, regarding the numerical stability of Algorithm \ref{alg:MXL}, the only possible source of arithmetic errors is the exponentiation/normalization step.
Indeed, if the eigenvalues of $\mathbf{Y}_{\play}$ are large, this step could incur an overflow where both the numerator and the denominator evaluate to machine infinity.
This potential instability can be fixed as follows:
If $y_{\play}$ denotes the largest eigenvalue of $\mathbf{Y}_{\play}$ and we let $\mathbf{Y}_{\play}' = \mathbf{Y}_{\play} - y_{\play} \mathbf{I}$, the algorithm's exponentiation step can be rewritten as:
\begin{equation}
\label{eq:exp-stable}
\mathbf{X}_{\play}
	\leftarrow \frac{\exp(\mathbf{Y}_{\play}')}{\exp(-y_{\play}) + \norm{\exp(\mathbf{Y}_{\play}')}}.
\end{equation}
Thanks to this shift, the elements of the numerator are now bounded from above by $1$ (because the largest diagonal element of $\mathbf{Y}_{\play}'$ is $e^{y_{\play} - y_{\play}} = 1$), so there is no danger of encountering a numerical indeterminacy of the type $\mathtt{Inf}/\mathtt{Inf}$.
Thus, to avoid computer arithmetic issues, we employ the stable expression \eqref{eq:exp-stable} in all numerical implementations of Algorithm \ref{alg:MXL}.

\paragraph{Asynchronous implementation}
\label{par:asynchronicity}

{
An implicit assumption in \eqref{eq:MXL} is that user updates \textendash\ albeit local \textendash\ are concurrent.
This can often be achieved via a global update timer that synchronizes the players' updates;
however, in a fully decentralized setting, even this degree of coordination may be challenging to attain.
Thus, to overcome this synchronicity limitation, we discuss below an asynchronous variant of Algorithm \ref{alg:MXL} where each player updates his action based on an individual, independent schedule.}


Specifically, assume that each player $\play\in\mathcal{\players}$ has an individual timer that triggers an update event, i.e. a request for gradient feedback and, subsequently, an update of $\mathbf{X}_{\play}$.
Of course, the players' gradient estimates $\hat \bV_{\play}$ will then be subject to delays and asynchronicities (in addition to noise), so the update structure of Algorithm \ref{alg:MXL} must be modified appropriately.
To that end, let $\playset_{n} \subseteq \playset$ be the set of players that update their actions at the $n$-th overall update event (typically $\abs{\playset_{n}}=1$ if players update at random times),
and let $d_{\play}(n)$ be the corresponding number of epochs that have elapsed since the last update of player $\play$.
We then obtain the following asynchronous variant of \eqref{eq:MXL}:
\begin{equation}
\label{eq:MXL-async}
\begin{aligned}
\mathbf{Y}_{\play}(n+1)
	&= \mathbf{Y}_{\play}(n) + \gamma_{n_{\play}} \one\braces{\play\in\playset_{n}} \cdot \hat\bV_{\play}(n),
	\\
\mathbf{X}_{\play}(n+1)
	&= A_{\play} \frac{\exp(\mathbf{Y}_{\play}(n+1))}{1 + \norm{\exp(\mathbf{Y}_{\play}(n+1))}},
\end{aligned}
\end{equation}
where
$n_{\play}$ denotes the number of updates performed by player $\play$ up to epoch $n$ while the (asynchronous) estimate $\hat\bV_{\play}(n)$ satisfies
\begin{equation}
\label{eq:grad-multi}
\exof{\hat\bV_{\play}(n)}
	= \bV_{\play}(\mathbf{X}_{1}(n - d_{1}(n)),\dotsc,\mathbf{X}_{N}(n - d_{N}(n))).
\end{equation}
Assuming that the delays $d_{\play}(n)$ are bounded and the players' update rates are strictly positive (in the sense that $\liminf_{n\to\infty} n_{\play}/n > 0$), it is possible to employ the analysis of \cite{CGM15} to show that Theorems \ref{thm:stationary}\textendash\ref{thm:local} remain true under the asynchronous variant \eqref{eq:MXL-async}.
However, a detailed analysis would take us too far afield, so we relegate it to future work.






\subsection{Contention-based medium access}

Contention-based \acf{MAC} aims to provide an efficient means for accessing and sharing a wireless channel in the presence of several competing wireless users that interfere with each other.
To model this, consider a set of wireless users indexed by $\playset = \{1,\dotsc,N\}$,
each updating their individual channel access probability $x_{\play}$ in response to the amount of contention in the network \cite{Chen2010}.
In practice, wireless users cannot be assumed to know the exact channel access probabilities of other users, so user $i$ infers the level of wireless contention via an aggregate \emph{contention measure} $q_{\play}(\mathbf{x}_{-\play})$ which is determined as a (symmetric) function of the access probability profile $\mathbf{x}_{-\play} = (x_{1},\dotsc,x_{i-1},x_{i+1}\dotsc,x_{N})$ of all other users.%
\footnote{For instance, a standard contention measure of this type is the conditional collision probability
$q_{\play}({\bf x}_{-\play}) = 1 - \prod_{\playalt\neq\play} (1-x_{\playalt})$. \cite{Chen2010}.}

With this in mind, the objective of each user is to select their individual channel access probability so as to maximize the benefit derived from acessing the channel more often minus the induced contention $x_{\play} q_{\play}(\mathbf{x}_{-\play})$ incurred by all other users.
This leads to the utility function formulation
\begin{equation}
u_{\play}(x_{\play};\mathbf{x}_{-\play})
	= U_{\play}(x_{\play}) - x_{\play} q_{\play}(\mathbf{x}_{-\play})
\end{equation}
where $U_{\play}$ is a continuous and nondecreasing function representing the utility of user $\play$ when there are no other users in the channel.
Thus, in economic terms, $u_{\play}$ simply represents the user's net gain from channel access, discounted by the associated contention cost.

We thus obtain the multi-agent, multi-objective formulation
\begin{flalign}
\label{eq:game-MAC}
\text{for all $\play\in\playset$}
	\quad
	\begin{cases}
	\textrm{maximize}
		&\quad
		u_{\play} (x_{\play};\mathbf{x}_{-\play}),
		\\
	\textrm{subject to}
		&\quad
		x_{\play} \in [0,1].
	\end{cases}
\end{flalign}
whose solutions can be analyzed through the specification of the utility function $U_{i}(x_{i})$ and the choice of the contention measure $q_{i}$. If $U_{i}(x_{i})$ is assumed to be continuously differentiable and strictly concave (see for example \cite{Cui08,Chen2010} and references therein), \eqref{eq:game-MAC} is a special case of \eqref{eq:game} with $M=1$ and $A_{\play}=1$.
We further note the unilateral optimization problems above are interconnected because of the dependence of the contention measure $q_{i}(\mathbf{x}_{-i})$ on the access probabilities of all users;
as a result, the theory of noncooperative games theory arises naturally as the most suitable framework to analyze and solve \eqref{eq:game-MAC} \cite{Cui08}.



\subsection{Metric learning for image similarity search}
\label{sec:image}


A key challenge in content-based image retrieval is to design an automatic procedure capable of retrieving documents from a large database based on their similarity to a given request, often distributed over several computing cores in a massively parallel computing grid (or cloud) \cite{Negrel-2013,Smeulders-2000}.
To formalize this, an image is typically represented via its \emph{signature}, i.e. a real $d$-dimensional vector $\mathbf{i}\in\mathbb{R}^{d}$ that collects and encodes the most distinctive features of said image.
Given a database $\mathcal{D}$ of such signatures, each image $\mathbf{i}\in\mathcal{D}$ is further associated with a set $\mathcal{S}_{\mathbf{i}}\subseteq\mathcal{D}$ of \emph{similar} images and a set $\mathcal{U}_{\mathbf{i}}\subseteq\mathcal{D}$ of \emph{dissimilar} ones, based on each image's content.
Accordingly, the goal is to design a distance metric which is minimized between similar images and is maximized between dissimilar ones.


A widely used distance measure between image signatures is the \emph{Mahalanobis distance}, defined as
\begin{equation}
\label{eq:Mahalanobis}
d_{\mathbf{X}} (\mathbf{i},\mathbf{j}) = (\mathbf{i}-\mathbf{j})^{\top}\, \mathbf{X} \, (\mathbf{i} - \mathbf{j}),
\end{equation}
where $\mathbf{X} \succcurlyeq 0$ is a $d\times d$ positive-definite matrix.%
\footnote{The baseline case $\mathbf{X} = \mathbf{I}$ corresponds to the Euclidean distance, which is often unsuitable for image discrimination purposes \cite{Bellet-2013}.}
This so-called ``precision matrix'' must then be chosen by the optimizer so that $d_{\mathbf{X}} (\mathbf{i},\mathbf{j}) < d_{\mathbf{X}} (\mathbf{i},\mathbf{k})$ whenever $\mathbf{i}$ is similar to $\mathbf{j}$ but dissimilar to $\mathbf{k}$.
With this in mind, we obtain the minimization objective
\begin{equation}
\label{eq:metric-all}
F(\mathbf{X};\mathcal{T})
	= \sum_{\mathclap{(\mathbf{i}, \mathbf{j},\mathbf{k}) \in\mathcal{T}}} C\big(d_{\mathbf{X}}(\mathbf{i},\mathbf{j}) -  d_{\mathbf{X}}(\mathbf{i},\mathbf{k}) - \varepsilon\big)
	+ \norm{\mathbf{X} - \mathbf{I}}_{F}^{2},
\end{equation}
where:
\begin{inparaenum}
[\itshape a\hspace*{1pt}\upshape)]
\item
$\mathcal{T}$ is the set of all triples $(\mathbf{i},\mathbf{j},\mathbf{k})$ such that $\mathbf{j}\in\mathcal{S}_{\mathbf{i}}$ and $\mathbf{k}\in\mathcal{U}_{\mathbf{i}}$;
\item
$C$ is a convex, nondecreasing function that penalizes matrices $\mathbf{X}$ that do not capture similarities and dissimilarities between images;
\item
the parameter $\varepsilon > 0$ reinforces this penalty;
and
\item
$\norm{\cdot}_{F}$ denotes the ordinary Frobenius ($\ell_{2}$) norm.%
\footnote{This last regularization term is included in order to maximize predictive accuracy by reducing the effects of over-fitting to training data \cite{Law-2014}.}
\end{inparaenum}

An additional requirement in the above is to employ a low-rank precision matrix $\mathbf{X}$ so as to
reduce model complexity and computational costs \cite{Negrel-2013},
enable distributed storing and retrieval \cite{Lim-2013},
and better exploit correlations between features that further reduce over-fitting effects \cite{Law-2014}.
A computationally efficient way to achieve this is to include a trace constraint of the form $\tr\mathbf{X} \leq c$ for some $c\ll d$,%
\footnote{By contrast, an $\ell_{0}$ rank constraint of the form $\rank(\mathbf{X}) \leq c$ generically leads to an untractable NP-hard problem formulation.}
leading to the feasible region
\begin{equation}
\label{eq:metric-feasible}
\boldsymbol{\mathcal{X}}
	= \setdef{\mathbf{X}\in\mathbb{R}^{d\times d}}{\mathbf{X} \succcurlyeq 0 \;\text{and}\; \tr\mathbf{X} \leq c}.
\end{equation}

Combining \eqref{eq:metric-all} and \eqref{eq:metric-feasible}, we see that metric learning in image retrieval is a special case of the general problem \eqref{eq:game} with $N=1$ optimizers.
However, when $\mathcal{D}$ is large, the number of variables involved is computationally prohibitive:
for instance, \eqref{eq:metric-all} may contain up to $10^{9}$\textendash$10^{11}$ terms for a modest database with $10^{4}$ images.
To circumvent this obstacle, a common approach is to replace $\mathcal{T}$ with a smaller, randomly drawn population sample $\mathcal{W}\subseteq\mathcal{T}$, and then take the average over all such samples.
In so doing, we obtain the stochastic optimization problem \cite{Bottou-2012}:
\begin{equation}
\label{eq:metric-stochastic-optim}
\begin{aligned}
\textrm{minimize}
	&\quad
	\exof{F(\mathbf{X}; \mathcal{W})}
	\hspace{1ex}

	\\
\textrm{subject to}
	&\quad
	\mathbf{X} \succcurlyeq 0, \tr(\mathbf{X}) \leq c,
\end{aligned}
\end{equation} 
where the expectation is taken over the random samples $\mathcal{W}$.

The benefit of this formulation is that, at each realization, only a small-size, tractable data sample $\mathcal{W}$ is used for calculations at each computing node.
On the flip side however, the expectation in \eqref{eq:metric-stochastic-optim} cannot be calculated, so the optimizer only has access to information on the random, realized gradients $\nabla_{\mathbf{X}} F(\mathbf{X};\mathcal{W})$.
This type of uncertainty is typical of stochastic optimization problems, and the proposed \ac{MXL} method has been designed precisely with this feedback structure in mind.



\subsection{Energy efficiency maximization}
\label{sec:EE}

Consider the problem of \acl{EE} maximization in multi-user, multiple-carrier \ac{MIMO} networks \textendash\ see e.g. \cite{BSHD15,ZSBJD16,MB16} and references therein.
Here, wireless connections are established over transceiver pairs with $M$ (resp. $N$) antennas at the transmitter (resp. receiver), and communication takes place over a set of $S$ orthogonal subcarriers.
Accordingly, the $\play$-th transmitter is assumed to control his individual input signal covariance matrix $\mathbf{Q}_{\plays}$ over each subcarrier $s=1,\dotsc,S$, subject to the constraints:
\begin{inparaenum}
[\itshape a\upshape)]
\item
$\mathbf{Q}_{\plays} \succcurlyeq 0$ (since each $\mathbf{Q}_{\plays}$ is a signal covariance matrix);
and
\item
$\tr\mathbf{Q}_{\play} \leq P_{\!\max}$,
where
$\mathbf{Q}_{\play} = \diag(\mathbf{Q}_{\plays})_{s=1}^{S}$ is the covariance profile of the $\play$-th transmitter,
$\tr(\mathbf{Q}_{\play})$ represents the user's transmit power over all subcarriers,
and
$P_{\!\max}$ denotes the user's \emph{maximum} transmit power.
\end{inparaenum}

\begin{figure}[t]
\centering
\input{Figures/NetworkFigure}
\caption{A multi-user \ac{MIMO} system with $\plays=3$ connections.}
\vspace{-2ex}
\end{figure}

Assuming Gaussian input and \ac{SUD} at the receiver, each user's Shannon-achievable throughput is given by the well-known expression
\begin{equation}
\label{eq:rate}
R_{\play}(\mathbf{Q})
	= \log\det\big( \mathbf{W}_{-\play} + \mathbf{H}_{\play\play}\mathbf{Q}_{\play}\mathbf{H}_{\play\play}^{\dag} \big)
	- \log\det\big( \mathbf{W}_{-\play} \big),
\end{equation}
where
$\mathbf{H}_{\playalt\play}$ denotes the channel matrix profile between the $\playalt$-th transmitter and the $\play$-th receiver over all subcarriers,
$\mathbf{Q} = (\mathbf{Q}_{1},\dotsc,\mathbf{Q}_{N})$ denotes the users' transmit profile and
$\mathbf{W}_{-\play}
	\equiv \mathbf{W}_{-\play}(\mathbf{Q})
	= \mathbf{I} + \sum\nolimits_{\playalt\neq\play} \mathbf{H}_{\playalt\play} \mathbf{Q}_{\playalt} \mathbf{H}_{\playalt\play}^{\dag}$
is the \ac{MUI} covariance matrix at receiver $\play$.
The users' transmit \acdef{EE} is then defined as the ratio between their Shannon rate and the total consumed power, i.e.
\begin{equation}
\label{eq:EE}
\ee_{\play}(\mathbf{Q})
	= \frac{R_{\play}(\mathbf{Q})}{P_{\!c} + \tr(\mathbf{Q}_{\play})},
\end{equation}
where $P_{\!c}>0$ represents the total power consumed by circuit components at the transmitter \cite{BSHD15}.

The energy efficiency function above is not concave w.r.t the covariance matrix $\mathbf{Q}_{\play}$ of user $i$,
but it can be recast as such via a suitable Charnes-Cooper transformation \cite{CC62}.
To this aim, following \cite{MB16}, consider the adjusted control variables
\begin{equation}
\label{eq:Xdef}
\mathbf{X}_{\play}
	= \frac{P_{\!c} + P_{\!\max}}{P_{\!\max}} \frac{\mathbf{Q}_{\play}}{P_{\!c} + \tr(\mathbf{Q}_{\play})},
\end{equation}
where the normalization constant $(P_{\!c} + P_{\!\max})/P_{\!\max}$ implies that $\tr(\mathbf{X}_{\play}) \leq 1$ with equality if and only if $\tr(\mathbf{Q}) = P_{\!\max}$.
The action set of user $\play$ is thus given by
\begin{equation}
\label{eq:new_strat}
\textstyle
\boldsymbol{\mathcal{X}}_{\play}
	= \setdef{\diag(\mathbf{X}_{\plays})_{s=1}^{S}}{\mathbf{X}_{\plays} \mgeq0 \text{ and } \sum\nolimits_{s} \tr \mathbf{X}_{\plays} \leq 1},
\end{equation}
and, using \eqref{eq:Xdef},
the \acl{EE} expression \eqref{eq:EE} yields the utility function
\begin{align}
\label{eq:new_EE}
u_{\play}(\mathbf{X}_{i};\mathbf{X}_{-i})
	= \frac{ P_{\!c} + (1 - \tr\mathbf{X}_{\play})\,P_{\!\max}}{P_{\!c}(P_{\!c} + P_{\!\max})}
	\log\det\left({
	\mathbf{I} + \frac{ P_{\!c} P_{\!\max}\, \wilde\bH_{\play} \mathbf{X}_{\play} \wilde\bH_{\play}^{\dag}}{P_{\!c} + (1 - \tr\mathbf{X}_{\play}) P_{\!\max}}
	}\right),
\end{align}
where $\wilde\bH_{\play}  = \cramped{\mathbf{W}_{-\play}^{-1/2}} \mathbf{H}_{\play \play}$ denotes the \emph{effective channel matrix} of user $\play$.

We thus see that the resulting \acl{EE} maximization problem is of the general form \eqref{eq:game} with $A_{\play} = 1$.%
\footnote{Strictly speaking, the block-diagonal constraint in \eqref{eq:new_strat} does not appear in \eqref{eq:feasible}, but it is satisfied automatically by the learning method presented in Sec.~\ref{sec:learning}.}
In this context, randomness and uncertainty stem from the noisy estimation of the users' \ac{MUI} covariance matrices (which depend on the other users' behavior),
the noise being due to the scarcity of perfect \ac{CSIT}, random measurement errors, etc.
To the best of our knowledge, the \ac{MXL} method discussed in Section \ref{sec:learning}
comprises the first distributed solution scheme for \acl{EE} maximization in general multi-user/multi-antenna/multi-carrier networks with local, causal \textendash\ and possibly imperfect \textendash\ \ac{CSI} feedback at the transmitter.

\subsection{Motivation and contributions}
\subsection{Related work}
\label{sec:related}

The \ac{MXL} method studied in this paper has strong ties to matrix regularization \cite{KSST12} and mirror descent methods \cite{Nes09,SS11} for (online) convex optimization.
In particular, the important special case of real vector variables on the simplex (real diagonal $\mathbf{X}_{\play}$ with constant trace) is closely related to the well-known \ac{EW} learning algorithm for multi-armed bandit problems \cite{Vov90}.
More recently, \ac{MXL} schemes were also proposed for single-user regret minimization in the context of online power control \cite{SMT15} and throughput maximization \cite{MB14} for dynamic \ac{MIMO} systems.
The goal there was to show that, in the presence of temporal variabilities, the long-run performance of \acl{MXL} matches that of the best fixed policy in hindsight, a property known as ``no regret'' \cite{SS11,CBL06}.
Nonetheless, in a multi-user, game-theoretic setting, a no-regret strategy may end up assigning positive weight \emph{only} to strictly dominated actions \cite{VZ13}.
As a result, (external) regret minimization is neither necessary nor sufficient to ensure convergence to a \acl{NE} (or other game-theoretic solution concept).

In \cite{MBM12,MM16,YN15}, a special case of \eqref{eq:game} was studied in the context of \ac{MIMO} throughput maximization in the presence of stochasticity, in both single-user \cite{YN15} and multiple access channels \cite{MM16,MBM12}.
In both cases, the problem boils down to a semidefinite optimization problem, possibly distributed over several agents \cite{MM16}.
The existence of a single objective function greatly facilitates the analysis;
however, many cases of practical interest (such as the examples we discuss in the following section) cannot be modeled as problems of this type, so a potential-based analysis is often unsuitable for practical applications.

In this paper, we derive the convergence properties of \acl{MXL} in concave $N$-player games, and we investigate the algorithm's long-term behavior under feedback errors and uncertainty.
Specifically, building on earlier work by Rosen \cite{Ros65}, we consider a broad class of games that satisfy a local \emph{variational stability} condition which ensures that \aclp{NE} are locally isolated.
In a series of recent papers, Scutari et al. used a stronger, global variant of this condition to establish the convergence of a class of Gauss\textendash Seidel, best-response methods,
and successfully applied these algorithms to a wide range of communication problems \textendash\ for a panoramic survey, see \cite{SFPP10} and references therein.
However, in the presence of noise and uncertainty, the convergence of best-response methods is often compromised:
as an example, in the case of throughput maximization with imperfect feedback, \acl{IWF} (a standard best-response scheme) \cite{YRBC04} fails to provide any perceptible gains over crude, uniform power allocation policies \cite{MM16}.
Thus, given that randomness, uncertainty and feedback imperfections are ubiquitous in practical systems, we focus throughout on attaining convergence results that are robust to learning impediments of this type.


\subsection{Main contributions and paper outline}
\label{sec:outline}
The main contribution of this work is the derivation of the convergence properties of \acl{MXL} for games played over bounded regions of positive-semidefinite matrices in the presence of feedback noise and randomness.
To put our theoretical analysis in context, we first discuss three examples from wireless networks and computer vision in Section \ref{sec:examples}.
More specifically, we illustrate
\begin{inparaenum}
[\itshape i\upshape)]
\item
a general game-theoretic framework for \emph{contention-based} \acdef{MAC};
\item
a stochastic optimization formulation for \emph{content-based image retrieval}
{(a key ``Big Data'' signal processing problem);}
and
\item
the multi-agent problem of transmit \acdef{EE} maximization in multi-carrier \ac{MIMO} networks (a critical design feature of emerging green multi-cellular networks).
\end{inparaenum}

In Section \ref{sec:model}, we revisit our core game-theoretic framework, and we discuss the notion of a \emph{stable} \acl{NE}.
Subsequently, in Section \ref{sec:learning}, we introduce our \acl{MXL} scheme, and we detail our assumptions for the stochasticity affecting the players' objectives and observations.
Our main results are then presented in Section \ref{sec:analysis} and can be summarized as follows:
Under fairly mild assumptions on the underlying stochasticity (zero-mean feedback errors with finite conditional variance), we show that
\begin{inparaenum}
[\upshape(\itshape i\hspace*{1pt}\upshape)]
\item
the algorithm's only termination states are \aclp{NE};
\item
if the game admits a globally (locally) stable equilibrium, then the algorithm converges globally (locally) to said equilibrium;
and
\item
on average, the algorithm converges to an $\varepsilon$-neighborhood of an extreme (resp. interior) strongly stable equilibrium within $\bigoh(1/\varepsilon)$ iterations (resp. $\bigoh(1/\varepsilon^{2})$).
\end{inparaenum}

The above results greatly extend and generalize the recent analysis of \cite{MM16} for throughput maximization in \ac{MIMO} \ac{MAC} systems.
To further validate our results in \ac{MIMO} environments, our theoretical analysis is supplemented with extensive numerical simulations for \acl{EE} maximization in multi-carrier \ac{MIMO} wireless networks in Section \ref{sec:numerics}.
To streamline the flow of the paper, technical proofs have been relegated to a series of appendices at the end.

\smallskip
\subsubsection*{Notation}
{
In what follows, the profile $\mathbf{X} = (\mathbf{X}_{1},\dotsc,\mathbf{X}_{N})$ is identified with the block-diagonal matrix $\diag(\mathbf{X}_{\play})_{\play=1}^{N}$, and we use the game-theoretic shorthand $(\mathbf{X}_{\play};\mathbf{X}_{-\play})$ when we need to focus on the action $\mathbf{X}_{\play}$ of player $\play$ against that of all other players.
Also, given $\mathbf{A}\in\mathbb{C}^{d\times d}$, we write $\norm{\mathbf{A}} = \sum_{j=1}^{d} \abs{\mathrm{eig}_{j}(\mathbf{A})}$ for the nuclear (trace) norm of $\mathbf{A}$ and $\dnorm{\mathbf{A}} = \max_{j} \abs{\mathrm{eig}_{j}(\mathbf{A})}$ for its (dual) singular norm.
}


\section{Introduction}
\label{sec:intro}
\input{Introduction}



\section{Motivation and Examples}
\label{sec:examples}
\input{Examples}



\section{Elements from Game Theory}
\label{sec:model}
\input{Model}



\section{Learning under Uncertainty}
\label{sec:learning}
\input{Learning}



\section{Convergence Analysis and Implementation}
\label{sec:analysis}
\input{Analysis}



\section{Numerical Results}
\label{sec:numerics}
\input{Numerics}



\section{Conclusions and Perspectives}
\label{sec:conclusions}
\input{Conclusions}



\appendices
\numberwithin{equation}{section}
\numberwithin{theorem}{section}
\numberwithin{proposition}{section}
\numberwithin{lemma}{section}


\section{The Exponentiation Step}
\label{app:entropy}
\input{App-Entropy}



\section{Stationarity Analysis}
\label{app:stationary}
\input{App-Stationary}



\section{Global Convergence}
\label{app:global}
\input{App-Global}



\section{Local Convergence}
\label{app:local}
\input{App-Local}



\section{Rates of Convergence}
\label{app:rates}
\input{App-Rates}




\bibliographystyle{IEEEtran}



