# Distributed Random Block-Coordinate descent methods for ill-posed composite
convex optimisation problems††thanks: Preliminary results reported in this
paper have been presented at the _61st IEEE Conference on Decision and Control
2022_ , [6].

Mathias Staudigl Maastricht University, Department of Advanced Computing
Sciences, P.O. Box 616, NL–6200 MD Maastricht, The Netherlands
<EMAIL_ADDRESS>Paulin Jacquot EDF R&D, OSIRIS, Palaiseau,
France
<EMAIL_ADDRESS>
###### Abstract

We develop a novel randomised block coordinate descent primal-dual algorithm
for a class of non-smooth ill-posed convex programs. Lying in the midway
between the celebrated Chambolle-Pock primal-dual algorithm and Tseng’s
accelerated proximal gradient method, we establish global convergence of the
last iterate as well optimal $O(1/k)$ and $O(1/k^{2})$ complexity rates in the
convex and strongly convex case, respectively, $k$ being the iteration count.
Motivated by distributed and data-driven control of power systems, we test the
performance of our method on a second-order cone relaxation of an AC-OPF
problem. Distributed control is achieved via the distributed locational
marginal prices (DLMPs), which are obtained dual variables in our optimisation
framework.

##  1 Introduction

In this paper we study non-smooth convex composite optimisation problems of
the form

$\min_{\bm{x}\in\mathcal{X}}\\{\Psi(\bm{x}):=\Phi(\bm{x})+R(\bm{x})\\},$ (P)

where
$\mathcal{X}:=\arg\min\left\\{\frac{1}{2}\|\bm{A}\bm{z}-\bm{b}\|^{2}|\bm{z}\in\mathbb{R}^{m}\right\\},$
and the functions $\Phi:\mathbb{R}^{m}\to\mathbb{R}$ and
$R:\mathbb{R}^{m}\to(-\infty,\infty]$ are convex and additively decomposable
of the form $\Phi(\bm{x}):=\sum_{i=1}^{d}\Phi_{i}(\bm{x}_{i})$ and
$R(\bm{x}):=\sum_{i=1}^{d}r_{i}(\bm{x}_{i})$. We assume that each function
$\phi_{i}:\mathbb{R}^{m_{i}}\to\mathbb{R}$ is smooth, whereas $r_{i}(\cdot)$
is only required to be proper convex and lower semi-continuous. We typically
think of the smooth component $\phi_{i}(\cdot)$ as a convex loss function,
whereas $r_{i}(\cdot)$ can take over to role of a regulariser or penalty
function. In particular, $r_{i}(\cdot)$ can represent an indicator function of
a closed convex set $\mathcal{K}_{i}\subset\mathbb{R}^{m_{i}}$, representing
individual membership constraints of the decision variable
$\bm{x}_{i}\in\mathbb{R}^{m_{i}}$. Thus, problem (P) can include certain
separable block constraints in addition to the non-separable constraints
embodied in the set $\mathcal{X}$. This setting gains relevance in distributed
optimisation problems with joint coupling constraints, such as multi-agent
optimisation problems [12]. Other examples include convex but non-smooth
penalty functions like the $\ell_{1}$-norm on $\mathbb{R}^{m_{i}}$ is also a
prominent choice in applications. The matrix
$\bm{A}=[\bm{A}_{1};\ldots;\bm{A}_{d}]$ is decomposed in $q\times m_{i}$
matrices $\bm{A}_{i}\in\mathbb{R}^{q\times m_{i}}$. Accordingly, we use the
notation $\bm{x}=[\bm{x}_{1};\ldots;\bm{x}_{d}]$ with each
$\bm{x}_{i}\in\mathbb{R}^{m_{i}}$ to represent the blocks of coordinates of
$\bm{x}$, and $m=m_{1}+\ldots+m_{d}$. The joint restriction
$\bm{x}\in\mathcal{X}$ state the a feasible decision variable is a solution to
a linear least squares problem and could be equivalently described as the set
of solutions to the normal equations
$\mathcal{X}=\\{\bm{x}\in\mathbb{R}^{m}|\bm{A}^{\top}\bm{A}\bm{x}=\bm{A}^{\top}\bm{b}\\}$.
When $\bm{b}$ is in the range of $\bm{A}$, then we can solve the linear system
$\bm{A}\bm{x}=\bm{b}$ exactly, and the optimisation problem reduces to a
linearly constrained convex non-smooth minimisation problem

$\min_{\bm{x}}\\{\Psi(\bm{x})=\Phi(\bm{x})+R(\bm{x})\\}\quad\text{s.t.:
}\bm{A}\bm{x}=\bm{b}.$ (1)

We call this the _consistent case_. Problems of type (1) are very general and
include all generic conic optimisation problems, such as linear programming,
second-order cone optimisation and semi-definite programming. Problems of this
kind are prominent in large-scale optimisation, with important applications in
inverse problems and signal processing, machine learning, and engineering. In
particular, partitioning the matrix $\bm{A}$ into suitably defined blocks,
problem (1) is a canonical formulation of various distributed optimisation
problems, as the next examples illustrate.

######  Example 1 (Consensus).

Consider the problem

$\min_{\bm{x}\in\mathbb{R}^{p}}\sum_{i=1}^{n}r_{i}(\bm{x})$

where $r_{1},\ldots,r_{d}$ are closed convex functions on $\mathbb{R}^{p}$.
This problem is equivalent to

$\min_{\bm{x}_{1},\ldots,\bm{x}_{d}}\sum_{i=1}^{n}r_{i}(\bm{x}_{i})\quad\text{
s.t.: }\bm{x}_{1}=\bm{x}_{2}=\ldots=\bm{x}_{n}$

This can be written as a linear constrained optimisation problem with matrix
$\bm{A}$ corresponding to the Laplacian of a connected undirected graph and
linear constraint $\bm{A}\bm{x}=0$.

######  Example 2 (Distributed model fitting).

Problem (P) can also cover composite minimisation problems which canonically
arise in machine learning. Consider the problem

$\min_{\bm{u}\in\mathbb{R}^{p}}\ell(\bm{K}\bm{u}-\bm{b})+\rho(\bm{u}),$

where $\ell(\cdot)$ is a convex and smooth statistical loss function of the
form

$\ell(\bm{K}\bm{u}-\bm{b})=\sum_{i=1}^{q}\ell_{i}(\bm{K}_{i}^{\top}\bm{u}-b_{i})$

where $\ell_{i}:\mathbb{R}\to\mathbb{R}$ is the loss for the $i$-th training
example, $\bm{K}_{i}\in\mathbb{R}^{p}$ is the feature vector for example $i$,
and $b_{i}$ is the output or response for example $i$. Define the variable
$\bm{x}=\left[\begin{array}[]{c}\bm{u}\\\
\bm{v}\end{array}\right]\in\mathbb{R}^{p+q}=\mathbb{R}^{m}$ and the linear
operator $\bm{A}\bm{x}=\bm{K}\bm{u}-\bm{v}$. Define the functions
$\phi_{i}:\mathbb{R}\to\mathbb{R}$ and $r_{i}:\mathbb{R}\to(-\infty,\infty]$
by

$\phi_{i}(\bm{x}_{i}):=\left\\{\begin{array}[]{ll}\ell_{i}(v_{i})&\text{if
}p+1\leq i\leq p+d\\\ 0&\text{ if }1\leq i\leq p\end{array}\right.\quad
r_{i}(\bm{x}_{i}):=\left\\{\begin{array}[]{ll}\rho_{i}(u_{i})&\text{if }1\leq
i\leq p\\\ 0&\text{ if }p+1\leq i\leq m\end{array}\right.$

Setting $\Phi(\bm{x})=\sum_{i=1}^{m}\phi_{i}(\bm{x}_{i})$ and
$R(\bm{x})=\sum_{i=1}^{m}r_{i}(\bm{x}_{i})$ yields the convex optimisation
problem

$\min_{\bm{x}}\Phi(\bm{x})+R(\bm{x})\quad\text{s.t.: }\bm{A}\bm{x}=\bm{b}.$

Problem (P) is more general as we do not insist on $\bm{b}$ belonging to the
range of $\bm{A}$, which means that the optimisation problems of interest to
us are non-smooth convex problems with ill-posed linear restrictions under
which the exact condition $\bm{A}\bm{x}=\bm{b}$ is replaced by the approximate
condition $\bm{A}\bm{x}\approx\bm{b}$. This assumption gains relevance in
inverse problems and PDE-constrained optimisation, where problems of the from
(P) appear in the method of residuals (Morozov regularisation) [21]. Another
instance where ill-posed linear restrictions appear is studied in an
application to power systems in Section 5 of this work. Motivated by all these
different optimisation problems, this paper derives a unified random block-
coordinate descent method which solves problem (P) under very general
assumptions on the problem data.

###  1.1 Related Methods

Our algorithms are closely related to randomised coordinate descent methods,
primal–dual coordinate update methods, and accelerated primal–dual methods. In
this subsection, let us briefly review the three classes of methods and
discuss their relations to our work.

#### Linear constrained minimisation

From the theoretical standpoint, one could formulate problem (P) as a linear
constrained optimisation problem of the type (1), with linear restriction
$\bm{A}^{\top}\bm{A}\bm{x}=\bm{A}^{\top}\bm{b}$ (the normal equations). Hence,
one could approach problem (P) via primal-dual techniques directly. While we
will show that our main algorithmic scheme (Algorithm 1) is indeed equivalent
to a suitably defined primal-dual method, it can be argued that this
connection is not always a recommended solution strategy in practice. First, a
primal-dual implementation has to deal with the symmetric matrix
$\bm{A}^{\top}\bm{A}$, whose dimension is $m\times m$. If $q$ is much smaller
than $m$ then primal-dual methods would have to process many more data points
than direct approaches. Second, if $\bm{A}$ is a sparse matrix, then
$\bm{A}^{\top}\bm{A}$ is usually not sparse anymore, which leads to heavy use
of numerical linear algebra techniques. Finally, it is often not known
a-priori whether the linear system is de facto consistent with the given
inputs.This is in particular the case in the application motivating the
development of the numerical scheme to be presented in this paper, which is
concerned with the distributed control of an electrical distribution system
within an AC-optimal power flow framework [46, 5]. A basic stability
desideratum on a numerical solution method (acting as a decentralised
coordination mechanism in our application) is therefore that overall system
convergence is guaranteed even if the linear constraints are not satisfied
with equality. Our method is exactly achieving this. Section 5 illustrates the
performance of our method on a 15-bus AC-OPF problem taken from [36].

#### Primal penalty methods

An alternative and popular approach to solve (P) is the penalty method. It
consists in solving a sequence of unconstrained optimisation problems
$\min_{\bm{x}}\Psi(\bm{x})+\frac{\sigma_{k}}{2}\|\bm{A}\bm{x}-\bm{b}\|^{2},$
where $\sigma_{k}$ is a positive and increasing penalty parameter sequence.
Intuitively it is clear that, by choosing $\sigma_{k}$ larger, the more
importance the optimisation gives to the constraints. Since penalty methods
are entirely primal, they do not use duality arguments, and hence they are in
principle able to solve the inconsistent case as well. However, their
implementation usually involves two loops: an inner loop solving the
minimisation problem for a fixed parameter $\sigma_{k}$ to some desired
accuracy, followed by an outer loop describing how the penalty parameter is
updated. Viewed from this angle, our algorithm is performing these operations
in a single-loop fashion.

#### Randomised coordinate descent methods

In the absence of the linear constraints, our algorithm specialises to
randomised coordinate descent (RCD), which was first proposed in [33], and
later generalised to the non-smooth case in [25, 40]. It was shown that RCD
features sublinear rates of convergence with rate $O(1/k)$, $k$ being the
iteration counter. Acceleration to a $O(1/k^{2})$ complexity and even linear
rates for strongly convex problems have been obtained. Extensions to parallel
computations, important for large-scale optimization problems, were first
proposed in [41].

#### Primal-dual coordinate update methods

To cope with linear constraints, a very popular approach is the alternating
direction method of multipliers (ADMM). Originally, ADMM [20, 16] was proposed
for two-block structured problems with separable objective. The convergence
and complexity analysis of this method is well documented in the literature
[22]; see [11] for a survey. Direct extensions of the ADMM to multi-block
settings such as (P) is not straightforward, and indeed even may fail to
converge [9]. Very recently, [18] proposed a randomised primal–dual coordinate
(RPDC) update method, whose asynchronous parallel version was then studied in
[19]. It was shown that RPDC converges with rate $O(1/k)$ under convexity
assumption.

#### Accelerated primal-dual methods

It is possible to accelerate the rate of convergence from $O(1/k)$ to
$O(1/k^{2})$ for gradient type methods. The first acceleration result was
shown by [32] for solving smooth unconstrained problems. The technique has
been generalised to accelerate gradient-type methods on possibly non-smooth
convex programs [3, 34]. Primal–dual methods on solving linearly constrained
problems can also be accelerated by similar techniques. Under convexity
assumption, the augmented Lagrangian method (ALM) is accelerated from $O(1/k)$
to $O(1/k^{2})$ in [23].

###  1.2 Contribution

#### Methodological contributions

We propose a block-coordinate implementation of the method developed by [28]
for linearly constrained optimisation, lying midway between the celebrated
Chambolle-Pock primal-dual splitting algorithm [8] and Tseng’s accelerated
proximal gradient method [45]. Specifically, our proposed method is a
distributed interpretation of the primal-dual algorithm of [8] which operates
on randomly selected coordinate blocks. A parallel between the Chambolle-Pock
method [8] and the accelerated proximal gradient of [45] was already drawn in
[27, 28]. Reducing the primal-dual algorithm to an implementation of Tseng’s
method enabled [27] to derive new convergence results based on primal
arguments, thus departing from strong duality requirements and the ergodic
rates typically issued for primal-dual methods. Our developments revisit the
coordinate-descent implementation proposed in [26] for the basic algorithm of
[28], to a block-coordinate descent setting featuring Nesterov type of
acceleration. In particular, in the strongly convex case, we derive a new step
size policy that achieves an accelerated rate of $O(k^{-2})$. Besides the
recent preprint [43], we are not aware of a similar algorithm achieving the
accelerated convergence rates $O(k^{-2})$ in a fully distributed computational
setting. Thus, our main contributions in this paper can be summarised as
follows:

  * (i)

In the convex case, provided problem (P) possesses a saddle point (defined in
Section 2.5) our main result is Theorem 1 which establishes an $O(1/k)$
iteration complexity in terms of the objective function gap and convergence of
the last iterate.

  * (ii)

In the strongly convex case and uniform sampling of coordinate blocks, our
main result is Theorem 2, which proves an accelerated $O(k^{-2})$ convergence
rate in the primal objective function values.

We remark that RCD methods have been shown to exhibit linear convergence rates
in the strongly convex regime [33]. Such fast convergence is, however, not to
be expected in the presence of linear constraints. While strong convexity of
the primal objective ensures smoothness of the Lagrangian dual function, but
not its strong concavity. Hence, in general, we do not expect to see linear
convergence rates by only assuming strong convexity in the primal. However, we
note that [47] obtain linear convergence rates in the well-posed regime if
there is one block variable that is independent of all others in the objective
(but coupled in the linear constraint) and also the corresponding component
function is smooth.
Related to this paper is also the very recent work [44]. They consider a
larger class of convex optimisation problems with linear constraints and
design a new randomized primal-dual algorithm with single block activation in
each iteration and similar complexity results as reported in the present work.
However, our method is able to solve ill-posed convex optimisation problems
with general sampling techniques.

###  1.3 Organisation of this paper

This paper is organised as follows. Section 2 describes our block coordinate
descent framework. Section 3 explains in detail our algorithmic approach.
Section 4 contains all details for the asymptotic convergence and complexity
statements. Section 5 describes a challenging application of our algorithm to
a distributed optimisation formulation of an AC-OPF problem formulate the
distribution grid model, based on the second-order cone relaxation of [13,
37]. Preliminary numerical results are reported to show the applicability of
our method using a 15 bus network studied in [36] as a concrete example.

###  1.4 Notation

We work in the space $\mathbb{R}^{p}$ composed by column vectors. For
$\bm{x},\bm{u}\in\mathbb{R}^{p}$ denote the standard Euclidean inner product
$\langle{\bm{x}},{\bm{u}}\rangle=\bm{x}^{\top}\bm{u}$ and the Euclidean norm
$\|\bm{x}\|=\langle{\bm{x}},{\bm{x}}\rangle^{1/2}$. We let
$\mathbb{S}^{p}_{+}:=\\{\bm{B}\in\mathbb{R}^{p\times
p}|\bm{B}^{\top}=\bm{B},\bm{B}\succeq 0\\}$ the space of positive definite
matrices. Given $\Lambda\in\mathbb{S}^{p}_{+}$, we let
$\|\bm{x}\|_{\Lambda}:=\langle{\Lambda\bm{x}},{\bm{x}}\rangle^{1/2}$ for
$\bm{x}\in\mathbb{R}^{p}$. The identity matrix of dimension $p$ is denoted as
$\bm{I}_{p}$. Whenever the dimension is clear from the context, we will just
write $\bm{I}$. We denote by $\lambda_{\max}(\bm{A})$ the largest eigenvalue
of a square $p\times p$ matrix $\bm{A}$. We call by
$\Gamma_{0}(\mathbb{R}^{p})$ the set of proper convex, lower semi-continuous
functions $f:\mathbb{R}^{p}\to(-\infty,+\infty]$. For such a function
$f\in\Gamma_{0}(\mathbb{R}^{p})$ the effective domain is defined as
$\mathrm{dom}(f):=\\{\bm{x}|f(\bm{x})<\infty\\}$. For $d\in\mathcal{N}$, we
set $[d]:=\\{1,\ldots,d\\}$. For a vector $\nu\in\mathbb{R}^{p}_{++}$, we let
$\nu^{-1}$ the vector of reciprocal values. For
$\Gamma\in\mathbb{S}^{m}_{++}$, define the weighted proximal operator
$\operatorname{\texttt{prox}}^{\Gamma}_{r}(\bm{x}):=\arg\min_{\bm{u}\in\mathbb{R}^{m}}\\{r(\bm{u})+\frac{1}{2}\|\bm{u}-\bm{x}\|_{\Gamma}^{2}\\}.$
If $\Gamma:=\operatorname{\texttt{blkdiag}}(\Gamma_{1},\ldots,\Gamma_{d})$,
then the proximal operator decomposes accordingly
$\operatorname{\texttt{prox}}^{\Gamma}_{r}(\bm{x})=\left(\operatorname{\texttt{prox}}^{\Gamma_{1}}_{r_{1}}(\bm{x}_{1}),\ldots,\operatorname{\texttt{prox}}_{r_{d}}^{\Gamma_{d}}(\bm{x}_{d})\right).$
If $f:\mathbb{R}^{m}\to\mathbb{R}$ is differentiable, we denote the gradient
of $f$ at $\bm{x}\in\mathbb{R}^{m}$ by $\nabla f(\bm{x})\in\mathbb{R}^{m}$.

##  2 Preliminaries

###  2.1 Block structure

We first describe the block setup which has become standard in the analysis of
block coordinate descent methods [40, 15, 41]. The block structure of (P) is
given by a decomposition of $\mathbb{R}^{m}$ into $d$ subspaces
$\mathbb{R}^{m_{i}},1\leq i\leq d$, so that
$\mathbb{R}^{m}=\mathbb{R}^{m_{1}}\times\cdots\times\mathbb{R}^{m_{d}}$. Let
$\bm{U}=[\bm{U}_{1};\ldots;\bm{U}_{d}]$ be the $m\times m$ identity matrix,
decomposed into column submatrices $\bm{U}_{i}\in\mathbb{R}^{m\times m_{i}}$.
For $\bm{x}\in\mathbb{R}^{m}$, let $\bm{x}_{i}=\bm{U}^{\top}_{i}\bm{x}$ be the
block of coordinates corresponding to the columns of $\bm{U}_{i}$. Any vector
$\bm{s}\in\mathbb{R}^{m}$ can be written as
$\bm{s}=\sum_{i=1}^{d}\bm{U}_{i}\bm{s}_{i}$. For $\emptyset\neq
I\subseteq[d]$, we write

$\bm{s}_{[I]}:=\sum_{i\in I}\bm{U}_{i}\bm{s}_{i}.$

We denote the $\ell^{2}$-norm on $\mathbb{R}^{m_{i}}$ as $\|\cdot\|_{i}$. If
$\bm{Q}=\operatorname{\texttt{blkdiag}}[\bm{Q}_{1};\ldots,\bm{Q}_{d}]$ is a
block-diagonal matrix with $\bm{Q}_{i}\in\mathbb{S}^{m_{i}}_{++}$, we define a
weighted norm on $\mathbb{R}^{m}$ by

$\|\bm{s}\|_{\bm{Q}}:=\sum_{i=1}^{d}\|\bm{s}_{i}\|_{\bm{Q}_{i}}\qquad\forall\bm{s}\in\mathbb{R}^{m}.$

###  2.2 Smoothness of $\Phi$

We assume throughout the paper that $\Phi:\mathbb{R}^{m}\to\mathbb{R}$ is
strongly and convex and possesses a Lipschitz continuous partial gradient.
Specifically, we assume that for each block $i\in[d]$ there exists a matrix
$\Lambda_{i}\in\mathbb{S}^{m_{i}}_{++}$ so that

$0\leq\phi_{i}(\bm{x}_{i}+\bm{t}_{i})-\phi_{i}(\bm{x}_{i})-\langle{\nabla\phi_{i}(\bm{x})},{\bm{t}_{i}}\rangle\leq\frac{1}{2}\|\bm{t}_{i}\|^{2}_{\Lambda_{i}},$
(2)

A typical situation is that $\Lambda_{i}=L_{i}\bm{I}_{m_{i}}$ for some
$L_{i}>0$, so that the gradient $\nabla\phi_{i}$ is Lipschitz continuous with
modulus $L_{i}>0$. Collecting all the matrices $\Lambda_{i}$ in one block
diagonal matrix
$\Lambda:=\operatorname{\texttt{blkdiag}}[\Lambda_{1};\ldots;\Lambda_{d}]$, it
follows

$\Phi(\bm{x}^{\prime})-\Phi(\bm{x})-\langle{\nabla\Phi(\bm{x})},{\bm{x}^{\prime}-\bm{x}}\rangle\leq\frac{1}{2}\|\bm{x}^{\prime}-\bm{x}\|^{2}_{\Lambda}.$

###  2.3 Properties of $R$

We assume that $R:\mathbb{R}^{m}\to(-\infty,\infty]$ and is block separable
$R(\bm{x})=\sum_{i=1}^{d}r_{i}(\bm{x}_{i})$, where the functions
$r_{i}\in\Gamma_{0}(\mathbb{R}^{m_{i}})$ are $\mu_{i}$-strongly convex and
closed, with $\mu_{i}\geq 0$. i.e. the shifted function $\bm{x}_{i}\mapsto
r_{i}(\bm{x}_{i})-\frac{\mu_{i}\|\bm{x}_{i}\|^{2}_{i}}{2}$ is convex. Calling
$\Upsilon_{i}:=\mu_{i}\bm{I}_{m_{i}}$, this gives

$r(\bm{x}^{\prime}_{i})\geq
r_{i}(\bm{x}_{i})+\langle{\xi_{i}},{\bm{x}^{\prime}_{i}-\bm{x}_{i}}\rangle+\frac{1}{2}\|\bm{x}^{\prime}_{i}-\bm{x}_{i}\|^{2}_{\Upsilon_{i}}\qquad\forall\bm{x}^{\prime}_{i}\in\mathbb{R}^{m_{i}},\xi\in\partial
r_{i}(\bm{x}_{i}).$ (3)

Typical examples for the regulariser $r_{i}$ are indicator functions of closed
convex sets, i.e. $r_{i}(\bm{x}_{i})=\delta_{\mathcal{K}_{i}}(\bm{x}_{i})$,
for $\mathcal{K}_{i}\subset\mathbb{R}^{m_{i}}$ convex and closed, or
structure-imposing regularisers like the $L_{p}$-penalty
$r_{i}(\bm{x}_{i})=\|\bm{x}_{i}\|^{p}_{L_{p}(\mathbb{R}^{m_{i}})}$ ($p\geq
1$), prominent in distributed estimation of high-dimensional signals and
neural networks. We let
$\Upsilon:=\operatorname{\texttt{blkdiag}}[\Upsilon_{1};\ldots;\Upsilon_{d}]$
be the $m\times m$ matrix collecting all strong-convexity parameters of the
individual blocks.

###  2.4 Properties of the quadratic penalty function

Define the function

$h(\bm{x}):=\frac{1}{2}\|\bm{A}\bm{x}-b\|^{2}.$ (4)

Let $h^{*}=\min_{\bm{x}\in\mathbb{R}^{d}}h(\bm{x})$, so that
$\mathcal{X}=\arg\min_{\mathbb{R}^{d}}h(\bm{x})$. Clearly, $h^{\ast}\geq 0$,
with equality if and only if the linear system $\bm{A}\bm{x}=b$ is consistent.
Since $h$ is quadratic, we have for all
$\bm{u},\bm{w},\bm{x}\in\mathbb{R}^{m}$

$h(\bm{u})-h(\bm{x})=\langle{\nabla
h(\bm{w})},{\bm{u}-\bm{x}}\rangle+\frac{1}{2}\|\bm{A}(\bm{u}-\bm{w})\|^{2}-\frac{1}{2}\|\bm{A}(\bm{x}-\bm{w})\|^{2}.$
(5)

In particular, if $\bm{x}^{*}\in\mathcal{X}$, the above implies for
$\bm{x}=\bm{w}=\bm{x}^{*}$

$h(\bm{x})-h(\bm{x}^{*})=\frac{1}{2}\|\bm{A}(\bm{x}-\bm{x}^{*})\|^{2}.$

We also have

$\displaystyle h(\bm{x}+\bm{U}_{i}\bm{t}_{i})$
$\displaystyle=h(\bm{x})+\langle{\nabla
h(\bm{x})},{\bm{U}_{i}\bm{t}_{i}}\rangle+\frac{1}{2}\|\bm{A}_{i}\bm{t}_{i}\|^{2}$
$\displaystyle=h(\bm{x})+\langle{\bm{U}_{i}^{\top}\nabla
h(\bm{x})},{\bm{t}_{i}}\rangle+\frac{1}{2}\|\bm{t}_{i}\|^{2}_{\bm{A}_{i}^{\top}\bm{A}_{i}}$
$\displaystyle\leq h(\bm{x})+\langle{\bm{U}_{i}^{\top}\nabla
h(\bm{x})},{\bm{t}_{i}}\rangle+\frac{\lambda_{i}}{2}\|\bm{t}_{i}\|^{2}_{i}$

where
$\lambda_{i}\equiv\lambda_{\max}(\bm{A}_{i}^{\top}\bm{A}_{i})\equiv\|\bm{A}_{i}\|_{2}$,
the spectral norm of the matrix $\bm{A}_{i}$, and Lemma 9 (proven in Appendix
A) immediately implies that for all $t\in[0,1]$

$h(t\bm{x}+(1-t)\bm{x}^{\prime})=th(\bm{x})+(1-t)h(\bm{x}^{\prime})-\frac{t(1-t)}{2}\|\bm{A}(\bm{x}-\bm{x}^{\prime})\|^{2}.$
(6)

###  2.5 On Saddle points

Problem (P) can be equivalently expressed as the linearly constrained
optimisation problem

$\min_{\bm{x}=(\bm{x}_{1},\ldots,\bm{x}_{d})}\Psi(\bm{x})=\Phi(\bm{x})+R(\bm{x})\quad\text{s.t.:
}\bm{A}^{\top}\bm{A}\bm{x}=\bm{A}^{\top}b$ (7)

The Lagrangian associated with this non-smooth convex optimization problem is

$\mathcal{L}(\bm{x},\bm{y})=\Psi(\bm{x})+\langle{\bm{y}},{\bm{A}^{\top}b-\bm{A}^{\top}\bm{A}\bm{x}}\rangle.$

######  Definition 1.

A pair $(\bm{x}^{\ast},\bm{y}^{\ast})$ is called a _saddle-point_ if

$0\in\partial\Psi(\bm{x}^{\ast})-\bm{A}^{\top}\bm{A}\bm{y}^{\ast},\text{ and
}\bm{A}^{\top}\bm{A}\bm{x}^{\ast}-\bm{A}^{\top}b=0.$ (8)

For convex programs, the conditions (8) are sufficient for $\bm{x}^{\ast}$ to
be an optimal solution of (P). They are also necessary if a constraint
qualification condition holds (e.g. the Slater condition, stating that there
exists $\bm{x}$ in the interior of the domain of $\Psi$ such that
$\bm{A}^{\top}\bm{A}\bm{x}=\bm{A}^{\top}b$).

###  2.6 Random Sampling

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space. A _sampling_ is
a random element with values in $2^{[d]}$. We will call the random variable
$\mathcal{I}:\Omega\to 2^{[d]}$ a _random sampling_. $\mathcal{I}(\omega)$
defines the selection of blocks employed in a single iteration of our method.
The set $\mathcal{I}(\Omega)=\Sigma$ is the set of all possible realisations
of the random selection mechanism. Let $\\{\iota_{k}\\}_{k\in\mathbb{N}}$
represent the stochastic process on $2^{[d]}$ in which each random variable
$\iota_{k}$ is an i.i.d copy of $\mathcal{I}$. We refer to such a sampling as
an _i.i.d sampling_. We assume that the sampling $\mathcal{I}$ is _proper_ :
There exists a vector $\bm{\pi}=[\pi_{1},\ldots,\pi_{d}]\in\mathbb{R}^{d}$
with

$\pi_{i}=\mathbb{P}(i\in\mathcal{I})\in(0,1)\qquad\forall i\in[d].$

With the sampling $\mathcal{I}$, we associate the matrix
$\bm{\Pi}\in\mathbb{R}^{d\times d}$ defined as

$(\bm{\Pi})_{ij}:=\mathbb{P}(\\{i,j\\}\subseteq\mathcal{I})\qquad\forall i\neq
j,\text{ and }(\bm{\Pi})_{ii}=\pi_{i}.$

We note that $\bm{P}\succ 0$ [39, Thm 3.1]. We further define the weighting
matrix

$\bm{P}:=\operatorname{\texttt{blkdiag}}[\frac{1}{\pi_{1}}\bm{I}_{m_{1}};\ldots;\frac{1}{\pi_{d}}\bm{I}_{m_{d}}].$

We emphasise that the random sampling model we adopt here is capable of
capturing many stationary randomised activation mechanisms. To illustrate
this, consider the following activation mechanisms (see [38, 39] for further
discussions):

  * •

_Single coordinate activation:_ at each iteration, one coordinate block is
activated. This means that $\mathcal{I}(\omega)$ takes values in the discrete
set $[d]$ only, i.e. $\Sigma=[d]$. In this case, we necessarily have
$\sum_{i=1}^{d}\pi_{i}=1$.

  * •

_Uniform Sampling:_ for all $i\in[d]$ it holds
$\mathbb{P}(i\in\mathcal{I})=\mathbb{P}(j\in\mathcal{I})$. This implies
$\pi_{i}=\frac{\mathbb{E}[\lvert\mathcal{I}\rvert]}{d}$ for all $i\in[d]$. A
special case of a uniform sampling is the popular class of $m$-nice samplings,
arising under the specification where $\Sigma$ is the set of all subsets of
$[d]$ containing exactly $m$ elements, each of which is activated with uniform
probability. Clearly, in this case one has $\pi_{i}=\frac{m}{d}$ for all
$i\in[d]$.

  * •

_Full Sampling:_ $\Sigma=\\{1,\ldots,d\\}$, which means that all coordinates
are updated in parallel.

##  3 Parallel distributed coordinate-descent algorithm

Parameters :
$(\pi_{i})_{i=1}^{d},\bm{Q}^{k}_{i}=\frac{1}{\pi_{i}\tau^{k}}\bm{B}_{i},(\sigma_{k})_{k\geq
0},(\tau_{k})_{k\geq 0}$

Initialization : $\bm{x}^{0}=\bm{w}^{0}\in\mathbb{R}^{m},S_{-1}=0,S_{0}=1$

for  _$k=0,1,2,\dots$_ do

1

2$\bm{z}^{k}=\frac{1}{S_{k}}(S_{k-1}\bm{w}^{k}+\sigma_{k}\bm{x}^{k})$

3 draw $\iota_{k}\subseteq[d]$ as an i.i.d copy of the sampling $\mathcal{I}$.

4 Update
$\bm{x}_{i}^{k+1}=\left\\{\begin{array}[]{ll}\operatorname{\mathsf{T}}^{k}_{i}(\bm{x}^{k})\text{
defined in eq.}\eqref{eq:FB}&\text{if }i\in\iota_{k},\\\
\bm{x}_{i}^{k}&\text{otherwise.}\end{array}\right.$
$\bm{w}^{k+1}=\bm{z}^{k}+\frac{\sigma_{k}}{S_{k}}\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})$

5 $S_{k+1}=S_{k}+\sigma_{k+1}$

6

7

Algorithm 1 Distributed Accelerated Proximal Gradient Algorithm

Our random block-coordinate descent algorithm for solving (P) recursively
updates three sequences $\\{(\bm{z}^{k},\bm{w}^{k},\bm{x}^{k})\\}_{k\geq 0}$.
Let $(\sigma_{k})_{k\geq 0}$ be a given sequence of positive numbers. At each
iteration $k\geq 0$, we define a weight sequence $(S_{k})_{k\geq 0}$
recursively by

$S_{0}=1,\quad
S_{k}=S_{k-1}+\sigma_{k},\quad\theta_{k}=\frac{\sigma_{k}}{S_{k}}.$

Consider the extrapolated point

$\bm{z}^{k}=(1-\theta_{k})\bm{w}^{k}+\theta_{k}\bm{x}^{k}=\bm{w}^{k}+\theta_{k}(\bm{x}^{k}-\bm{w}^{k}),$
(9)

together with the sequence

$\bm{w}^{k+1}=\bm{z}^{k}+\theta_{k}\bm{P}(\bm{x}^{k+1}-\bm{x}^{k}).$ (10)

This reads in coordinates
$\bm{w}_{i}^{k+1}=\bm{z}_{i}^{k}+\frac{\theta_{k}}{\pi_{i}}(\bm{x}_{i}^{k+1}-\bm{x}_{i}^{k})$
for all $i\in[d]$. To evaluate $\bm{w}^{k+1}$, we need the primal update
$\bm{x}^{k+1}$ which is obtained by a weighted forward-backward step involving
the first-order signal
$g^{k}_{i}=\nabla\phi_{i}(\bm{x}^{k}_{i})+S_{k}\nabla_{i}h(\bm{z}^{k}).$
Specifically, given a block-specific scaling matrix
$\bm{B}_{i}\in\mathbb{S}^{m_{i}}_{++}$, we define
$\bm{Q}_{i}^{k}:=\frac{1}{\pi_{i}\tau_{k}}\bm{B}_{i}$, and the _weighted
forward-backward operator_ by
$\operatorname{\mathsf{T}}^{k}_{i}(\bm{x}^{k})=\operatorname{\texttt{prox}}_{r_{i}}^{\bm{Q}^{k}_{i}}(\bm{x}^{k}_{i}-(\bm{Q}^{k}_{i})^{-1}g^{k}_{i})$,
which reads explicitly as

$\operatorname{\mathsf{T}}^{k}_{i}(\bm{x}^{k})=\arg\min_{\bm{u}_{i}\in\mathbb{R}^{m_{i}}}\left\\{r_{i}(\bm{u}_{i})+\langle{g^{k}_{i}},{\bm{u}_{i}-\bm{x}_{i}}\rangle+\frac{1}{2}\|\bm{u}_{i}-\bm{x}_{i}\|^{2}_{\bm{Q}^{k}_{i}}\right\\}.$
(11)

We will choose the matrices $\bm{B}_{i}$ later on to adapt for the strong
convexity present in the problem data. To obtain a compact representation of
the primal updates, let us introduce the iid Bernoulli process
$\epsilon_{i}^{k}:\Omega\to\\{0,1\\}$ by

$\epsilon_{i}^{k}(\omega)=\left\\{\begin{array}[]{cc}1&\text{if
}i\in\iota_{k}(\omega)\\\ 0&\text{else}\end{array}\right.$

and the random matrix
$\bm{E}_{k}=\operatorname{\texttt{blkdiag}}[\epsilon_{1}^{k}\bm{I}_{m_{1}},\ldots,\epsilon^{k}_{d}\bm{I}_{m_{d}}]\in\\{0,1\\}^{m\times
m}$. This matrix corresponds to the identity for those blocks that are
activated in round $k$, and zero for those not activated. Therefore, we call
it the _activation matrix_ of the Algorithm. By definition of the sampling
process, we have
$\mathbb{E}(\epsilon^{k}_{i})=\mathbb{P}(i\in\iota_{k})=\pi_{i}$ for all
$i\in[d]$ and $k\geq 0$, so that $\mathbb{E}[\bm{E}_{k}]=\bm{P}^{-1}$.

###  3.1 Connections among the iterates

Let
$\epsilon^{k}:=(\epsilon_{1}^{k},\ldots,\epsilon_{d}^{k})\in\\{0,1\\}^{d}$,
and $\mathcal{F}_{k}:=\sigma(\bm{x}^{0},\epsilon^{0},\ldots,\epsilon^{k-1})$
denote the history of the process up to step $k$. We denote by
$\mathbb{E}_{k}[\cdot]:=\mathbb{E}[\cdot|\mathcal{F}_{k}]$ the conditional
expectation at stage $k$. Denoting by
$\hat{\bm{x}}^{k+1}=\operatorname{\mathsf{T}}^{k}(\bm{x}^{k})$, we obtain the
compact representation of the primal update as

$\bm{x}^{k+1}=\bm{x}^{k}+\bm{E}_{k}(\hat{\bm{x}}^{k+1}-\bm{x}^{k}).$ (12)

Furthermore,
$\mathbb{E}_{k}(\bm{x}^{k+1})=\bm{P}^{-1}\hat{\bm{x}}^{k+1}+(\bm{I}-\bm{P}^{-1})\bm{x}^{k},$
so that

$\mathbb{E}_{k}[\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})]=\hat{\bm{x}}^{k+1}.$
(13)

######  Lemma 1.

For all $k\geq 0$, we have

$\begin{split}\frac{1}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})\|^{2}&=\frac{1}{2}\mathbb{E}_{k}\left[\|\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\bm{w}^{k})\|^{2}\right]\\\
&-\frac{1}{2}\mathbb{E}_{k}\left[\|\bm{A}(\bm{P}\bm{E}_{k}-\bm{I})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\|^{2}\right].\end{split}$
(14)

###### Proof.

Note that for any random variable $X:\Omega\to\mathbb{R}^{q}$, it holds true
$\|\mathbb{E}[X]\|^{2}=\mathbb{E}\left[\|X\|^{2}\right]-\mathbb{E}\left[\|X-\mathbb{E}X\|^{2}\right]$.
Set $X:=\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\bm{w}^{k})$. Then,
by eq. (13), we see $\mathbb{E}_{k}[X]=\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})$.
Furthermore,

$X-\mathbb{E}_{k}[X]=\bm{A}\left((\bm{I}-\bm{P})\bm{x}^{k}+\bm{P}\bm{x}^{k+1}-\hat{\bm{x}}^{k+1}\right)=\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\hat{\bm{x}}^{k+1}).$
(15)

Using (12), we get

$\bm{P}\bm{x}^{k+1}-\hat{\bm{x}}^{k+1}=(\bm{P}\bm{E}_{k}-\bm{I})\hat{\bm{x}}^{k+1}+\bm{P}(\bm{I}-\bm{E}_{k})\bm{x}^{k},$

which implies

$\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\hat{\bm{x}}^{k+1})=X-\mathbb{E}_{k}[X]=\bm{A}(\bm{P}\bm{E}_{k}-\bm{I})(\hat{\bm{x}}^{k+1}-\bm{x}^{k}).$
(16)

Collecting all therms above, one obtains

$\displaystyle\frac{1}{2}\mathbb{E}_{k}\left[\|\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\bm{w}^{k})\|^{2}\right]$
$\displaystyle=\frac{1}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})\|^{2}$
$\displaystyle+\frac{1}{2}\mathbb{E}_{k}\left[\|\bm{A}(\bm{P}\bm{E}_{k}-\bm{I})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\|^{2}\right].$

∎

The next Lemma provides a compact expression for the correction term in eq.
(14).

######  Lemma 2.

Let
$\bm{\Xi}=\mathbb{E}[\bm{E}_{k}\bm{P}\bm{A}^{\top}\bm{A}\bm{P}\bm{E}_{k}]$,
which is a symmetric time-independent matrix, thanks to iid sampling. We have

$\mathbb{E}_{k}\left[\|\bm{A}(\bm{P}\bm{E}_{k}-\bm{I})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\|^{2}\right]=\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\Xi-\bm{A}^{\top}\bm{A}}.$
(17)

###### Proof.

By definition of the sampling procedure, we have
$\mathbb{E}_{k}[\bm{P}\bm{E}_{k}]=\bm{I}_{m}$. Some simple algebra shows then
that

$\displaystyle\mathbb{E}_{k}\left[\|\bm{A}(\bm{P}\bm{E}_{k}-\bm{I})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\|^{2}\right]$
$\displaystyle=\mathbb{E}_{k}\left[(\hat{\bm{x}}^{k+1}-\bm{x}^{k})^{\top}\bm{M}_{k}(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\right],$

where
$\bm{M}_{k}:=\bm{E}_{k}\bm{P}\bm{A}^{\top}\bm{A}\bm{P}\bm{E}_{k}-\bm{A}^{\top}\bm{A}\bm{P}\bm{E}_{k}-\bm{E}_{k}\bm{P}\bm{A}^{\top}\bm{A}+\bm{A}^{\top}\bm{A}$
is a symmetric matrix, satisfying
$\mathbb{E}_{k}\bm{M}_{k}=\bm{\Xi}-\bm{A}^{\top}\bm{A}$. Hence,

$\displaystyle\mathbb{E}_{k}\left[(\hat{\bm{x}}^{k+1}-\bm{x}^{k})^{\top}\bm{M}_{k}(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\right]$
$\displaystyle=\mathbb{E}_{k}\left[\operatorname{Tr}\left(\bm{M}_{k}(\hat{\bm{x}}^{k+1}-\bm{x}^{k})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})^{\top}\right)\right]$
$\displaystyle=\operatorname{Tr}\left(\mathbb{E}_{k}[\bm{M}_{k}](\hat{\bm{x}}^{k+1}-\bm{x}^{k})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})^{\top}\right)$
$\displaystyle=\operatorname{Tr}\left((\Xi-\bm{A}^{\top}\bm{A})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})^{\top}\right)$
$\displaystyle=\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|_{\bm{\Xi}-\bm{A}^{\top}\bm{A}}.$

∎

######  Remark 1.

The matrix $\bm{\Xi}$ can be given a simple expression in terms of the
probability matrix $\bm{\Pi}$. A direct computation shows that $\bm{\Xi}$ is a
block-symmetric matrix with $d^{2}$ blocks $\bm{\Xi}[i,j],1\leq i,j\leq d,$
each block having the form

$\bm{\Xi}[i,j]=\frac{(\bm{\Pi})_{ij}}{\pi_{i}\pi_{j}}\bm{A}^{\top}_{i}\bm{A}_{j}^{\top}\text{
if }i\neq j,\text{ and
}\bm{\Xi}[i,i]=\frac{1}{\pi_{i}}\bm{A}^{\top}_{i}\bm{A}_{i}\quad\forall
i\in[d].$

In special instances this matrix can be simplified even further:

  1. 1.

If $\mathbb{P}(\lvert\mathcal{I}\rvert=1)=1$ the random sampling is a single
coordinate activation mechanism. In this case the associated matrix $\bm{\Pi}$
is diagonal with entries $\pi_{i}$. It follows that
$\bm{\Xi}[i,i]=\frac{1}{\pi_{i}}\bm{A}_{i}^{\top}\bm{A}_{i}$ for all $i\in[d]$
and $\bm{\Xi}[i,j]=0$ for $i\neq j$.

  2. 2.

We say that the matrix $\bm{A}$ has _orthogonal design_ if
$\bm{A}_{i}^{\top}\bm{A}_{j}$ is the zero matrix for $i\neq j$. Also in this
case the matrix $\bm{\Xi}$ is block diagonal, with the same coordinate
expression as above.

###  3.2 Properties of the penalty function $h$

We collect some important identities involving the penalty function $h$ in
this subsection. From the definition of the iterate $\bm{z}^{k}$, we see

$\begin{split}\nabla
h(\bm{z}^{k})&=\bm{A}^{\top}(\bm{A}\bm{z}^{k}-\bm{b})=\frac{1}{S_{k}}\left(S_{k-1}\bm{A}^{\top}(\bm{A}\bm{w}^{k}-\bm{b})+\sigma_{k}\bm{A}^{\top}(\bm{A}\bm{x}^{k}-\bm{b})\right)\\\
&=(1-\theta_{k})\nabla h(\bm{w}^{k})+\theta_{k}\nabla
h(\bm{x}^{k}).\end{split}$ (18)

Furthermore, we note that the definition of the iterate $\bm{w}^{k+1}$ implies

$\begin{split}h(\bm{w}^{k+1})&=h\left((1-\theta_{k})\bm{w}^{k}+\theta_{k}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k}))\right)\\\
&\overset{\eqref{eq:hstrongconvex}}{=}(1-\theta_{k})h(\bm{w}^{k})+\theta_{k}h(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k}))\\\
&-\frac{\theta_{k}(1-\theta_{k})}{2}\|\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\bm{w}^{k})\|^{2}.\end{split}$
(19)

######  Lemma 3.

We have

$\mathbb{E}_{k}[h(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k}))]=h(\hat{\bm{x}}^{k+1})+\frac{1}{2}\mathbb{E}_{k}[\|\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k}))-\hat{\bm{x}}^{k+1}\|^{2}].$
(20)

###### Proof.

The Pythagorean identity gives

$\displaystyle h(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k}))$
$\displaystyle=\frac{1}{2}\|(\bm{A}\hat{\bm{x}}^{k+1}-\bm{b})+\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\hat{\bm{x}}^{k+1})\|^{2}$
$\displaystyle=h(\hat{\bm{x}}^{k+1})+\frac{1}{2}\|\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\hat{\bm{x}}^{k+1})\|^{2}$
$\displaystyle+\langle{\bm{A}\hat{\bm{x}}^{k+1}-\bm{b}},{\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\hat{\bm{x}}^{k+1})}\rangle.$

Taking conditional expectations on both sides, using eq. (13), establishes the
claimed identity. ∎

###  3.3 Relation to primal-dual methods

Our method is very similar to the recent block coordinate primal-dual update
[47], who focus on the consistent case and uniform samplings [33, 41, 39, 15].
We generalise this to arbitrary samplings and show that the sequences produced
by Algorithm 1 are equivalent to a primal-dual process in the spirit of [26],
formulated as Algorithm 2.

Parameters :
$(\pi_{i})_{i=1}^{d},\bm{Q}^{k}_{i}=\frac{1}{\pi_{i}\tau^{k}}\bm{B}_{i},(\sigma_{k})_{k\in\mathbb{N}_{0}}$

Initialization :
$\bm{x}^{0}\in\mathbb{R}^{m},\bm{y}^{0}=\sigma_{0}(\bm{A}\bm{x}^{0}-\bm{b}),\bm{u}^{0}=\bm{A}\bm{x}^{0}-\bm{b}$

for  _$k=0,1,2,\dots$_ do

1 draw block $\iota_{k}$ as an i.i.d copy of the sampling $\mathcal{I}$.

2 Update

43
$\bm{x}_{i}^{k+1}=\left\\{\begin{array}[]{ll}\operatorname{\texttt{prox}}^{\bm{Q}^{k}_{i}}_{r_{i}}\left(\bm{x}^{k}_{i}-(\bm{Q}^{k}_{i})^{-1}(\nabla\phi_{i}(\bm{x}^{k}_{i})+\bm{A}^{\top}_{i}\bm{y}^{k})\right)&\text{
if }i\in\iota_{k},\\\ \bm{x}_{i}^{k}&\text{else.}\end{array}\right.$
$\bm{u}^{k+1}=\bm{u}^{k}+\bm{A}(\bm{x}^{k+1}-\bm{x}^{k})\;$

5$\bm{y}^{k+1}=\bm{y}^{k}+\sigma_{k}\bm{A}\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})+\sigma_{k+1}\bm{u}^{k+1}$

6

7

Algorithm 2 Primal-Dual Block Coordinate Descent Algorithm

Let $\\{(\bm{x}^{k},\bm{w}^{k},\bm{z}^{k})\\}_{k\geq 0}$ denote the sequences
generated by running Algorithm 1. Let $S_{k}$ be the cumulative step-size
process $S_{k}=1+\sum_{t=1}^{k}\sigma_{t}.$ Introduce the sequence
$\bm{y}^{k}:=S_{k}(\bm{A}\bm{z}^{k}-\bm{b})$, so that

$S_{k}\nabla_{i}h(\bm{z}^{k})=\bm{A}_{i}^{\top}S_{k}(\bm{A}\bm{z}^{k}-\bm{b})=\bm{A}_{i}^{\top}\bm{y}^{k}.$

In terms of this new dual variable $\bm{y}^{k}$, we can reorganize the primal
update so that

$\displaystyle\hat{\bm{x}}^{k+1}_{i}$
$\displaystyle=\arg\min_{\bm{x}_{i}\in\mathbb{R}^{m_{i}}}\left\\{r_{i}(\bm{x}_{i})+\langle{\nabla\phi_{i}(\bm{x}^{k}_{i})+S_{k}\nabla_{i}h(\bm{z}^{k})},{\bm{x}_{i}-\bm{x}^{k}_{i}}\rangle+\frac{1}{2}\|\bm{x}_{i}-\bm{x}_{i}\|^{2}_{\bm{Q}^{k}_{i}}\right\\}.$
$\displaystyle=\arg\min_{\bm{x}_{i}\in\mathbb{R}^{m_{i}}}\left\\{r_{i}(\bm{x}_{i})+\langle{\nabla\phi_{i}(\bm{x}^{k}_{i})+\bm{A}_{i}^{\top}\bm{y}^{k}},{\bm{x}_{i}-\bm{x}^{k}_{i}}\rangle+\frac{1}{2}\|\bm{x}_{i}-\bm{x}^{k}_{i}\|^{2}_{\bm{Q}^{k}_{i}}\right\\}$
$\displaystyle=\operatorname{\texttt{prox}}^{\bm{Q}_{i}^{k}}_{r_{i}}\left(\bm{x}^{k}_{i}-(\bm{Q}^{k}_{i})^{-1}(\nabla\phi_{i}(\bm{x}^{k}_{i})+\bm{A}^{\top}_{i}\bm{y}^{k})\right).$

The next iterate $\bm{x}^{k+1}$ is obtained by the block-coordinate update
rule (12). This gives line 3 of Algorithm 2.

Next, observe that

$\begin{split}S_{k}(\bm{A}\bm{w}^{k+1}-\bm{b})&=S_{k}(\bm{A}(\bm{z}^{k}+\theta_{k}\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\bm{b}))\\\
&=\bm{y}^{k}+\theta_{k}S_{k}\bm{A}\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})\\\
&=\bm{y}^{k}+\sigma_{k}\bm{A}\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})\end{split}$ (21)

Hence, after introducing the residual $\bm{u}^{k}=\bm{A}\bm{x}^{k}-\bm{b}$,
satisfying

$\bm{u}^{k+1}-\bm{u}^{k}=\bm{A}(\bm{x}^{k+1}-\bm{x}^{k})=\sum_{i\in\iota_{k}}\bm{A}_{i}(\hat{\bm{x}}^{k+1}_{i}-\bm{x}^{k}_{i}),$
(22)

we obtain line 4 of Algorithm 2, as well as

$\displaystyle\bm{y}^{k+1}$
$\displaystyle=S_{k+1}\left(\bm{A}(\bm{w}^{k+1}+\theta_{k+1}(\bm{x}^{k+1}-\bm{w}^{k+1}))-\bm{b}\right)$
$\displaystyle=S_{k+1}(\bm{A}\bm{w}^{k+1}-\bm{b})+\sigma_{k+1}(\bm{A}(\bm{x}^{k+1}-\bm{w}^{k+1})$
$\displaystyle=S_{k}(\bm{A}\bm{w}^{k+1}-\bm{b})+\sigma_{k+1}(\bm{A}\bm{x}^{k+1}-\bm{b})$
$\displaystyle=\bm{y}^{k}+\sigma_{k}\bm{A}\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})+\sigma_{k+1}\bm{u}^{k+1},$

where we have used in the last step the identity (21). Thereby, we obtain line
5 of Algorithm 2. This completes the verification that the sequence
$\\{(\bm{x}^{k},\bm{w}^{k},\bm{z}^{k})\\}_{k\geq 0}$ generated by Algorithm 1
is equivalent to the just constructed sequence
$\\{(\bm{x}^{k},\bm{u}^{k},\bm{y}^{k})\\}_{k\geq 0}$ corresponding to the
iterates of Algorithm 2.

######  Remark 2.

The implementation of Algorithm 2 is fully parallelisable, involving a
computational architecture with $d$ agents and a single central coordinator.
The agents manage the coordinate blocks $\bm{x}_{i}$ in a fully decentralised
way, using information about the centrally updated dual variable $\bm{y}^{k}$
only. A practical implementation of the computational scheme is as follows:

  1. 1.

Given the current data $(\bm{x}^{k},\bm{u}^{k},\bm{y}^{k})$ the coordinator
realises a sampling $\iota_{k}$.

  2. 2.

All agents in $\iota_{k}$ receive the order to update their control variables
in parallel, given their current position $\bm{x}_{i}^{k}$, the data matrix
$\bm{A}_{i}$ and the penalty of resource utilisation $\bm{y}^{k}$.

  3. 3.

Once all active agents executed their computation, they report to the vector
$\bm{A}_{i}(\bm{x}^{k+1}_{i}-\bm{x}^{k}_{i})$ to the central coordinator.

  4. 4.

The central coordinator updates the penalty parameter $\bm{y}^{k}$ by
executing the dual update in line 5 of the Algorithm 2

Distributed primal-dual methods as the one described have received enormous
attention in the control and machine learning community; see e.g [7, 31, 4,
47, 14, 24].

######  Remark 3.

Consider the special case with $\pi_{i}=1/d$ and $d=1$, as well as
$\sigma_{k}\equiv\sigma$. Then, Algorithm 2 coincides with the primal-dual
method of Chambolle-Pock [8]. In fact, in this case it follows that
$\bm{u}^{k}=\bm{A}\bm{x}^{k}-\bm{b}$, and
$\bm{y}^{k+1}=\bm{y}^{k}+\sigma(\bm{A}(2\bm{x}^{k+1}-\bm{x}^{k})-b)$.

##  4 Convergence Analysis

This section is concerned with the convergence properties of Algorithm 1. We
start with a basic descent property of the primal forward-backward step. This
will involve the identification of a Lyapunov function to obtain energy decay
estimates in a variable metric. Building on this result, we investigate two
important scenarios in isolation. First, we consider the merely convex case,
which is obtained when $\Upsilon=0$. If $\Upsilon\succ 0$, then accelerated
rates in the primal iterates can be obtained. This, however, requires a
delicate analysis of a step-size policy, which is relegated to Appendix B.

###  4.1 Lyapunov function and key estimates

We start the analysis with a small extension of ”Property 1” stated in [45]
for Bregman proximal gradient algorithms.

######  Lemma 4.

For all $k\geq 0$ and $\bm{x}\in\mathbb{R}^{m}$, define

$\zeta^{k}(\bm{x})=\Phi(\bm{x}^{k})+\langle{\nabla\Phi(\bm{x}^{k})},{\bm{x}-\bm{x}^{k}}\rangle+S_{k}\langle{\nabla
h(\bm{z}^{k})},{\bm{x}-\bm{z}^{k}}\rangle+\frac{1}{2}\|\bm{x}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}}$
(23)

Then, for all $\bm{x}\in\mathbb{R}^{m}$ it holds true that

$R(\hat{\bm{x}}^{k+1})+\zeta^{k}(\hat{\bm{x}}^{k+1})\leq
R(\bm{x})+\zeta^{k}(\bm{x})-\frac{1}{2}\|\bm{x}-\hat{\bm{x}}^{k+1}\|^{2}_{\bm{Q}^{k}+\Upsilon},$
(24)

###### Proof.

Collect the forward-backward operator in coordinates
$\operatorname{\mathsf{T}}^{k}(\bm{x}):=[\operatorname{\mathsf{T}}^{k}_{1}(\bm{x}_{1});\ldots;\operatorname{\mathsf{T}}^{k}_{d}(\bm{x}_{d})]$,
and set $\hat{\bm{x}}^{k+1}=\operatorname{\mathsf{T}}^{k}(\bm{x}^{k})$. From
the definition of the forward-backward operator (11) we see that

$0\in\partial
r_{i}(\hat{\bm{x}}^{k+1})+\nabla_{i}\zeta_{k}(\hat{\bm{x}}^{k+1})\qquad\forall
i\in[d].$

Hence, for all $\bm{x}\in\mathbb{R}^{m}$ we get from eq. (3)

$r_{i}(\bm{x}_{i})\geq
r_{i}(\hat{\bm{x}}^{k+1}_{i})+\langle{-\nabla_{i}\zeta^{k}(\hat{\bm{x}}^{k+1})},{\bm{x}_{i}-\hat{\bm{x}}^{k+1}_{i}}\rangle+\frac{1}{2}\|\bm{x}_{i}-\hat{\bm{x}}^{k+1}_{i}\|^{2}_{\Upsilon_{i}},$

for all $i\in[d]$. Summing over all blocks $i\in[d]$, it follows

$R(\bm{x})\geq
R(\hat{\bm{x}}^{k+1})+\langle{-\nabla\zeta^{k}(\hat{\bm{x}}^{k+1})},{\bm{x}-\hat{\bm{x}}^{k+1}}\rangle+\frac{1}{2}\|\bm{x}-\hat{\bm{x}}^{k+1}\|^{2}_{\Upsilon},$

Furthermore, it is easy to see that $\bm{x}\mapsto\zeta^{k}(\bm{x})$ is
1-strongly convex in the norm $\|\cdot\|_{\bm{Q}^{k}}$. Hence,

$\zeta^{k}(\bm{x})\geq\zeta^{k}(\hat{\bm{x}}^{k+1})+\langle{\nabla\zeta^{k}(\hat{\bm{x}}^{k+1})},{\bm{x}-\hat{\bm{x}}^{k+1}}\rangle+\frac{1}{2}\|\bm{x}-\hat{\bm{x}}^{k+1}\|^{2}_{\bm{Q}^{k}}.$

Adding these two inequalities and rearranging gives the claimed result. ∎

Lemma 1, together with Lemma 2 shows that

$\displaystyle\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{w}^{k}\|^{2}_{\bm{A}^{\top}\bm{A}}$
$\displaystyle=\frac{1}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})\|^{2}$
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:Ex1}}}{{=}}\frac{1}{2}\mathbb{E}_{k}\left[\|\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\bm{w}^{k})\|^{2}\right]$
$\displaystyle-\frac{1}{2}\mathbb{E}_{k}\left[\|\bm{A}(\bm{P}\bm{E}_{k}-\bm{I})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\|^{2}\right]$
$\displaystyle=\frac{1}{2}\mathbb{E}_{k}\left[\|\bm{A}(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})-\bm{w}^{k})\|^{2}\right]-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|_{\Xi-\bm{A}^{\top}\bm{A}}$

Applying eq. (15), (16), (19) and (20), as well as the identity
$\theta_{k}=\frac{\sigma_{k}}{S_{k}}$, the above becomes

$\begin{split}\frac{1}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})\|^{2}&\overset{\eqref{eq:hwk}}{=}\frac{S_{k}}{\sigma_{k}}h(\bm{w}^{k})+\frac{S_{k}}{S_{k-1}}\mathbb{E}_{k}[h(\bm{x}^{k}+\bm{P}(\bm{x}^{k+1}-\bm{x}^{k}))]\\\
&-\frac{S^{2}_{k}}{\sigma_{k}S_{k-1}}\mathbb{E}_{k}[h(\bm{w}^{k+1})]-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\Xi-\bm{A}^{\top}\bm{A}}\\\
&\overset{\eqref{eq:hforward},\eqref{eq:centerd2}}{=}\frac{S_{k}}{\sigma_{k}}h(\bm{w}^{k})+\frac{S_{k}}{2S_{k-1}}\mathbb{E}_{k}[\|\bm{A}(\bm{P}\bm{E}_{k}-\bm{I})(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\|^{2}]\\\
&-\frac{S^{2}_{k}}{\sigma_{k}S_{k-1}}\mathbb{E}_{k}[h(\bm{w}^{k+1})]-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\Xi-\bm{A}^{\top}\bm{A}}+\frac{S_{k}}{S_{k-1}}h(\hat{\bm{x}}^{k+1})\\\
&\overset{\eqref{eq:Ex2}}{=}\frac{S_{k}}{\sigma_{k}}h(\bm{w}^{k})+\frac{\sigma_{k}}{2S_{k-1}}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\Xi-\bm{A}^{\top}\bm{A}}-\frac{S^{2}_{k}}{\sigma_{k}S_{k-1}}\mathbb{E}_{k}[h(\bm{w}^{k+1})]\\\
&+\frac{S_{k}}{S_{k-1}}h(\hat{\bm{x}}^{k+1}).\end{split}$ (25)

Finally, we need a characterisation of the sequence $(\bm{w}^{k})_{k\geq 0}$,
which is quite standard in the analysis of random coordinate descent methods
[15, 38, 39]. We therefore skip the straightforward proof.

######  Lemma 5.

Let $(\bm{x}^{k},\bm{w}^{k})_{k\geq 0}$ be the iterates of Algorithm 1. Then,
for all $k\geq 1$, we have

$\bm{w}^{k}_{i}=\sum_{t=0}^{k}\gamma_{i}^{k,t}\bm{x}^{t}_{i},$ (26)

where for each $i\in[d]$, the coefficients $(\gamma_{i}^{k,t})_{t=0}^{k}$ are
defined recursively by setting
$\gamma_{i}^{0,0}=1,\gamma_{i}^{1,0}=1-\frac{\theta_{0}}{\pi_{i}},\gamma_{i}^{1,1}=\frac{\theta_{0}}{\pi_{i}}$,
and for all $k\geq 1$

$\gamma^{k+1,t}_{i}:=\left\\{\begin{array}[]{ll}(1-\theta_{k})\gamma_{i}^{k,t}&\text{
if }t=0,1,\ldots,k-1,\\\
(1-\theta_{k})\gamma_{i}^{k,k}+\theta_{k}(1-\pi_{i}^{-1})&\text{ if }t=k,\\\
\theta_{k}/\pi_{i}&\text{if }t=k+1.\end{array}\right.$ (27)

Moreover, for all $k\geq 0$, the following identity holds

$\gamma_{i}^{k+1,k}+\gamma^{k+1,k+1}_{i}=\theta_{k}+(1-\theta_{k})\gamma_{i}^{k,k}.$
(28)

Moreover, if $\theta_{0}\in(0,\min_{i\in[d]}\pi_{i}]$ and $(\theta_{k})_{k\geq
0}$ is a decreasing sequence, then for all $k\geq 0$ and $i\in[d]$, the
coefficients $(\gamma_{i}^{k,t})_{t=0}^{k}$ are all positive and add up to 1.

Let $\hat{\bm{x}}^{k+1}=\operatorname{\mathsf{T}}_{k}(\bm{x}^{k})$. Using (2),
we see

$\displaystyle\Psi(\hat{\bm{x}}^{k+1})$
$\displaystyle=\Phi(\hat{\bm{x}}^{k+1})+R(\hat{\bm{x}}^{k+1})$
$\displaystyle\leq
R(\hat{\bm{x}}^{k+1})+\Phi(\bm{x}^{k})+\langle{\nabla\Phi(\bm{x}^{k})},{\hat{\bm{x}}^{k+1}-\bm{x}^{k}}\rangle+\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\Lambda}$

Using (23), we see

$\displaystyle\Phi(\bm{x}^{k})+\langle{\nabla\Phi(\bm{x}^{k})},{\hat{\bm{x}}^{k+1}-\bm{x}^{k}}\rangle=\zeta^{k}(\hat{\bm{x}}^{k+1})-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}}-S_{k}\langle{\nabla
h(\bm{z}^{k})},{\hat{\bm{x}}^{k+1}-\bm{z}^{k}}\rangle.$

Consequently, for $\bm{x}^{\ast}\in\mathcal{X}=\arg\min_{\bm{x}}h(\bm{x})$ as
a reference point, the following bounds are obtained:

$\displaystyle\Psi(\hat{\bm{x}}^{k+1})$ $\displaystyle\leq
R(\hat{\bm{x}}^{k+1})+\zeta^{k}(\hat{\bm{x}}^{k+1})-S_{k}\langle{\nabla
h(\bm{z}^{k})},{\hat{\bm{x}}^{k+1}-\bm{z}^{k}}\rangle-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda}$
$\displaystyle\overset{\eqref{eq:ProxInequality}}{\leq}R(\bm{x}^{*})+\zeta^{k}(\bm{x}^{*})-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}}-S_{k}\langle{\nabla
h(\bm{z}^{k})},{\hat{\bm{x}}^{k+1}-\bm{z}^{k}}\rangle$
$\displaystyle-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda}$
$\displaystyle=\Psi(\bm{x}^{*})-\Phi(\bm{x}^{*})+\zeta^{k}(\bm{x}^{*})-S_{k}\langle{\nabla
h(\bm{z}^{k})},{\hat{\bm{x}}^{k+1}-\bm{z}^{k}}\rangle$
$\displaystyle-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda}$
$\displaystyle\overset{\eqref{eq:zeta}}{=}\Psi(\bm{x}^{*})-\Phi(\bm{x}^{*})+\Phi(\bm{x}^{k})+\langle{\nabla\Phi(\bm{x}^{k})},{\bm{x}^{*}-\bm{x}^{k}}\rangle+S_{k}\langle{\nabla
h(\bm{z}^{k})},{\bm{x}^{*}-\hat{\bm{x}}^{k+1}}\rangle$
$\displaystyle+\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}+\Upsilon}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda}$
$\displaystyle\overset{\eqref{eq:phiconvex}}{=}\Psi(\bm{x}^{*})+S_{k}\langle{\nabla
h(\bm{z}^{k})},{\bm{x}^{*}-\hat{\bm{x}}^{k+1}}\rangle+\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}+\Upsilon}$
$\displaystyle-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda}.$

Via eq. (18), we obtain

$S_{k}\langle{\nabla
h(\bm{z}^{k})},{\bm{x}^{*}-\hat{\bm{x}}^{k+1}}\rangle=S_{k-1}\langle{\nabla
h(\bm{w}^{k})},{\bm{x}^{*}-\hat{\bm{x}}^{k+1}}\rangle+\sigma_{k}\langle{\nabla
h(\bm{x}^{k})},{\bm{x}^{*}-\hat{\bm{x}}^{k+1}}\rangle.$

Next, we apply identity (5) to each inner product separately, in order to
conclude

$\displaystyle\langle{\nabla
h(\bm{w}^{k})},{\bm{x}^{*}-\hat{\bm{x}}^{k+1}}\rangle$
$\displaystyle=h(\bm{x}^{\ast})-h(\hat{\bm{x}}^{k+1})-\frac{1}{2}\|\bm{A}(\bm{x}^{*}-\bm{w}^{k})\|^{2}+\frac{1}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})\|^{2},$
$\displaystyle\langle{\nabla
h(\bm{x}^{k})},{\bm{x}^{*}-\hat{\bm{x}}^{k+1}}\rangle$
$\displaystyle=h(\bm{x}^{\ast})-h(\hat{\bm{x}}^{k+1})-\frac{1}{2}\|\bm{A}(\bm{x}^{*}-\bm{x}^{k})\|^{2}+\frac{1}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\|^{2}.$

Combined with the previous display, this shows for
$\bm{x}^{\ast}\in\arg\min_{\bm{x}\in\mathbb{R}^{m}}h(\bm{x})$,

$\displaystyle\Psi(\hat{\bm{x}}^{k+1})-\Psi(\bm{x}^{*})$ $\displaystyle\leq-
S_{k}[h(\hat{\bm{x}}^{k+1})-h(\bm{x}^{*})]-\frac{S_{k-1}}{2}\|\bm{A}(\bm{x}^{*}-\bm{w}^{k})\|^{2}+\frac{S_{k-1}}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})\|^{2}$
$\displaystyle-\frac{\sigma_{k}}{2}\|\bm{A}(\bm{x}^{*}-\bm{x}^{k})\|^{2}+\frac{\sigma_{k}}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{x}^{k})\|^{2}$
$\displaystyle+\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}+\Upsilon}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda}$
$\displaystyle=-\left(S_{k}h(\hat{\bm{x}}^{k+1})+S_{k-1}h(\bm{w}^{k})+\sigma_{k}h(\bm{x}^{k})-2S_{k}h(\bm{x}^{*})\right)+\frac{S_{k-1}}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})\|^{2}$
$\displaystyle+\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}+\Upsilon}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda-\sigma_{k}\bm{A}^{\top}\bm{A}}.$

Finally, applying identity (25), one sees

$\displaystyle-\left(S_{k}h(\hat{\bm{x}}^{k+1})+S_{k-1}h(\bm{w}^{k})+\sigma_{k}h(\bm{x}^{k})-2S_{k}h(\bm{x}^{*})\right)+\frac{S_{k-1}}{2}\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{w}^{k})\|^{2}$
$\displaystyle=-\left(S_{k}h(\hat{\bm{x}}^{k+1})+S_{k-1}h(\bm{w}^{k})+\sigma_{k}h(\bm{x}^{k})-2S_{k}h(\bm{x}^{*})\right)$
$\displaystyle+S_{k-1}\left(\frac{S_{k}}{\sigma_{k}}h(\bm{w}^{k})+\frac{\sigma_{k}}{2S_{k-1}}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\Xi-\bm{A}^{\top}\bm{A}}-\frac{S^{2}_{k}}{\sigma_{k}S_{k-1}}\mathbb{E}_{k}[h(\bm{w}^{k+1})]+\frac{S_{k}}{S_{k-1}}h(\hat{\bm{x}}^{k+1})\right)$
$\displaystyle=\frac{S^{2}_{k-1}}{\sigma_{k}}h(\bm{w}^{k})-\frac{S^{2}_{k}}{\sigma_{k}}\mathbb{E}_{k}[h(\bm{w}^{k+1})]+2S_{k}h(\bm{x}^{*})-\sigma_{k}h(\bm{x}^{k})+\frac{\sigma_{k}}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\Xi-\bm{A}^{\top}\bm{A}}$
$\displaystyle=\frac{S^{2}_{k-1}}{\sigma_{k}}(h(\bm{w}^{k})-h(\bm{x}^{*}))-\frac{S^{2}_{k}}{\sigma_{k}}\mathbb{E}_{k}[h(\bm{w}^{k+1})-h(\bm{x}^{*})]-\sigma_{k}(h(\bm{x}^{k})-h(\bm{x}^{*}))+\frac{\sigma_{k}}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\Xi-\bm{A}^{\top}\bm{A}},$

where the last step uses the identity
$2S_{k}=\frac{S^{2}_{k}}{\sigma_{k}}-\frac{S^{2}_{k-1}}{\sigma_{k}}+\sigma_{k}$.
Plugging this expression into the penultimate display, we arrive that

$\begin{split}\Psi(\hat{\bm{x}}^{k+1})-\Psi(\bm{x}^{*})&\leq\frac{S^{2}_{k-1}}{\sigma_{k}}(h(\bm{w}^{k})-h(\bm{x}^{*}))-\frac{S^{2}_{k}}{\sigma_{k}}\mathbb{E}_{k}[h(\bm{w}^{k+1})-h(\bm{x}^{*})]-\sigma_{k}(h(\bm{x}^{k})-h(\bm{x}^{*}))\\\
&+\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}+\Upsilon}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda-\sigma_{k}\Xi}.\end{split}$
(29)

Define vector-valued function
$\vec{\psi}(\bm{x}):=[\phi_{1}(\bm{x}_{1})+r_{1}(\bm{x}_{1}),\ldots,\psi_{d}(\bm{x}_{d})+r_{d}(\bm{x}_{d})]^{\top}\in\mathbb{R}^{d}$.
Let

$\hat{\Psi}_{k}:=\sum_{t=0}^{k}\sum_{i=1}^{d}\gamma_{i}^{k,t}\psi_{i}(\bm{x}^{t}_{i})=\sum_{t=0}^{k}\gamma^{k,t}\cdot\vec{\psi}(\bm{x}^{t}),$
(30)

where
$\gamma^{k,t}\cdot\vec{\psi}(\bm{x}^{t})=\sum_{i=1}^{d}\gamma_{i}^{k,t}\psi_{i}(\bm{x}^{t}_{i})$.
Thanks to Lemma 5, we have $\hat{\Psi}_{k}\geq\Psi(\bm{w}^{k})$ for all $k\geq
0$.

######  Lemma 6.

Let $\bm{M}:=\operatorname{\texttt{blkdiag}}[\bm{M}_{1};\ldots;\bm{M}_{d}]$
with $\bm{M}_{i}\in\mathbb{S}^{m_{i}}_{++}$. We have

$\displaystyle\mathbb{E}_{k}[\|\bm{x}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{P}\bm{M}}]$
$\displaystyle=\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{M}}+\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{(\bm{P}-\bm{I})\bm{M}},$
(31) $\displaystyle\mathbb{E}_{k}[\hat{\Psi}_{k+1}]$
$\displaystyle=(1-\theta_{k})\hat{\Psi}_{k}+\theta_{k}\Psi(\hat{\bm{x}}^{k+1}).$
(32)

###### Proof.

Eq. (31) can be easily seen by observing

$\displaystyle\mathbb{E}_{k}[\|\bm{x}^{k+1}_{i}-\bm{x}^{*}_{i}\|^{2}_{\frac{1}{\pi_{i}}\bm{M}_{i}}]$
$\displaystyle=\pi_{i}\|\hat{\bm{x}}^{k+1}_{i}-\bm{x}^{*}_{i}\|^{2}_{\frac{1}{\pi_{i}}\bm{M}_{i}}+(1-\pi_{i})\|\bm{x}^{k}_{i}-\bm{x}^{\ast}_{i}\|^{2}_{\frac{1}{\pi_{i}}\bm{M}_{i}}$
$\displaystyle=\|\hat{\bm{x}}^{k+1}_{i}-\bm{x}^{*}_{i}\|^{2}_{\bm{M}_{i}}+(\pi^{-1}_{i}-1)\|\bm{x}^{k}_{i}-\bm{x}^{*}_{i}\|^{2}_{\bm{M}_{i}}.$

Summing over all $i\in\\{1,\ldots,d\\}$ gives the result. To prove (32), we
first observe

$\mathbb{E}_{k}[\psi_{i}(\bm{x}^{k+1}_{i})]=\pi_{i}\psi_{i}(\hat{\bm{x}}^{k+1}_{i})+(1-\pi_{i})\psi_{i}(\bm{x}^{k}_{i}),$

so that
$\mathbb{E}_{k}[\vec{\psi}(\bm{x}^{k+1})]=\bm{P}^{-1}\vec{\psi}(\hat{\bm{x}}^{k+1})+(\bm{I}-\bm{P}^{-1})\vec{\psi}(\bm{x}^{k})$.
It follows

$\displaystyle\hat{\Psi}_{k+1}$
$\displaystyle=\sum_{t=0}^{k-1}\gamma^{k+1,t}\cdot\vec{\psi}(\bm{x}^{t})+\gamma^{k+1,k}\cdot\vec{\psi}(\bm{x}^{k})+\gamma^{k+1,k+1}\cdot\vec{\psi}(\bm{x}^{k+1})$
$\displaystyle=\sum_{t=0}^{k-1}(1-\theta_{k})\gamma^{k,t}\cdot\vec{\psi}(\bm{x}^{t})+\left((1-\theta_{k})\gamma^{k,k}+\theta_{k}(1-\pi^{-1})\right)\cdot\vec{r}(\bm{x}^{k})+\theta_{k}\pi^{-1}\cdot\vec{\psi}(\bm{x}^{k+1})$
$\displaystyle=(1-\theta_{k})\hat{\Psi}_{k}+\theta_{k}\left((1-\pi^{-1})\cdot\vec{\psi}(\bm{x}^{k})+\pi^{-1}\cdot\vec{\psi}(\bm{x}^{k+1})\right).$

Taking conditional expectations on both sides, it follows

$\displaystyle\mathbb{E}_{k}[\hat{\Psi}_{k+1}]$
$\displaystyle=(1-\theta_{k})\hat{\Psi}_{k}+\theta_{k}\left((1-\pi^{-1})\cdot\vec{\psi}(\bm{x}^{k})+\Psi(\hat{\bm{x}}^{k+1})+(\pi^{-1}-1)\cdot\vec{\psi}(\bm{x}^{k})\right)$
$\displaystyle=(1-\theta_{k})\hat{\Psi}_{k}+\theta_{k}\Psi(\hat{\bm{x}}^{k+1}).$

∎

The next result characterises a suitable Lyapunov function in an adapted
variable metric. We set
$\bm{B}:=\operatorname{\texttt{blkdiag}}[\bm{B}_{1};\ldots;\bm{B}_{d}]$.

######  Lemma 7.

Assume that $\theta_{0}\in(0,\min_{i\in[d]}\pi_{i}]$ and that the sequences
$(\sigma_{k})_{k\geq 0},(\tau_{k})_{k\geq 0}$ satisfy the matrix inequalities

$\displaystyle\bm{P}\bm{B}\succeq\tau_{k}(\Lambda+\sigma_{k}\Xi),$ (33)
$\displaystyle\bm{P}^{2}\bm{B}+\tau_{k}\bm{P}\Upsilon\succeq\frac{\tau_{k}\sigma_{k+1}}{\tau_{k+1}\sigma_{k}}(\bm{P}^{2}\bm{B}-\tau_{k+1}(\bm{I}-\bm{P})\Upsilon).$
(34)

Define the matrix-valued sequence $(\bm{W}_{k})_{k\geq
0}\subset\mathbb{S}^{m}_{++}$ by
$\bm{W}_{k}=\frac{\sigma_{k}}{\tau_{k}}\bm{P}^{2}\bm{B}+\sigma_{k}(\bm{P}-\bm{I})\Upsilon$,
and introduce the functions

$\displaystyle
F_{k}:=\hat{\Psi}_{k}+S_{k-1}(h(\bm{w}^{k})-h(\bm{x}^{*})),\text{ and }$ (35)
$\displaystyle
V_{k}(\bm{x}):=\frac{1}{2}\|\bm{x}^{k}-\bm{x}\|^{2}_{\bm{W}_{k}}+S_{k-1}(F_{k}-\Psi(\bm{x})),$
(36)

We have for all $\bm{x}^{*}\in\mathcal{X}$,

$\mathbb{E}_{k}\left[V_{k+1}(\bm{x}^{*})\right]\leq
V_{k}(\bm{x}^{*})-\sigma^{2}_{k}(h(\bm{x}^{k})-h(\bm{x}^{*}))-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\frac{\sigma_{k}}{\tau_{k}}(\bm{P}\bm{B}-\tau_{k}(\sigma_{k}\Xi+\Lambda))}$
(37)

###### Proof.

Using identity (31) with $\bm{M}=\bm{Q}_{k}+\Upsilon$, we get

$\displaystyle\mathbb{E}_{k}\left[\frac{1}{2}\|\bm{x}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{P}\bm{Q}^{k}+\bm{P}\Upsilon}\right]=\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}+\Upsilon}+\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{(\bm{P}-\bm{I})(\bm{Q}^{k}+\Upsilon)}$

Furthermore, eq. (32) yields

$\displaystyle\frac{S_{k}}{\sigma_{k}}\mathbb{E}_{k}[\hat{\Psi}_{k+1}-\Psi(\bm{x}^{*})]=\frac{S_{k-1}}{\sigma_{k}}\left(\hat{\Psi}_{k}-\Psi(\bm{x}^{*})\right)+\left(\Psi(\hat{\bm{x}}^{k+1})-\Psi(\bm{x}^{*})\right),$

and (29) shows

$\displaystyle\frac{S^{2}_{k}}{\sigma_{k}}\mathbb{E}_{k}[h(\bm{w}^{k+1})-h(\bm{x}^{*})]$
$\displaystyle\leq\frac{S^{2}_{k-1}}{\sigma_{k}}[h(\bm{w}^{k})-h(\bm{x}^{*})]-[\Psi(\hat{\bm{x}}^{k+1})-\Psi(\bm{x}^{*})]-\sigma_{k}[h(\bm{x}^{k})-h(\bm{x}^{*})]$
$\displaystyle+\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{Q}^{k}+\Upsilon}-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-\Lambda-\sigma_{k}\Xi}.$

Adding these three expression together, we obtain the recursion

$\begin{split}&\mathbb{E}_{k}\left[\frac{1}{2}\|\bm{x}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{P}\bm{Q}^{k}+\bm{P}\Upsilon}+\frac{S_{k}}{\sigma_{k}}\left(\hat{\Psi}_{k+1}-\Psi(\bm{x}^{*})+S_{k}(h(\bm{w}^{k+1})-h(\bm{x}^{*}))\right)\right]\\\
&\leq\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\bm{P}\bm{Q}^{k}+(\bm{P}-\bm{I})\Upsilon}+\frac{S_{k-1}}{\sigma_{k}}\left(\hat{\Psi}_{k}-\Psi(\bm{x}^{*})+S_{k-1}(h(\bm{w}^{k})-h(\bm{x}^{*}))\right)\\\
&-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{Q}^{k}-(\Lambda+\sigma_{k}\Xi)}-\sigma_{k}(h(\bm{x}^{k})-h(\bm{x}^{*})).\end{split}$
(38)

Recall that $\bm{Q}^{k}=\frac{1}{\tau_{k}}\bm{P}\bm{B}$. Condition (33)
guarantees that $\bm{Q}^{k}\succeq\Lambda+\sigma_{k}\Xi$. Multiplying both
sides of (38) by $\sigma_{k}$, we obtain

$\displaystyle\mathbb{E}_{k}\left[\frac{1}{2}\|\bm{x}^{k+1}-\bm{x}^{*}\|^{2}_{\frac{\sigma_{k}}{\tau_{k}}\bm{P}^{2}\bm{B}+\sigma_{k}\bm{P}\Upsilon}+S_{k}\left(\hat{\Psi}_{k+1}-\Psi(\bm{x}^{*})+S_{k}(h(\bm{w}^{k})-h(\bm{x}^{*}))\right)\right]$
$\displaystyle\leq\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\frac{\sigma_{k}}{\tau_{k}}\bm{P}^{2}\bm{B}+\sigma_{k}(\bm{P}-\bm{I})\Upsilon}+S_{k-1}\left(\hat{\Psi}_{k}-\Psi(\bm{x}^{*})+S_{k-1}(h(\bm{w}^{k})-h(\bm{x}^{*}))\right)$
$\displaystyle-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\frac{\sigma_{k}}{\tau_{k}}\bm{P}\bm{B}-\sigma_{k}(\Lambda+\sigma_{k}\Xi)}-\sigma^{2}_{k}(h(\bm{x}^{k})-h(\bm{x}^{*})).$

Under condition (34), it follows that

$\|\bm{x}^{k+1}-\bm{x}^{*}\|^{2}_{\bm{W}_{k+1}}\leq\|\bm{x}^{k+1}-\bm{x}^{*}\|^{2}_{\frac{\sigma_{k}}{\tau_{k}}\bm{P}^{2}\bm{B}+\sigma_{k}\bm{P}\Upsilon}.$

Exploiting this relation in the penultimate display, using definitions (35)
and (36), one readily obtains

$\mathbb{E}_{k}\left[V_{k+1}(\bm{x}^{*})\right]\leq
V_{k}(\bm{x}^{*})-\sigma^{2}_{k}(h(\bm{x}^{k})-h(\bm{x}^{*}))-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\frac{\sigma_{k}}{\tau_{k}}(\bm{P}\bm{B}-\tau_{k}(\sigma_{k}\Xi+\Lambda))}$

∎

###  4.2 The Convex case

Suppose that $\Upsilon=0$, so that
$\bm{W}_{k}=\frac{\sigma_{k}}{\tau_{k}}\bm{P}^{2}\bm{B}$. In this case, the
matrix inequality (34) is satisfied for any sequences $(\sigma_{k})_{k\geq
0},(\tau_{k})_{k\geq 0}$ satisfying
$\frac{\sigma_{k}}{\tau_{k}}\geq\frac{\sigma_{k+1}}{\tau_{k+1}}$. In
particular, this can be realised with the policy $\sigma_{k}=\sigma\tau_{k}$
for all $k\geq 0$, and some constant $\sigma\in(0,\min_{i\in[d]}\pi_{i}]$.
This implies $\theta_{k}=\frac{\sigma}{1+k\sigma}$ for all $k\geq 0$. This
specification satisfies all conditions needed for Lemma 5 to hold. We further
set $\tau_{k}\equiv 1$. It only remains to see if matrix inequality (33) is
satisfied. In fact, this condition gives us a restriction on $\sigma$, reading
as

$\bm{P}\bm{B}\succeq\sigma\Xi+\Lambda.$ (39)

We show later in this section that this matrix inequality is achievable. In
this regime, we will show that we get an $O(1/k)$ iteration complexity
estimate of the averaged sequence $(\bm{w}^{k})_{k}$, in terms of the
objective function value. This extends the results reported in [27, 28, 26] to
general random block-coordinate activation schemes. The following Assumptions
shall be in place throughout this section:

######  Assumption 1.

Problem (P) admits a saddle point, i.e. a pair
$(\bm{x}^{\ast},\bm{y}^{\ast})\in\mathbb{R}^{m}\times\mathbb{R}^{m}$
satisfying (8).

######  Assumption 2.

The solution set of problem (P) $\mathcal{X}^{\ast}$ is nonempty.

Thanks to eq. (26), we know that
$\bm{w}^{k}\in\texttt{Conv}(\bm{x}^{0},\ldots,\bm{x}^{k})\subseteq\texttt{dom}(r)$.
Therefore,

$\displaystyle r(\bm{w}^{k})-r(\bm{x}^{\ast})$
$\displaystyle\geq\langle{-\nabla\phi(\bm{x}^{*})-\bm{A}^{\top}\bm{A}\bm{y}^{*}},{\bm{w}^{k}-\bm{x}^{*}}\rangle$
$\displaystyle\geq\phi(\bm{x}^{*})-\phi(\bm{w}^{k})+\langle{\bm{A}\bm{y}^{*}},{\bm{A}(\bm{x}^{*}-\bm{w}^{k})}\rangle$

Whence,

$\Psi(\bm{w}^{k})-\Psi(\bm{x}^{*})\geq-\|\bm{A}\bm{y}^{*}\|\cdot\|\bm{A}(\bm{w}^{k}-\bm{x}^{*})\|=-\delta\sqrt{h(\bm{w}^{k})-h(\bm{x}^{*})},$
(40)

where $\delta:=\sqrt{2}\|A\bm{y}^{*}\|$. We assume that $\tau_{k}\equiv 1$ and
$\sigma_{k}=\sigma>0$. The energy functions (35) and (36) take the form
$F_{k}=\hat{\Psi}_{k}+(1+(k-1)\sigma)(h(\bm{w}^{k})-h(\bm{x}^{*}))$, and
$V_{k}(\bm{x}):=\frac{1}{2}\|\bm{x}^{k}-\bm{x}\|^{2}_{\bm{P}^{2}\bm{B}}+(1+(k-1)\sigma)(F_{k}-\Psi(\bm{x}))$.

######  Lemma 8.

Suppose that Assumption 1 applies. Then, the process
$\left(S_{k-1}(F_{k}-\Psi(\bm{x}^{*}))\right)_{k\geq 0}$ is almost surely
bounded.

###### Proof.

Note that

$\displaystyle S_{k-1}(F_{k}-\Psi(\bm{x}^{*}))$
$\displaystyle=S_{k-1}(\hat{\Psi}_{k}-\Psi(\bm{x}^{*}))+S^{2}_{k-1}(h(\bm{w}^{k})-h(\bm{x}^{*}))$
$\displaystyle\geq
S_{k-1}(\Psi(\bm{w}^{k})-\Psi(\bm{x}^{*}))+S^{2}_{k-1}(h(\bm{w}^{k})-h(\bm{x}^{*}))$
$\displaystyle\geq-\delta
S_{k-1}\sqrt{h(\bm{w}^{k})-h(\bm{x}^{*})}+S^{2}_{k-1}(h(\bm{w}^{k})-h(\bm{x}^{*}))$
$\displaystyle=\frac{S^{2}_{k-1}}{2}(h(\bm{w}^{k})-h(\bm{x}^{*}))+\left(\frac{S^{2}_{k-1}}{2}(h(\bm{w}^{k})-h(\bm{x}^{*}))-\delta
S_{k-1}\sqrt{h(\bm{w}^{k})-h(\bm{x}^{*})}\right)$

The convex function $t\mapsto\frac{S^{2}}{2}t^{2}-\delta St$ attains the
global minimum at the value $-\frac{\delta^{2}}{2}$. Hence,

$S_{k-1}(F_{k}-\Psi(\bm{x}^{*}))\geq\frac{S^{2}_{k-1}}{2}(h(\bm{w}^{k})-h(\bm{x}^{*}))-\frac{\delta^{2}}{2}\geq-\frac{\delta^{2}}{2}.$
(41)

The last inequality uses the fact that $\bm{x}^{\ast}\in\mathcal{X}$, so that
$h(\bm{w}^{k})\geq h(\bm{x}^{\ast})$. This completes the proof. ∎

We are now in the position to give the proof of the first main result of this
paper.

######  Theorem 1.

Suppose that Assumptions 1 and 2 apply. Consider the sequence
$(\bm{z}^{k},\bm{w}^{k},\bm{x}^{k})_{k\geq 0}$ generated by Algorithm 1 with
$\sigma_{k}\equiv\sigma$ and $\tau_{k}\equiv 1$ for all $k\geq 1$. Then, the
following holds:

  * (i)

We have $h(\bm{x}^{k})-h(\bm{x}^{*})=o(1/k)$ and
$h(\bm{w}^{k})-h(\bm{x}^{*})=O(1/k^{2})$ a.s.

  * (ii)

$(\bm{x}^{k})_{k}$ and $(\bm{w}^{k})_{k}$ converge a.s. to a random variable
with values in $\mathcal{X}^{\ast}$.

  * (iii)

$\lvert\Psi(\bm{w}^{k})-\Psi(\bm{x}^{*})\rvert=O(1/k)$ a.s.

###### Proof.

(i) From Lemma 8, it follows that the process
$\mathcal{E}_{k}:=V_{k}(\bm{x}^{*})+\frac{\delta^{2}}{2}$ is non-negative and
satisfies the recursion

$\mathbb{E}_{k}[\mathcal{E}_{k+1}]\leq\mathcal{E}_{k}-\sigma^{2}(h(\bm{x}^{k})-h(\bm{x}^{*}))-\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\sigma(\bm{P}\bm{B}-(\sigma\Xi+\Lambda))}.$
(42)

Lemma 10 implies that $(\mathcal{E}_{k})_{k\geq 0}$ converges a.s. to a
limiting random variable $\mathcal{E}_{\infty}\in[0,\infty)$ and

$\sum_{k=0}^{\infty}\left[\sigma^{2}(h(\bm{x}^{k})-h(\bm{x}^{*}))+\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\sigma(\bm{P}\bm{B}-(\sigma\Xi+\Lambda))}\right]<\infty\quad\text{a.s.}$

Since $h(\bm{x}^{k})-h(\bm{x}^{*})\geq 0$, it follows
$\lim_{k\to\infty}\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\sigma(\bm{P}\bm{B}-(\sigma\Xi+\Lambda))}=0$.
Consequently, the sequence
$\left(\frac{1}{2}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\sigma(\bm{P}\bm{B}-(\sigma\Xi+\Lambda))}\right)_{k\geq
0}$ is bounded, and therefore $h(\bm{x}^{k})-h(\bm{x}^{*})=o(1/k)$. By
definition of the energy function $V_{k}(\bm{x}^{*})$, the sequence
$(\bm{x}^{k})_{k\geq 0}$ is bounded and (41) yields
$h(\bm{w}^{k})-h(\bm{x}^{*})\geq-\frac{\delta^{2}}{S^{2}_{k-1}}.$ Hence,
$h(\bm{w}^{k})-h(\bm{x}^{*})=O(1/k^{2})$ a.s. Moreover, it is easy to see that
$\lim_{k\to\infty}\sum_{t=0}^{k}\gamma_{i}^{k,t}=1$ and
$\lim_{k\to\infty}\gamma_{i}^{k,t}=0$ for all $i\in[d]$. Hence, by Lemma 5 and
the Silverman-Toeplitz theorem, the sequence $(\bm{w}^{k})_{k}$ converges to
the same limit point as $(\bm{x}^{k})_{k}$, a.s.

By definition of the energy function $V_{k}$, the sequence
$(\bm{x}^{k})_{k\geq 0}$ is a.s. bounded. We show next that all accumulation
bounds are contained in the set of solutions $\mathcal{X}^{\ast}$ of (P).

Let $\Omega_{1}$ be a set of probability 1 on which the subsequence
$(\bm{x}^{k_{j}}(\omega))_{j\in\mathbb{N}}$ converges to a limiting random
variable $\bar{\bm{x}}(\omega)$ for all $\omega\in\Omega_{1}$. It follows from
$\lim_{j\to\infty}\|\hat{\bm{x}}^{k_{j}}(\omega)-\bm{x}^{k_{j}}(\omega)\|=0$
on $\Omega_{1}$ that $(\hat{\bm{x}}^{k_{j}})_{j\in\mathbb{N}}$ converges to
$\bar{\bm{x}}$ as well on $\Omega_{1}$. Furthermore,

$\displaystyle h(\bm{z}^{k})-h(\bm{x}^{\ast})$
$\displaystyle=\theta_{k}[h(\bm{w}^{k})-h(\bm{x}^{\ast})]+(1-\theta_{k})[h(\bm{x}^{k})-h(\bm{x}^{\ast})]$
$\displaystyle-\frac{\theta_{k}(1-\theta_{k})}{2}\|\bm{A}(\bm{x}^{k}-\bm{w}^{k})\|^{2}$
$\displaystyle=O(k^{-2}),$

and for all $\bm{x}^{\ast}\in\mathcal{X}$,

$\displaystyle\langle{\nabla
h(\bm{z}^{k})},{\hat{\bm{x}}^{k+1}-\bm{x}^{\ast}}\rangle$
$\displaystyle=\langle{\bm{A}^{\top}\bm{A}(\bm{z}^{k}-\bm{x}^{\ast})},{\hat{\bm{x}}^{k+1}-\bm{x}^{\ast}}\rangle$
$\displaystyle\leq\|\bm{A}(\bm{z}^{k}-\bm{x}^{\ast})\|\cdot\|\bm{A}(\hat{\bm{x}}^{k+1}-\bm{x}^{\ast})\|$
$\displaystyle\leq
2\sqrt{h(\bm{z}^{k})-h(\bm{x}^{\ast})}\cdot\sqrt{h(\hat{\bm{x}}^{k+1})-h(\bm{x}^{\ast})}=o(k^{-1})$

Hence, $\lim_{k\to\infty}S_{k}\langle{\nabla
h(\bm{z}^{k})},{\hat{\bm{x}}^{k+1}-\bm{x}^{\ast}}\rangle=0$ for all
$\bm{x}^{\ast}\in\mathcal{X}$ a.s. Using Lemma 4 at the point
$\bm{x}^{\ast}\in\mathcal{X}^{\ast}\subset\mathcal{X}$ shows that

$\displaystyle
R(\hat{\bm{x}}^{k_{j+1}})+\langle{\nabla\Phi(\bm{x}^{k_{j}})},{\hat{\bm{x}}^{k_{j+1}}-\bm{x}^{k_{j}}}\rangle+S_{k}\langle{\nabla
h(\bm{z}^{k_{j}})},{\hat{\bm{x}}^{k_{j+1}}-\bm{z}^{k_{j}}}\rangle+\frac{1}{2}\|\hat{\bm{x}}^{k_{j+1}}-\bm{x}^{k_{j}}\|^{2}_{\bm{Q}}$
$\displaystyle\leq
R(\bm{x}^{\ast})+\langle{\nabla\Phi(\bm{x}^{k_{j}})},{\bm{x}^{\ast}-\bm{x}^{k_{j}}}\rangle+S_{k}\langle{\nabla
h(\bm{z}^{k_{j}})},{\bm{x}^{\ast}-\bm{z}^{k_{j}}}\rangle+\frac{1}{2}\|\hat{\bm{x}}^{\ast}-\bm{x}^{k_{j}}\|^{2}_{\bm{Q}}-\frac{1}{2}\|\bm{x}^{\ast}-\hat{\bm{x}}^{k_{j+1}}\|_{\bm{Q}}.$

Let $j\to\infty$, and use the lower semi-continuity of the function
$R(\cdot)$, one arrives that the inequality

$R(\bar{\bm{x}})\leq
R(\bm{x}^{\ast})+\langle{\nabla\Phi(\bar{\bm{x}})},{\bm{x}^{\ast}-\bar{\bm{x}}}\rangle.$

Since $h(\bm{x}^{k_{j}})-h(\bm{x}^{\ast})\to 0$, the limit point
$\bar{\bm{x}}$ is an element of $\mathcal{X}$. Convexity of $\Phi$ in turn
yields
$\Phi(\bm{x}^{\ast})\geq\Phi(\bar{\bm{x}})+\langle{\nabla\Phi(\bar{\bm{x}})},{\bm{x}^{\ast}-\bar{\bm{x}}}\rangle$.
We conclude that $\Psi(\bar{\bm{x}})\leq\Psi(\bm{x}^{\ast})$ and
$\bar{\bm{x}}\in\mathcal{X}$. Therefore, $\bar{\bm{x}}\in\mathcal{X}^{\ast}$.
(ii) We next show that cluster points are unique and thereby demonstrate
almost sure convergence of the last iterate. We argue by contradiction.
Suppose there are converging subsequences $(\bm{x}^{k})_{k\in\mathcal{K}_{1}}$
and $(\bm{x}^{k})_{k\in\mathcal{K}_{2}}$ with limit points $\bm{x}^{\prime}$
and $\bm{x}^{\prime\prime}$, respectively. Hence,
$\Psi(\bm{x}^{\prime})=\Psi(\bm{x}^{\prime\prime})\equiv\Psi^{\ast}$.
Following the above argument, we see that
$\bm{x}^{\prime},\bm{x}^{\prime\prime}\in\mathcal{X}^{\ast}$ a.s. The point
$\bm{x}^{\ast}$ chosen in the previous argument was arbitrary, so we can
replace it by $\bm{x}^{\prime}$. Let us simplify notation by setting
$a_{k}:=S_{k}(F_{k}-\Psi^{\ast})$. Since $V_{k}(\bm{x}^{\prime})$ converges
almost surely, we see

$\displaystyle\lim_{k\to\infty}V_{k}(\bm{x}^{\prime})$
$\displaystyle=\lim_{k\to\infty,k\in\mathcal{K}_{1}}V_{k}(\bm{x}^{\prime})=\lim_{k\to\infty,k\in\mathcal{K}_{1}}\left(\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{\prime}\|^{2}_{\bm{W}}+a_{k}\right)$
$\displaystyle=\lim_{k\to\infty,k\in\mathcal{K}_{1}}a_{k}.$

Similarly,

$\displaystyle\lim_{k\to\infty}V_{k}(\bm{x}^{\prime})$
$\displaystyle=\lim_{k\to\infty,k\in\mathcal{K}_{2}}V_{k}(\bm{x}^{\prime})=\lim_{k\to\infty,k\in\mathcal{K}_{2}}\left(\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{\prime}\|^{2}_{\bm{W}}+a_{k}\right)$
$\displaystyle=\|\bm{x}^{\prime\prime}-\bm{x}^{\prime}\|_{\bm{W}}+\lim_{k\to\infty,k\in\mathcal{K}_{2}}a_{k}.$

We conclude that

$\lim_{k\to\infty,k\in\mathcal{K}_{1}}a_{k}=\|\bm{x}^{\prime\prime}-\bm{x}^{\prime}\|_{\bm{W}}+\lim_{k\to\infty,k\in\mathcal{K}_{2}}a_{k}.$

Repeating the same analysis for $\bm{x}^{\prime\prime}$ instead of
$\bm{x}^{\prime}$, we get

$\lim_{k\to\infty,k\in\mathcal{K}_{2}}a_{k}=\|\bm{x}^{\prime\prime}-\bm{x}^{\prime}\|_{\bm{W}}+\lim_{k\to\infty,k\in\mathcal{K}_{1}}a_{k}.$

Combining these two equalities, we see that
$\bm{x}^{\prime}=\bm{x}^{\prime\prime}$. Therefore, the whole sequence
$(\bm{x}^{k})_{k}$ converges pointwise almost surely to a solution. Lemma 5
yields the same assertion for $(\bm{w}^{k})_{k}$.
(iii) Taking expectations on both sides of (37) and iterating the thus
obtained expression gives

$\mathbb{E}\left[\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{*}\|^{2}_{\bm{W}}+S_{k-1}(\hat{\Psi}_{k}-\Psi(\bm{x}^{*}))+S^{2}_{k-1}(h(\bm{w}^{k})-h(\bm{x}^{*}))\right]\leq\frac{1}{2}\|\bm{x}^{0}-\bm{x}^{*}\|^{2}_{\bm{W}}.$

From here, we deduce

$S_{k-1}\mathbb{E}[\hat{\Psi}_{k}-\Psi(\bm{x}^{*})]\leq\frac{1}{2}\|\bm{x}^{0}-\bm{x}^{*}\|^{2}_{\bm{W}}$

which gives $\hat{\Psi}_{k}-\Psi(\bm{x}^{*})\leq O(1/k)$ a.s. Since
$\hat{\Psi}_{k}\geq\Psi(\bm{w}^{k})$, it follows
$\Psi(\bm{w}^{k})-\Psi(\bm{x}^{*})\leq O(1/k)$ a.s. To obtain a lower bound,
it suffices to apply eq. (40) and use the results from (i) to obtain

$\Psi(\bm{w}^{k})-\Psi(\bm{x}^{\ast})\geq-\delta\sqrt{h(\bm{w}^{k})-h(\bm{x}^{\ast})}=O(1/k)\quad\text{a.s.}$

∎

We conclude this section with a concrete example showing that the matrix
inequality (39) can be satisfied.

######  Example 3.

Assume that
$\bm{\Xi}=\operatorname{\texttt{blkdiag}}(\pi_{1}^{-1}\bm{A}_{1}^{\top}\bm{A}_{1},\ldots,\pi_{d}^{-1}\bm{A}_{d}^{\top}\bm{A}_{d})$.
In this case eq. (39) can be decomposed to the block specific condition
$\pi^{-1}_{i}\bm{B}_{i}\succeq\Lambda_{i}+\frac{\sigma}{\pi_{i}}\bm{A}_{i}^{\top}\bm{A}_{i}$.
Let us assume further that $\Lambda_{i}=L_{i}\bm{I}_{m_{i}}$, for a scalar
$L_{i}>0$. If we choose $\bm{B}_{i}=\pi_{i}(\lambda_{i}+L_{i})\bm{I}_{m_{i}}$,
where $\lambda_{i}=\|\bm{A}_{i}\|_{2}$. Then, a sufficient condition for
satisfying the matrix inequality is $\lambda_{i}+L_{i}\geq
L_{i}+\frac{\sigma}{\pi_{i}}\lambda_{i}$ for all $i\in[d]$. This inequality
can be satisfied by choosing $\sigma=\min_{i\in[d]}\pi_{i}$, which is the
largest value for a given set of activation probabilities, that is compatible
with Lemma 5.

###  4.3 The strongly convex case

In this section we study the performance of Algorithm 1 in the strongly convex
regime.

######  Assumption 3.

$\Upsilon\succ 0$.

The main challenge we need to overcome is to design step size sequences which
obey the matrix inequalities (33) and (34). Let us focus on (33) first and
show how it can be restated in a more symmetric way. To that end, define the
matrix $\bm{B}=\bm{P}^{-2}\Upsilon$ and
$\bm{M}^{1/2}:=\bm{P}^{-1/2}\Upsilon^{1/2}$, so that
$\bm{M}=\bm{P}^{-1}\Upsilon$. In terms of this new parameter matrix, we see

$\bm{P}^{-1}\Upsilon\succeq\tau_{k}\Lambda+\sigma_{k}\tau_{k}\bm{\Xi}\iff\frac{1}{\tau_{k}}\bm{I}\succeq\bm{M}^{-1/2}\Lambda\bm{M}^{-1/2}+\sigma_{k}\bm{M}^{-1/2}\bm{\Xi}\bm{M}^{-1/2}.$
(43)

It is easy to see that $\bm{\Xi}$ is symmetric. Furthermore, for
$\bm{u}\in\mathbb{R}^{m}\setminus\\{0\\}$, we see
$\bm{u}^{\top}\bm{\Xi}\bm{u}=\mathbb{E}[\|\bm{A}\bm{t}\|^{2}]\geq 0$, where
$\bm{t}:=\bm{P}\bm{E}_{k}\bm{u}$ is a random vector in $\mathbb{R}^{m}$ with
mean $\bm{u}$. It follows $\bm{\Xi}\in\mathbb{S}^{m}_{+}$. This suggests to
relate the step size parameters to the spectra of the involved matrices. The
rest of this section assumes the following parametric model on the coupling
between $\sigma_{k}$ and $\tau_{k}$.

######  Assumption 4.

The sequence $(\tau_{k})_{k}$ is non-increasing and positive. They are related
by the coupling equation

$\sigma_{k}\tau_{k}=\alpha-\beta\tau_{k}\qquad\forall k\geq 0,$ (44)

where $\alpha$ and $\beta$ are positive numbers.

We note that this coupling relation automatically means that
$\sigma_{k}=\alpha/\tau_{k}-\beta$ is non-decreasing. A sufficient condition
to satisfy eq. (43) is

$\frac{1}{\tau_{k}}\geq\lambda_{\max}(\bm{M}^{-1/2}\Lambda\bm{M}^{-1/2})+\sigma_{k}\lambda_{\max}(\bm{M}^{-1/2}\Xi\bm{M}^{-1/2}).$

Using the model (44), this boils down to

$\frac{1}{\tau_{k}}(1-\alpha\lambda_{\max}(\bm{M}^{-1/2}\bm{\Xi}\bm{M}^{-1/2}))\geq\lambda_{\max}(\bm{M}^{-1/2}\Lambda\bm{M}^{-1/2})-\beta\lambda_{\max}(\bm{M}^{-1/2}\bm{\Xi}\bm{M}^{-1/2}).$

Choosing

$\alpha:=\lambda_{\max}(\bm{M}^{-1/2}\bm{\Xi}\bm{M}^{-1/2})^{-1},$ (45)

we get

$\beta\geq\frac{\lambda_{\max}(\bm{M}^{-1/2}\Lambda\bm{M}^{-1/2})}{\lambda_{\max}(\bm{M}^{-1/2}\bm{\Xi}\bm{M}^{-1/2})}=\alpha\lambda_{\max}(\bm{M}^{-1/2}\Lambda\bm{M}^{-1/2}).$

We make the choice

$\beta=\alpha\lambda_{\max}(\bm{M}^{-1/2}(\Lambda+\Upsilon)\bm{M}^{-1/2}),$
(46)

so that the coupling equation (44) gives rise to a step size process
satisfying matrix inequality (33).

######  Remark 4.

To get some intuition behind the derived conditions for $\alpha$ and $\beta$,
consider the independent sampling case under which the operator $\bm{\Xi}$
factorises to
$\bm{\Xi}=\operatorname{\texttt{blkdiag}}(\frac{1}{\pi_{1}}\bm{A}_{1}^{\top}\bm{A}_{1},\ldots,\frac{1}{\pi_{d}}\bm{A}_{d}^{\top}\bm{A}_{d})$.
Further, suppose
$\Lambda=\operatorname{\texttt{blkdiag}}[L_{1}\bm{I}_{m_{1}};\ldots;L_{d}\bm{I}_{m_{d}}]$.
Then, the above choice says that

$\displaystyle\alpha=\left(\max_{i\in[d]}\frac{1}{\mu_{i}\pi_{i}^{2}}\|\bm{A}_{i}\|_{2}^{2}\right)^{-1},\text{and
}$ (47)
$\displaystyle\kappa:=\frac{\beta}{\alpha}=\max_{i\in[d]}\frac{L_{i}+\mu_{i}}{\pi_{i}\mu_{i}}.$
(48)

Clearly, $\kappa$ is related to the condition number of problem (P), and
satisfies $\kappa\geq\frac{1}{\pi_{i}}$ for all $i\in[d]$. We point out that
the explicit expression for $\kappa$ does not rely on the independent sampling
assumption. Only the quantification of the parameter $\alpha$ exploits this
special structure.

Let us fix these identified values for $\alpha$ and $\beta$, and continue our
derivation with matrix inequality (34). Let us choose
$\bm{B}=\bm{P}^{-2}\Upsilon=\operatorname{\texttt{blkdiag}}[\pi_{1}^{2}\mu_{1}\bm{I}_{m_{1}};\ldots;\pi_{d}^{2}\mu_{d}\bm{I}_{m_{d}}]$.
Using the block structure, and the specification of the sequence
$(\sigma_{k})_{k}$, it can be verified that (34) reduces to the scalar-valued
inequality

$\tau^{2}_{k+1}\left[(\alpha-\beta\tau_{k})(\pi_{i}+\tau_{k})+\beta(1-\pi_{i})\tau^{2}_{k}\right]-\tau_{k+1}\tau^{2}_{k}[\alpha(1-\pi_{i})-\beta\pi_{i}]-\tau^{2}_{k}\alpha\pi_{i}\geq
0.$ (49)

This defines a quadratic inequality in $\tau_{k+1}$ of the form
$c_{i,1}(\tau_{k})\tau^{2}_{k+1}+c_{i,2}(\tau_{k})\tau_{k+1}-c_{i,3}(\tau_{k})\geq
0$, with coefficients

$\displaystyle c_{i,1}(t):=(\alpha-\beta t)(\pi_{i}+t)+\beta(1-\pi_{i})t^{2},$
$\displaystyle c_{i,2}(t):=t^{2}[\beta\pi_{i}-\alpha(1-\pi_{i})],$
$\displaystyle c_{i,3}(t):=t^{2}\alpha\pi_{i}$

This suggests to use the recursive step-size policy

$\tau_{k+1}=\max_{i\in[d]}\frac{\frac{\tau^{2}_{k}}{2}(1/\pi_{i}-\kappa-1)+\tau_{k}\sqrt{\left(1+(1/\pi_{i}-\kappa)\frac{\tau_{k}}{2}\right)^{2}-\left(\frac{2-\pi_{i}}{\pi_{i}}+2\kappa\right)\frac{\tau^{4}_{k}}{4}}}{1+\tau_{k}(1/\pi_{i}-\kappa)-\kappa\tau^{2}_{k}}.$
(50)

We study the qualitative properties of the so produced sequence
$(\tau_{k})_{k}$ in detail in Appendix B, under the following assumptions:

######  Assumption 5.

For all $i\in[d]$ we have $\Lambda_{i}=L_{i}\bm{I}_{m_{i}}$.

######  Assumption 6.

The sampling is uniform: $\pi_{i}\equiv\pi\in(0,1)$ for all $i\in[d]$.

Uniform sampling is a very common sampling scheme employed in the literature.
It contains as special cases the single coordinate activation scheme (i.e.
$\pi_{i}=\frac{1}{d}$), and the uniform sampling scheme. In particular, the
frequently employed $m$-nice sampling (cf. Section 2.6) is contained in this
framework. Under this assumption we prove that if $\tau_{0}$ is chosen
sufficiently small, then the sequence $(\tau_{k})_{k}$ exhibits the same
qualitative behaviour as a Nesterov accelerated method [35, 10]. We therefore
can prove accelerated rates for our scheme in the strongly convex case.

######  Theorem 2.

Suppose that Assumptions 1-6 apply. Consider the sequence
$(\bm{z}^{k},\bm{w}^{k},\bm{x}^{k})_{k\geq 0}$ generated by Algorithm 1 with
the following step size policy constructed via eqs. (44), (45), (46). Then,
the following holds:

  * (i)

We have $h(\bm{x}^{k})-h(\bm{x}^{*})=o(k^{-2})$ and
$h(\bm{w}^{k})-h(\bm{x}^{*})=O(k^{-4})$ a.s.

  * (ii)

$(\bm{x}^{k})_{k}$ and $(\bm{w}^{k})_{k}$ converge a.s. to the solution
$\mathcal{X}^{\ast}=\\{\bar{\bm{x}}\\}$.

  * (iii)

$\lvert\Psi(\bm{w}^{k})-\Psi(\bm{x}^{*})\rvert=O(k^{-2})$ a.s.

###### Proof.

(i) If a Lagrange multiplier exists then Lemma 8 still applies and guarantees
that the function $V_{k}(\bm{x}^{\ast})+\frac{\delta^{2}}{2}\geq 0$ for all
$k\geq 0$, and the recursion (42) applies, which reads in the present case

$\mathbb{E}_{k}[\mathcal{E}_{k+1}]\leq\mathcal{E}_{k}-\sigma^{2}_{k}(h(\bm{x}^{k})-h(\bm{x}^{*}))-\frac{\sigma_{k}}{2\tau_{k}}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{P}^{-1}\Upsilon-\tau_{k}(\sigma_{k}\Xi+\Lambda)}.$

Since $\bm{P}^{-1}\Upsilon-\tau_{k}(\sigma_{k}\Xi+\Lambda)\succeq 0$, we can
upper bound the right-hand side of the above display as

$\mathbb{E}_{k}[\mathcal{E}_{k+1}]\leq\mathcal{E}_{k}-\sigma^{2}_{k}(h(\bm{x}^{k})-h(\bm{x}^{*}))-\frac{\sigma_{k}}{2\tau_{k}}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{P}^{-1}\Upsilon}.$
(51)

The supermartingale convergence theorem implies that

$\lim_{k\to\infty}\sigma^{2}_{k}(h(\bm{x}^{k})-h(\bm{x}^{*}))=0\text{ and
}\lim_{k\to\infty}\frac{\sigma_{k}}{2\tau_{k}}\|\hat{\bm{x}}^{k+1}-\bm{x}^{k}\|^{2}_{\bm{P}^{-1}\Upsilon}=0.$

Since $\sigma_{k}=O(k)$, it follows
$h(\bm{x}^{k})-h(\bm{x}^{\ast})=o(k^{-2})$. Furthermore, the sequence
$(\bm{x}^{k})_{k}$ is bounded and eq. (41) yields
$h(\bm{w}^{k})-h(\bm{x}^{\ast})=O(k^{-4})$, since $S_{k}=O(k^{2})$ as shown in
Appendix B. Moreover, $(\bm{w}^{k})_{k}$ and $(\bm{x}^{k})_{k}$ share the same
limit.
(ii) We can repeat the arguments of Theorem 1, to conclude that accumulation
points of $(\hat{\bm{x}}^{k})_{k}$ and $(\bm{x}^{k})_{k}$ converge to the same
random variable with values in $\mathcal{X}^{\ast}$. However, since the
problem is strongly convex, we must have
$\mathcal{X}^{\ast}=\\{\bar{\bm{x}}\\}$ for some $\bar{\bm{x}}\in\mathcal{X}$.
Hence, the sequences actually converge and the common limit point is
$\bar{\bm{x}}$ a.s. (iii) We next show convergence rates in terms of the
objective function gap. Iterating eq. (51), it follows
$\mathbb{E}[V_{k}(\bm{x}^{\ast})]\leq\frac{1}{2}\|\bm{x}^{0}-\bm{x}^{\ast}\|^{2}_{\bm{W}_{0}}$.
This is equivalent to

$\mathbb{E}\left[\frac{1}{2}\|\bm{x}^{k}-\bm{x}^{\ast}\|^{2}_{\bm{W}_{k}}+S_{k-1}(\hat{\Psi}_{k}-\Psi(\bm{x}^{\ast}))+S^{2}_{k-1}(h(\bm{w}^{k})-h(\bm{x}^{\ast}))\right]\leq\frac{1}{2}\|\bm{x}^{0}-\bm{x}^{\ast}\|^{2}_{\bm{W}_{0}}.$

In particular,

$\mathbb{E}[\hat{\Psi}_{k}-\Psi(\bm{x}^{\ast})]\leq\frac{1}{2S_{k-1}}\|\bm{x}^{0}-\bm{x}^{\ast}\|^{2}_{\bm{W}_{0}}=O(k^{-2}).$

Furthermore, $\hat{\Psi}_{k}\geq\Psi(\bm{w}^{k})$, so that via eq. (40) we
obtain

$\displaystyle\mathbb{E}[\Psi(\bm{w}^{k})-\Psi(\bm{x}^{\ast})]\leq
O(k^{-2})\text{ and
}\mathbb{E}[\Psi(\bm{w}^{k})-\Psi(\bm{x}^{\ast})]\geq-\delta\mathbb{E}\left[\sqrt{h(\bm{w}^{k})-h(\bm{x}^{\ast})}\right].$

Combining these bounds, it follows
$\mathbb{E}[\Psi(\bm{w}^{k})-\Psi(\bm{x}^{\ast})]=O(k^{-2})$. ∎

##  5 Application to power systems

In this section we apply Algorithm 1 to the distributed coordination of an
energy grid in an AC-optimal power flow formulation. Specifically, we
illustrate how our block-coordinate descent method provides an efficient way
to replicate the transmission-level locational marginal prices to the
distribution level.111A preliminary version of this application is studied in
our conference paper [6]. We refer to this paper for background and further
motivation.

###  5.1 Power flow model

Consider a radial network with $N$ buses $\mathcal{N}=\\{1,\dots,N\\}$,
excluding the slack node $0$. Hence, the distribution network assumes the
structure of a tree graph and we identify the transmission lines with the
label of the parental node. The network is optimised over a time window
$\mathcal{T}=\\{1,\dots,T\\}$. We use $\bm{p}_{n}=(p_{n,1},\dots,p_{n,T})$ and
$\bm{q}_{n}=(q_{n,1},\dots,q_{n,T})$ to denote active and reactive power net
consumption at bus $n$ at each time point $t\in\mathcal{T}$. We let
$\operatorname{\mathsf{i}}=\sqrt{-1}$ and write a complex number
$z\in\mathbb{C}$ as $\text{Re}(z)+\operatorname{\mathsf{i}}\text{Im}(z)$.
Thus, $p_{n,t}<0$ means that there is production of energy at bus $n$ at time
$t$. At the slack node $n=0$, we assume that power will only be generated and
there is no consumption, i.e. $p_{0,t}\leq 0$. At node $n$, we denote by
$f_{n}=(f_{n,t})_{t\in\mathcal{T}}$ and $g_{n}=(g_{n,t})_{t\in\mathcal{T}}$
the real and reactive power flows, so that
$f_{n}+\operatorname{\mathsf{i}}g_{n}$ is the complex power injected into node
$n$. By Ohm’s law, we have $g_{n}=\text{Im}(V_{n}\bar{I}_{n})$ and
$f_{n}=\text{Re}(V_{n}\bar{I}_{n})$. We let
$\ell_{n}=(\ell_{n,t})_{t\in\mathcal{T}}=\lvert V_{n}\rvert$ is the squared
current magnitude, $R_{n}$ and $X_{n}$ denotes the series resistance and
reactance. Hence, $z_{n}=R_{n}+\operatorname{\mathsf{i}}X_{n}$ is the series
impedance, and $y_{n}=z_{n}^{-1}$ the shunt admittance.

The load flow equations using the BFM [1, 30], without explicit consideration
of transmission tap ratios, are as follows:

$\displaystyle\bm{f}_{n}\,{-}\,\\!\\!\sum\limits_{m:m_{-}=n}\\!\\!(\bm{f}_{m}\,{-}\,R_{m}\bm{l}_{m})\,{+}\,\bm{p}_{n}\,{+}\,G_{n}\bm{v}_{n}\,{=}\,0$
$\displaystyle[\bm{y}^{\textup{p}}_{n}]$ (52a)
$\displaystyle\bm{g}_{n}\,{-}\,\\!\\!\sum\limits_{m:m_{-}=n}\\!\\!(\bm{g}_{m}\,{-}\,R_{m}\bm{l}_{m})\,{+}\,\bm{q}_{n}\,{-}\,B_{n}\bm{v}_{n}\,{=}\,0$
$\displaystyle[\bm{y}^{\textup{q}}_{n}]$ (52b)
$\displaystyle\bm{v}_{n}\,{-}\,2(R_{n}\bm{f}_{n}\,{+}\,X_{n}\bm{g}_{n})\,{+}\,(R_{n}^{2}\,{+}\,X_{n}^{2})\,\bm{l}_{n}\,{=}\,\bm{v}_{n_{-}}$
(52c) $\displaystyle f_{n,t}^{2}+g_{n,t}^{2}\leq v_{n,t}l_{n,t}$
$\displaystyle\forall t\in\mathcal{T}$ (52d) $\displaystyle
f_{n,t}^{2}+g_{n,t}^{2}\leq S_{n}^{2}$ $\displaystyle\forall t\in\mathcal{T}$
(52e) $\displaystyle(f_{n,t}-R_{n}l_{n,t})^{2}+(g_{n,t}-X_{n}l_{n,t})^{2}\leq
S_{n}^{2}$ $\displaystyle\forall t\in\mathcal{T}$ (52f)
$\displaystyle\underline{V}_{n}\leq v_{n,t}\leq\bar{V}_{n},$
$\displaystyle\forall t\in\mathcal{T}$ (52g)

where

  * •

$\bm{v}_{n}=(v_{n,1},\dots,v_{n,T})$ and $\bm{v}_{n_{-}}$ are the squared
voltage magnitudes at buses $n$ and $n_{-}$,

  * •

$\bm{l}_{n}$ is the squared current magnitude on branch $(n,n_{-})$,

  * •

$\bm{f}_{n}$ and $\bm{g}_{n}$ are the active and the reactive parts of the
power flow over line $(n,n_{-})$,

  * •

$R_{n}$ and $X_{n}$ are the resistance and the reactance of branch
$(n,n_{-})$,

  * •

$G_{n}$ and $B_{n}$ are the line conductance and susceptance at $n$.

Equation (52a) and (52b) are the active and reactive flow conservation
equations, (52c) is an expression of Ohm’s law for the branch $(n,n_{-})$, and
(52d) is a SOCP relaxation of the definition of the power flow, which implies
$f_{n,t}^{2}+g_{n,t}^{2}=v_{n,t}l_{n,t}$. There exist sufficient conditions
under which the optimization problem subject to (52) remains exact, see [13,
5, 17]. Equations (52e) and (52f) are limitations on the squared power flow
magnitude on $(n,n_{-})$, and (52g) gives lower and upper bounds on the
voltage at $n$. For the coupling flow conservation laws, dual variables are
attached, which are the DLMPs corresponding to active and reactive power.

For later reference, we point out that the network flow constraints
(52a)-(52b) can be compactly summarised ed as
$\bm{A}_{0}\bm{x}_{0}+\sum_{a}\bm{A}_{a}\bm{x}_{a}=\bm{b}$ for suitably
defined matrices $\bm{A}_{0},\bm{A}_{a}$ and right-hand side vector $\bm{b}$.
However, for our computational scheme to work, we do not need to assume that
the linear constraint holds exactly. Instead, we assume that the linear
relation $\bm{A}_{0}\bm{x}_{0}+\sum_{a}\bm{A}_{a}\bm{x}_{a}=\bm{b}$ holds
inexactly. Specifically, motivated by data-driven approaches for power flow
models [29], we assume that there exists a random variable $\xi$ with bounded
support such that

$\bm{A}_{0}\bm{x}_{0}+\sum_{a}\bm{A}_{a}\bm{x}_{a}=\bm{b}+\xi,\text{ and
}\frac{1}{2}\|\bm{A}_{0}\bm{x}_{0}+\sum_{a}\bm{A}_{a}\bm{x}_{a}-\bm{b}\|^{2}\leq\delta.$

In such inexact regimes, conventional deterministic OPF solvers are not
applicable, whereas Algorithm 3 is designed to handle such scenarios.

###  5.2 Load aggregators

The set of buses $\mathcal{N}$ is partitioned into a collection
$(\mathcal{N}_{a})_{a\in\mathcal{A}}$ of subsets, such that each node subset
$\mathcal{N}_{a}$ is managed by a Load Aggregator $a\in\mathcal{A}$. Each Load
Aggregator (LA) controls the flexible net power consumption ($p_{n,t}$) and
generation at each node $n\in\mathcal{N}_{a}$, given at time $t$ by

$p_{n,t}=p^{\operatorname{c}}_{n,t}-p^{\operatorname{p}}_{n,t},\qquad
q_{n,t}=q^{\operatorname{c}}_{n,t}-q^{\operatorname{p}}_{n,t},$ (53a) for all
$n\in\mathcal{N}$ and $t\in\mathcal{T}.$ $p^{\operatorname{c}}_{n,t}\geq 0$ is
the consumption part and $p^{\operatorname{p}}_{n,t}\geq 0$ is the production
part of the power profile. Power consumption and production at the nodes are
made flexible by the presence of deferrable loads (electric vehicles, water
heaters) and _Distributed Energy Resources_ (DERs). The consumption at each
node $n\in\mathcal{N}$ must satisfy a global energy demand $E_{n}$ over the
full time window, $\sum_{t\in\mathcal{T}}p^{\operatorname{c}}_{n,t}\geq
E_{n},\qquad\forall n\in\mathcal{N}.$ (53b) Consumption and production are
also constrained by power bounds and active to reactive power ratio:
$\displaystyle\underline{P}_{n,t}\leq
p^{\operatorname{c}}_{n,t}\leq\overline{P}_{n,t},$ $\displaystyle\forall
n\in\mathcal{N},\ \forall t\in\mathcal{T},$ (53c) $\displaystyle
q^{\operatorname{c}}_{n,t}=\tau^{\operatorname{c}}_{n}p^{\operatorname{c}}_{n,t},$
$\displaystyle\forall n\in\mathcal{N},\ \forall t\in\mathcal{T},$ (53d)
$\displaystyle 0\leq p^{\operatorname{p}}_{n,t}\leq\overline{P},$
$\displaystyle\forall n\in\mathcal{N},\ \forall t\in\mathcal{T},$ (53e)
$\displaystyle\underline{\rho}^{\operatorname{p}}_{n,t}p^{\operatorname{p}}_{n,t}\leq
q^{\operatorname{p}}_{n,t}\leq\overline{\rho}^{\operatorname{p}}_{n,t}p^{\operatorname{p}}_{n,t},$
$\displaystyle\forall n\in\mathcal{N},\ \forall t\in\mathcal{T}.$ (53f)

Constraints (53a)-(53f) define the feasible set $\mathcal{X}_{a}$ of LA
decisions, containing vectors
$\bm{x}_{a}=(\bm{p}_{n},\bm{q}_{n})_{n\in\mathcal{N}_{a}}$. Both, consumption
and production, must be scheduled by the LA, taking into account the current
spot market prices, and other specific local factors characterising the
private objectives of the LA. Formally, there is a convex cost function
$\phi_{a}(\bm{x}_{a})$ which the LA would like to unilaterally minimise,
subject to private feasibility $\bm{x}_{a}\in\mathcal{X}_{a}$.

###  5.3 The distribution system operator

In order to guarantee stability of the distribution network, the Distribution
System Operator (DSO) takes the individual aggregators’ decisions into account
and adjusts the power flows so that the flow conservation constraints
(52a)-(52b), together with the SOCP constraints (52c)-(52g), are satisfied.
Let $\bm{x}_{0}=(\bm{p}_{0},\bm{q}_{0},\bm{f},\bm{g},\bm{v},\bm{l})$ denote
the vector of the variables controlled by the DSO, and define the DSO’s
feasible set
$\mathcal{X}_{0}=\\{\bm{x}_{0}|\eqref{Ohm}-\eqref{voltagebound}\text{ hold for
}n\in\mathcal{N}\\}.$ Then, the set of DSO decision variables inducing a
physically meaningful network flow for a given tuple of LA decisions
$\bm{x}_{\mathcal{A}}$ is described as
$\mathcal{F}(\bm{x}_{\mathcal{A}})=\\{\bm{x}_{0}\in\mathcal{X}_{0}|\eqref{flowconservationactive}-\eqref{flowconservationreactive}\text{
hold for }\bm{x}_{\mathcal{A}}\\}.$ Denoting the DSO cost function
$\phi_{0}(\bm{x}_{0})$, we arrive at the DSO’s decision problem

$\Psi(\bm{x}_{\mathcal{A}})=\min\\{\phi_{0}(\bm{x}_{0})|\bm{x}_{0}\in\mathcal{F}(\bm{x}_{\mathcal{A}})\\},$

This represents the smallest costs to the DSO, given the profile of flexible
net consumption and generation at each affiliated node $n\in\mathcal{N}_{a}$.

###  5.4 A privacy-preserving DLMP solver

The privacy-preserving DLMP solver (PPDLMP), described in Algorithm 3, asks
the DSO to adjust DLMPs based on the prevailing plans reported by the LAs.
Once the price update is completed, a single LA is selected at random to adapt
the power profile within the subnetwork they manage. The local update of the
LA results in bid vector $w^{k}$, which will be fed into the DSO final
computational step to perform dispatch.

Parameters : $p=\lvert\mathcal{A}\rvert$, $\sigma>0$, $\bm{B}_{0}$,
$\bm{B}_{a}$

Initialization at  _each aggregator $a\in\mathcal{A}$_ :

$\bm{x}_{0}^{0}\in\mathcal{X}_{0}$

send bid $u_{a}=\bm{A}_{a}\bm{x}^{0}_{a}-b_{a}$ to the DSO

1

Initialization at  _the DSO_ :

$\bm{x}^{0}_{0}\in\mathcal{X}_{0},\;\bm{v}^{0}=\sigma\sum_{a\in\mathcal{A}}\bm{u}_{a},\;\bm{y}^{0}=\bm{v}^{0}+\sigma(\bm{A}_{0}\bm{x}^{0}_{0}-b_{0})$

2

3

for  _$k=0,1,2,\dots$_ do

at  _the DSO_ do

4
$\bm{x}^{k+1}_{0}=\arg\min_{\tilde{\bm{x}}_{0}\in\mathcal{X}_{0}}\\{\langle{\nabla\phi_{0}(\bm{x}_{0}^{k})+\bm{A}_{0}^{\top}\bm{y}^{k}},{\tilde{\bm{x}}_{0}}\rangle+\frac{1}{2}\|\tilde{\bm{x}_{0}}-\bm{x}^{k}_{0}\|^{2}_{\bm{B}_{0}}\\}$

5

6

at  _LA $a$ drawn uniformly at random_ do

7 receive DLMP $y^{k}$ from DSO

8
$\bm{x}^{k+1}_{a}=\arg\min_{\tilde{\bm{x}}_{a}\in\mathcal{X}_{a}}\\{\langle{\nabla\phi_{a}(\bm{x}^{k}_{a})+\bm{A}_{a}^{\top}\bm{y}^{k}},{\tilde{\bm{x}}_{a}}\rangle+\frac{p}{2}\|\tilde{\bm{x}}_{a}-\bm{x}^{k}_{a}\|^{2}_{\bm{B}_{a}}\\}$

9 $\bm{s}^{k}=\bm{A}_{a}(\bm{x}^{k+1}_{a}-\bm{x}^{k}_{a})$

10

11

at  _each other aggregator $a^{\prime}\neq a$_ do

12 $\bm{x}^{k+1}_{a^{\prime}}=\bm{x}^{k}_{a^{\prime}}$

13

14

at  _the DSO_ do

15 receive bid $\bm{s}^{k}$ from LA $a$

16
$\bm{y}^{k+1}=\bm{y}^{k}+\sigma[\bm{A}_{0}(2\bm{x}^{k+1}_{0}-\bm{x}^{k}_{0})-b_{0}]+\sigma(p+1)\bm{s}^{k}$

17 $\bm{v}^{k+1}=\bm{v}^{k}+\sigma\bm{s}^{k}$

18

19

20

21

Algorithm 3 privacy-preserving DLMP solver (PPDLMP)

From a practical point of view, it is important to point out that, while
executing PPDLMP, the bus-specific data (like cost function, power profiles,
etc.) remain private information. This applies equally to the DSO and the LA.
Coordination of the system-wide behaviour is achieved via exchanging
information about _dual variables_ only, describing the DLMPs and the
expressed bids of the LAs.

###  5.5 Convergence of PPDLMP

We deduce the convergence of PPDLMP via the analysis of Algorithm 1. In order
to recover the OPF problem, we identify each function $\phi_{i}$ with a cost
function of the DSO or LA, and $r_{i}$ is an indicator function of the
feasible set $\mathcal{X}_{a}$ and $\mathcal{X}_{0}$, respectively. The
sampling technology is an example of a uniform sampling of pairs involving the
DSO and one aggregator. Specifically, the set of realisations of the sampling
is
$\Sigma=\left\\{\\{0,a\\}:a\in\\{1,\ldots,|\mathcal{A}|\\}\equiv\\{1,\ldots,p\\}\right\\}.$
Each pair is realised with the uniform probability $(\Pi)_{0,a}=\frac{1}{p}$.
The marginal distribution of the sampling technology is thus given by
$\pi_{0}=1$ and $\pi_{a}=1/p$ for $a\in\\{1,\ldots,p\\}$. The optimisation
problem is characterised by $d=p+1$ blocks and weighting matrix
$\bm{P}=\operatorname{\texttt{blkdiag}}[\bm{I}_{m_{0}};\ldots;p\bm{I}_{m_{d}}]$,
where $m_{0}$ is the dimension of the feasible set of the DSO, and $m_{a}$ is
the dimension of the feasible set of aggregator $a\in\mathcal{A}$. Now, define
$R(\bm{x})=r_{0}(\bm{x}_{0})+\sum_{a\in\mathcal{A}}r_{a}(\bm{x}_{a})$, where
$r_{0}(\bm{x}_{0})=\delta_{\mathcal{X}_{0}}(\bm{x}_{0})$ and
$r_{a}(\bm{x}_{a})=\delta_{\mathcal{X}_{a}}(\bm{x}_{a})$, in which
$\delta_{C}$ is the indicator function of a set $C$. We now show that
Algorithm 3 is a special instantiation of Algorithm 2, which in turn we know
to be equivalent to Algorithm 1.
If the load aggregator $a$ is chosen at step $k$, Line 5 in Algorithm 2
becomes

$\displaystyle\bm{y}^{k+1}=\bm{y}^{k}+\sigma\bm{A}\bm{P}(\bm{x}^{k+1}-\bm{x}^{k})+\sigma\bm{u}^{k+1}$
$\displaystyle=\bm{y}^{k}+\sigma\bm{A}_{0}(\bm{x}_{0}^{k+1}-\bm{x}^{k}_{0})+\sigma
p\bm{A}_{a}(\bm{x}^{k+1}_{a}-\bm{x}^{k}_{a})+\sigma\bm{u}^{k+1}$
$\displaystyle=\bm{y}^{k}+\sigma[\bm{A}_{0}(2\bm{x}^{k+1}_{0}-\bm{x}^{k}_{0})-b_{0}]+\sigma
p\bm{A}_{a}(\bm{x}^{k+1}_{a}-\bm{x}^{k}_{a})+\bm{v}^{k+1}$

where we have used the identity $\bm{u}^{k+1}=\bm{A}\bm{x}^{k}-\bm{b}$ and
define
$\bm{v}^{k}=\sigma\bm{u}^{k}-\sigma(\bm{A}_{0}\bm{x}^{k}_{0}-b_{0})=\sigma\sum_{a\in\mathcal{A}}(\bm{A}_{a}\bm{x}^{k}_{a}-b_{a})$.
Exploiting Line 4 in Algorithm 2, we find that $\bm{v}^{k}$ can be computed
locally and inductively by choosing the initial condition
$\bm{v}^{0}=\sigma\sum_{a\in\mathcal{A}}(\bm{A}_{a}\bm{x}^{0}_{a}-b_{a})$
initially, and performing sequential updates

$\bm{v}^{k+1}=\bm{v}^{k}+\sigma\bm{A}_{a}(\bm{x}^{k+1}_{a}-\bm{x}^{k}_{a})=\bm{v}^{k}+\sigma\bm{s}^{k}$

where $\bm{s}^{k}=\bm{A}_{a}(\bm{x}^{k+1}_{a}-\bm{x}^{k}_{a})$. This shows

$\bm{y}^{k+1}=\bm{y}^{k}+\sigma[\bm{A}_{0}(2\bm{x}^{k+1}_{0}-\bm{x}^{k}_{0})-b_{0}]+\sigma(p+1)\bm{s}^{k}+\bm{v}^{k},$

which agrees with Lines 4,7 and 8 of Algorithm 3.

###  5.6 Numerical Results

We apply Algorithm 3 to a realistic 15-bus network example based on the
instance proposed in [36], over a time horizon $\mathcal{T}=\\{0,1\\}$. The
parameters $(R_{n},X_{n},S_{n},B_{n},V_{n})$ are those used in [36]. We
consider variable, flexible active and reactive loads (as opposed to fixed
ones, as in [36]): parameters
$(\bm{\underline{P}}_{n},\bm{\overline{P}}_{n},E_{n},\tau^{\mathrm{c}}_{n})_{n}$
are generated based on the values of [36].The underlying parameter values are
specified in Table 1: line parameters $R_{n},X_{n},S_{n},B_{n}$,
$n\in\mathcal{N}$ are taken from [36],
$(\bm{\underline{P}}_{n},\bm{\overline{P}}_{n},E_{n},\tau^{\mathrm{c}}_{n})_{n}$
for the flexible loads are generated based on the _fixed_ load values of [36].
As in [36], bus $11$ is the only bus to offer renewable production, with
$\bm{\overline{P}^{\mathrm{p}}}_{11}:=[0.438,0.201]$ and
$\underline{\rho^{\mathrm{p}}}=\overline{\rho^{\mathrm{p}}}=0$ (the renewable
production is fully active). The bounds $(\underline{V}_{n},\overline{V}_{n})$
are set to 0.81 and $1.21$ for each $n\in\mathcal{N}$, while $V_{0}=1.0$.

We consider a zero cost function for each LA ($\phi_{a}=0$ for each
$a\in\mathcal{A}$), meaning that LAs are indifferent to consumption profiles
for as long as their feasibility constraints are satisfied. This is a
reasonable assumption in practice for certain types of consumption
flexibilities (electric vehicles, batteries). We consider the DSO objective
$\phi(\bm{x})=\phi_{0}(\bm{x}_{0})=\sum_{t\in\mathcal{T}}c_{t}(p^{\mathrm{p}}_{0t})+k^{\text{loss}}\sum_{n,t}R_{n}\ell_{nt},$
with loss penalisation factor $k^{\text{loss}}=0.001$ and with:
$c_{0}:p\mapsto 2p+p^{2},\ c_{1}:p\mapsto p,$ giving an expensive time period
and a cheap one, which can be interpreted as peak and offpeak periods.

$n$ | $S_{n}$ | $R_{n}{\cdot}10^{3}$ | $X_{n}{\cdot}10^{3}$ | $B_{n}{\cdot}10^{3}$ | $\underline{P}_{n}$ | $\overline{P}_{n}$ | $E_{n}$ | $\tau^{c}_{n}$
---|---|---|---|---|---|---|---|---
1 | 2.000 | 1.0 | 120.0 | 1.1 | [0.593, 0.256] | [1.566, 1.539] | 2.213 | 0.234
2 | 0.256 | 88.3 | 126.2 | 2.8 | [0.000, 0.000] | [0.000, 0.000] | 0.000 | 0.000
3 | 0.256 | 138.4 | 197.8 | 2.4 | [0.003, 0.011] | [0.020, 0.035] | 0.047 | 0.418
4 | 0.256 | 19.1 | 27.3 | 0.4 | [0.015, 0.013] | [0.027, 0.019] | 0.033 | 0.249
5 | 0.256 | 17.5 | 25.1 | 0.8 | [0.021, 0.024] | [0.043, 0.053] | 0.072 | 0.251
6 | 0.256 | 48.2 | 68.9 | 0.6 | [0.017, 0.001] | [0.032, 0.037] | 0.039 | 0.251
8 | 0.256 | 40.7 | 58.2 | 1.2 | [0.021, 0.009] | [0.040, 0.039] | 0.049 | 0.251
7 | 0.256 | 52.3 | 74.7 | 0.6 | [-0.233,-0.210] | [-0.173,-0.115] | -0.352 | 0.000
9 | 0.256 | 10.0 | 14.3 | 0.4 | [0.008, 0.002] | [0.032, 0.028] | 0.015 | 0.620
10 | 0.256 | 24.1 | 34.5 | 0.4 | [0.004, 0.001] | [0.024, 0.040] | 0.013 | 0.300
11 | 0.256 | 10.3 | 14.8 | 0.1 | [0.010, 0.010] | [0.015, 0.024] | 0.028 | 0.250
12 | 0.600 | 1.0 | 120.0 | 0.1 | [0.243, 0.057] | [0.642, 0.625] | 0.895 | 0.208
13 | 0.204 | 155.9 | 111.9 | 0.2 | [0.001, 0.000] | [0.003, 0.003] | 0.003 | 0.571
14 | 0.204 | 95.3 | 68.4 | 0.1 | [0.015, 0.012] | [0.032, 0.042] | 0.042 | 0.371
Table 1: Parameters for the 15 buses network based on [36]

The solution obtained by Algorithm 3 after 2000 iterations is illustrated in
Figure 1, which displays the active flows directions as well as the
Distribution Locational Marginal Price (DLMP) values.

012.5cm $3.721$
$0.012$
22.5cm $3.603$
$0.048$
32.5cm $3.421$
$0.092$
42.5cm $3.429$
$0.093$
52.5cm $3.433$
$0.094$
62.5cm $3.439$
$0.096$
82.5cm $0.004$
$0.279$
72.5cm $-0.003$
$0.279$
92.5cm $0.003$
$0.279$
102.5cm $0.001$
$0.279$
112.5cm $-0.0$
$0.279$
122.5cm $3.719$
$0.001$
132.5cm $3.758$
$0.016$
142.5cm $3.781$
$0.024$
(a) t=0

012.5cm $1.004$
$0.005$
22.5cm $0.98$
$0.02$
32.5cm $0.942$
$0.041$
42.5cm $0.946$
$0.042$
52.5cm $0.949$
$0.043$
62.5cm $0.951$
$0.044$
82.5cm $0.003$
$0.115$
72.5cm $-0.0$
$0.115$
92.5cm $0.002$
$0.115$
102.5cm $0.001$
$0.115$
112.5cm $-0.0$
$0.115$
122.5cm $3.567$
$0.732$
132.5cm $3.583$
$0.737$
142.5cm $3.593$
$0.741$
(b) t=1

Figure 1: Directions of active flows $\bm{f}$ and DLMPs $(y^{p},y^{q})$ at the
solution given by Algorithm 3. Saturated lines are dashed.

The solutions show that the active (and reactive) DLMPs obtained for each time
period are close to the DLMPs at the root node
$(\bm{y}^{\textup{p}}_{0},\bm{y}^{\textup{q}}_{0})$, with the following
exceptions:

  * •

For the branch composed of nodes $8,7,9,10,11$, active DLMPs are close to
$0.0$ due to the presence of renewable production (at null cost) at node 11,
and of negative load at node $7$, which together fully compensate for the
demand on this branch. Since Line $(3,8)$ is saturated, no energy can be
exported further.

  * •

Active DLMPs on the branch composed of nodes $(12,13,14)$ at $t=1$ are much
larger than on other nodes: this is explained by the congestion of line
$(0,12)$.

  * •

The DLMP for node 7 and $t=0$ is _strictly negative_ : the (negative)
consumption for this node is at its upper bound
$p_{7,0}=\overline{P}_{7,0}=-0.173$. The negative DLMP suggests that the
system will be better off if less power is injected by node 7.

Figure 2: Convergence of DLMPs $y^{k}$ and the objective function value
$\Phi(\bm{x}^{k})$

Figure 3: Convergence of the average iterate $\bm{w}^{k}$

Convergence of Algorithm 3 for the 15-bus network is shown in Figures 2 and 3.
Figure 2 displays the convergence of the last iterate with respect to various
criteria:

  * (i)

convergence of $\phi(\bm{x}^{k})$ to the optimal cost $\phi(\bm{x}^{\ast})$,
where we computed $\bm{x}^{\ast}$ using the CvxOpt Python library;

  * (ii)

convergence to zero of the primal residuals $h(\bm{x}^{k})$ and
$\|A\bm{x}^{k}-b\|$;

  * (ii)

convergence to zero of the stopping criterion formed on the KKT residual
$R^{\text{KKT}}(\bm{x}^{k},y^{k})=\text{dist}_{\infty}(0,(\partial_{x},-\partial_{y})L(\bm{x}^{k},y^{k}))$,
where $L(\bm{x},y)=\Phi(\bm{x})+\langle{y},{A\bm{x}-b}\rangle$ denotes the
Lagrangian of (P) in the feasible case;

  * (iv)

convergence of the Figure 3 shows the convergence to zero of the primal
infeasiblity in the averaged trajectory $\bm{w}^{k}$, as predicted by the
theory.

##  6 Conclusion

In this paper we developed a novel random block-coordinate descent method for
ill-posed convex non-smooth optimisation problems. Our scheme gives optimal
iteration complexity guarantees in terms of the last iterate of the sequence
generated by the algorithm. As such, we directly generalise results from [26]
to a fully distributed optimisation setting. Motivated by the need to achieve
a distributed optimisation of the power system, we follow a data-driven
approach leading to a potentially inconsistent AC-OPF formulation. We show
that our Algorithmic scheme is immediately applicable to such a challenging
and important setting, and achieves distributed control of the power grid via
distributed locational marginal prices as price signals to independent agents
(i.e. aggregators). Future work should consist in extending the method to non-
convex optimisation problems, so that exact formulation of the power flow
constraints can be implemented. We also plan to conduct extensions of this
work where the electric network is exposed to stochastic uncertainty.

## Acknowledgements

This research benefited from the support of the FMJH Program PGMO and from the
support of EDF. We thank Yura Malitsky for insightful discussions on this
topic and Olivier Bilenne for his initial work on this paper.

##  Appendix A Auxiliary facts

The following relation will be useful in the analysis.

######  Lemma 9 ([2], Corollary 2.14).

For all $\bm{x},\bm{y}\in\mathbb{R}^{d}$ and $t\in[0,1]$, one has

$\|t\bm{x}+(1-t)\bm{y}\|^{2}=t\|\bm{x}\|^{2}+(1-t)\|\bm{y}\|^{2}-t(1-t)\|\bm{x}-\bm{y}\|^{2}.$

Let $(\Omega,\mathcal{F},(\mathcal{F}_{k})_{k\geq 0},\mathbb{P})$ be a
filtered probability space satisfying the usual conditions. We need the
celebrated Robbins-Siegmund Lemma for the convergence analysis [42].

######  Lemma 10.

For every $k\geq 0$, let $\theta_{k},u_{k},\zeta_{k}$ and $t_{k}$ be non-
negative $\mathcal{F}_{k}$-measurable random variables such that
$(\zeta_{k})_{n\geq 0}$ and $(t_{k})_{k\geq 0}$ are summable and for all
$k\geq 0$,

$\mathbb{E}[\theta_{k+1}|\mathcal{F}_{k}]\leq(1+t_{k})\theta_{k}+\zeta_{k}-u_{k}\qquad\mathbb{P}-\text{a.s.}$
(54)

Then $(\theta_{k})_{k\geq 0}$ converges and $(u_{k})_{k\geq 0}$ is summable
$\mathbb{P}$-a.s.

##  Appendix B Step size policy for the accelerated algorithm

Assume that parameter coupling equation (44) holds, where

$\alpha=\lambda_{\max}(\bm{M}^{-1/2}\Xi\bm{M}^{-1/2})^{-1}\text{ and
}\beta=\alpha\lambda_{\max}(\bm{M}^{-1/2}(\Lambda+\Upsilon)\bm{M}^{-1/2}).$

Set $\kappa=\frac{\beta}{\alpha}$. With the choice
$\bm{B}=\bm{P}^{-2}\Upsilon$, the matrix inequality (34) reads as

$[(\alpha-\beta\tau_{k})(\bm{I}+\tau_{k}\bm{P})-\beta\tau^{2}_{k}(\bm{I}-\bm{P})]\tau_{k+1}^{2}+\tau_{k}^{2}[\beta\bm{I}+\alpha(\bm{I}-\bm{P})]\tau_{k+1}-\alpha\tau_{k}^{2}\bm{I}\succeq
0.$

Exploiting the block structure of the involved matrices, we can equivalently
write this condition as a system of quadratic inequalities

$\tau^{2}_{k+1}\left[(\alpha-\beta\tau_{k})(\pi_{i}+\tau_{k})+\beta(1-\pi_{i})\tau^{2}_{k}\right]-\tau_{k+1}\tau^{2}_{k}[\alpha(1-\pi_{i})-\beta\pi_{i}]-\tau^{2}_{k}\alpha\pi_{i}\geq
0,$ (55)

holding for all $i\in[d]$ simultaneously. This defines a quadratic inequality
in $\tau_{k+1}$ of the form
$c_{i,1}(\tau_{k})\tau^{2}_{k+1}+c_{i,2}(\tau_{k})\tau_{k+1}-c_{i,3}(\tau_{k})\geq
0$, with coefficients

$\displaystyle c_{i,1}(y):=(\alpha-\beta
y)(\pi_{i}+y)+\beta(1-\pi_{i})y^{2}=\alpha\pi_{i}\left[1-\kappa
y^{2}-y(\kappa-1/\pi_{i})\right],$ $\displaystyle
c_{i,2}(y):=y^{2}[\beta\pi_{i}-\alpha(1-\pi_{i})]=\alpha\pi_{i}y^{2}(\kappa-1/\pi_{i}+1),$
$\displaystyle c_{i,3}(y):=y^{2}\alpha\pi_{i}$

The structure of the coefficients shows that we can eliminate the parameters
$\alpha$ and $\pi_{i}$, and just continue with the coefficients

$\displaystyle a_{i,1}(y):=1-\kappa
y^{2}-y(\kappa-1/\pi_{i}),\;a_{i,2}(y):=y^{2}(\kappa-1/\pi_{i}+1),\;a_{i,3}(y):=y^{2}.$

Define the function $F:\mathbb{R}^{2}\to\mathbb{R}$ by

$F(x,y):=\max_{i\in[d]}\\{a_{i,1}(y)x^{2}+a_{i,2}(y)x-a_{i,3}(y)\\}\equiv\max_{i\in[d]}F_{i}(x,y).$

This gives us an implicit relation for the step size policy (50) by
$F(\tau_{k+1},\tau_{k})=0$ from which we take the positive solution, if it
exists. To analyse this implicit relation, we proceed as follows.

######  Lemma 11.

The equation $F(x,0)=0$ has the unique solution $x=0$. Furthermore, the
equation $F(x,\alpha/\beta)=0$ has the unique positive solution
$x=\frac{\alpha}{\beta}=\frac{1}{\kappa}$.

###### Proof.

This is a direct numerical computation. ∎

By definition we have
$\kappa=\max_{i\in[d]}\frac{L_{i}+\mu_{i}}{\mu_{i}\pi_{i}}\geq\frac{1}{\pi_{i}}>1$
for all $i\in[d]$ as well. Consequently,

$\delta_{i}:=\kappa-\frac{1}{\pi_{i}}\geq 0\qquad\forall i\in[d],$

and we obtain the concise representation $a_{i,1}(y)=1-\kappa
y^{2}-y\delta_{i}$, as well as, $a_{i,2}(y)=y^{2}(\delta+1)$. It is clear that
$a_{i,1}(y)>0$ for
$0<y<\frac{-\delta_{i}+\sqrt{\delta_{i}^{2}+4\kappa}}{2\kappa}$. Since
$\pi_{i}\in(0,1)$, it can be readily verified that restricting
$y\in(0,1/\kappa)$ is a sufficient condition for $a_{i,1}(y)>0$. Hence, if
$\tau_{k}\in(0,1/\kappa)$, then the characterisation of the step size reduces
to finding the root of the polynomial equation

$\tau^{2}_{k+1}+\frac{a_{2}(\tau_{k})}{a_{1}(\tau_{k})}\tau_{k+1}-\frac{a_{3}(\tau_{k})}{a_{1}(\tau_{k})}=0.$

Since $\frac{a_{2}(y)}{a_{1}(y)}=y^{2}\frac{\delta+1}{1-\kappa
y^{2}-y\delta}$, and $\frac{a_{3}(y)}{a_{1}(y)}=y^{2}\frac{1}{1-\kappa
y^{2}-\delta y}$, we set

$b_{k}:=\frac{\delta+1}{1-\tau_{k}^{2}\kappa-\delta\tau_{k}}\text{ provided
that }\tau_{k}\in(0,1/\kappa),$

so that our polynomial becomes

$\tau^{2}_{k+1}+\tau^{2}_{k}b_{k}\tau_{k+1}-\tau^{2}_{k}\frac{b_{k}}{\delta+1}=0.$
(56)

This has the unique positive solution

$\tau_{k+1}=\frac{\tau_{k}\left[-b_{k}\tau_{k}+\sqrt{b_{k}^{2}\tau^{2}_{k}+4b_{k}/(\delta+1)}\right]}{2}=\frac{2}{\delta+1}\frac{1}{1+\sqrt{1+\frac{4}{b_{k}\tau^{2}_{k}(1+\delta)}}}.$
(57)

Setting $d_{k}:=1-\delta_{k}\tau_{k}-\kappa\tau^{2}_{k}$, and
$s_{k}:=(1+\delta)\tau_{k}$, this recursion can be equivalently stated as

$s_{k+1}=\frac{2s_{k}}{s_{k}+\sqrt{s^{2}_{k}+4d_{k}}}.$

To wit, we note that

$\displaystyle
s_{k+1}=\frac{2}{1+\sqrt{1+\frac{4d_{k}}{\tau^{2}_{k}(1+\delta)^{2}}}}=\frac{2}{1+\sqrt{1+\frac{4d_{k}}{s^{2}_{k}}}}=\frac{2s_{k}}{s_{k}+\sqrt{s^{2}_{k}+4d_{k}}}.$

Hence,
$\frac{s_{k+1}}{s_{k}}=\frac{2}{1+\sqrt{1+\frac{4d_{k}}{s^{2}_{k}}}}\leq 1$
using that $d_{k}>0$. It follows by induction that $s_{k}\leq
s_{0}=(1+\delta)\tau_{0}$. We deduce that the sequence $(\tau_{k})_{k}$ is
monotonically decreasing and non-negative. Therefore,
$\lim_{k\to\infty}\tau_{k}=\tau_{\infty}\geq 0$ exists. We argue that
$\tau_{\infty}=0$. Suppose that $\tau_{\infty}\geq t>0$. Then
$s_{\infty}=\lim_{k\to\infty}s_{k}=(1+\delta)\tau_{\infty}$, and
$\lim_{k\to\infty}(s_{k+1}-s_{k})=0$. But then

$s_{\infty}=\frac{2s_{\infty}}{s_{\infty}+\sqrt{s^{2}_{\infty}+4d_{\infty}}}<s_{\infty},$

a contradiction! We conclude $\tau_{\infty}=0$. If follows that
$\lim_{k\to\infty}b_{k}=1+\delta$, and $\lim_{k\to\infty}d_{k}=1$. Moreover,
$d_{k}\uparrow 1$. Hence,

$s_{k+1}=\frac{2s_{k}}{s_{k}+\sqrt{s^{2}_{k}+4d_{k}}}\geq\frac{2s_{k}}{s_{k}+\sqrt{s^{2}_{k}+4}}=\frac{2}{1+\sqrt{1+4/s_{k}^{2}}},$

which is equivalent to $\frac{1-s_{k+1}}{s^{2}_{k+1}}\leq\frac{1}{s^{2}_{k}}$.
It follows from [45] that $s_{k}\geq\frac{2s_{0}}{ks_{0}+2}$ for all $k\geq
0$, which implies

$\tau_{k}\geq\frac{2\tau_{0}}{(1+\delta)\tau_{0}k+2},\text{ provided
}\tau_{0}\in(0,1/\kappa).$

This implies via our parameter coupling $\sigma_{k}=\alpha/\tau_{k}-\beta$ for
$k\geq 1$, that

$\displaystyle S_{k}$
$\displaystyle=1+\sum_{i=1}^{k}\frac{\alpha}{\tau_{i}}-k\beta\leq
1+\frac{\alpha(1+\delta)}{2}\sum_{i=1}^{k}i+k\left(\frac{\alpha}{\tau_{0}}-\beta\right)$
$\displaystyle=1+\frac{\alpha(1+\delta)}{4}k(k+1)+k\sigma_{0}.$

Hence, $S_{k}=O(k^{2})$.

## References

  * [1] M.E. Baran and F.F. Wu.  Optimal capacitor placement on radial distribution systems.  IEEE Transactions on Power Delivery, 4(1):725–734, 1989.
  * [2] Heinz H. Bauschke and Patrick L. Combettes.  Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  Springer, 2016.
  * [3] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
  * [4] Pascal Bianchi, Walid Hachem, and Franck Iutzeler.  A coordinate descent primal-dual algorithm and application to distributed asynchronous optimization.  IEEE Transactions on Automatic Control, 61(10):2947–2957, 2015\.
  * [5] Daniel Bienstock.  Electrical transmission system cascades and vulnerability: an operations research viewpoint.  MOS-SIAM Series on Optimization, 2015.
  * [6] Olivier Bilenne, Paulin Jacquot, Nadia Oudjane, Mathias Staudigl, and Cheng Wan.  A privacy-preserving distributed computational approach for distributed locational marginal prices.  In 61st IEEE Conference on Decision and Control, 2022.
  * [7] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends®in Machine learning, 3(1):1–122, 2011.
  * [8] Antonin Chambolle and Thomas Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of Mathematical Imaging and Vision, 40(1):120–145, May 2011\.
  * [9] Caihua Chen, Bingsheng He, Yinyu Ye, and Xiaoming Yuan.  The direct extension of admm for multi-block convex minimization problems is not necessarily convergent.  Mathematical Programming, 155(1):57–79, 2016.
  * [10] Alexandre d’Aspremont, Damien Scieur, and Adrien Taylor.  Acceleration methods.  Foundations and Trends in Optimization, 5(1-2):1–245, 2021.
  * [11] Jonathan Eckstein and Wang Yao.  Understanding the convergence of the alternating direction method of multipliers: Theoretical and computational perspectives.  Pac. J. Optim., 11(4):619–644, 2015.
  * [12] Francisco Facchinei and Christian Kanzow.  Generalized nash equilibrium problems.  4OR - A Quarterly Journal of Operations Research, 5:173–210, 2007\.
  * [13] M. Farivar and S. H. Low.  Branch flow model: Relaxations and convexification—part i.  IEEE Transactions on Power Systems, 28(3):2554–2564, 2013.
  * [14] Olivier Fercoq and Pascal Bianchi.  A coordinate-descent primal-dual algorithm with large step size and possibly nonseparable functions.  SIAM Journal on Optimization, 29(1):100–134, 2019.
  * [15] Olivier Fercoq and Peter Richtárik.  Accelerated, parallel, and proximal coordinate descent.  SIAM Journal on Optimization, 25(4):1997–2023, 2015.
  * [16] Daniel Gabay and Bertrand Mercier.  A dual algorithm for the solution of nonlinear variational problems via finite element approximation.  Computers & Mathematics with Applications, 2(1):17–40, 1976.
  * [17] L. Gan, N. Li, U. Topcu, and S. H. Low.  Exact convex relaxation of optimal power flow in radial networks.  IEEE Transactions on Automatic Control, 60(1):72–87, 2015.
  * [18] Xiang Gao, Yang-Yang Xu, and Shu-Zhong Zhang.  Randomized primal–dual proximal block coordinate updates.  Journal of the Operations Research Society of China, 7(2):205–250, 2019.
  * [19] Xiang Gao and Shu-Zhong Zhang.  First-order algorithms for convex optimization with nonseparable objective and coupled constraints.  Journal of the Operations Research Society of China, 5(2):131–159, 2017.
  * [20] Roland Glowinski and Americo Marroco.  Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualitéd’une classe de problèmes de dirichlet non linéaires.  Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique, 9(R2):41–76, 1975.
  * [21] Markus Grasmair, Markus Haltmeier, and Otmar Scherzer.  The residual method for regularizing ill-posed problems.  Applied Mathematics and Computation, 218(6):2693–2710, 2011.
  * [22] Deren Han, Defeng Sun, and Liwei Zhang.  Linear rate convergence of the alternating direction method of multipliers for convex composite programming.  Mathematics of Operations Research, 43(2):622–637, 2022/08/10 2017\.
  * [23] Myeongmin Kang, Myungjoo Kang, and Miyoun Jung.  Inexact accelerated augmented lagrangian methods.  Computational Optimization and Applications, 62(2):373–404, 2015\.
  * [24] Puya Latafat, Nikolaos M Freris, and Panagiotis Patrinos.  A new randomized block-coordinate primal-dual proximal algorithm for distributed optimization.  IEEE Transactions on Automatic Control, 64(10):4050–4065, 2019\.
  * [25] Zhaosong Lu and Lin Xiao.  On the complexity analysis of randomized block-coordinate descent methods.  Mathematical Programming, 152(1):615–642, 2015.
  * [26] D. Russell Luke and Yura Malitsky.  Block-Coordinate Primal-Dual Method for Nonsmooth Minimization over Linear Constraints, pages 121–147.  Springer International Publishing, Cham, 2018.
  * [27] Yura Malitsky.  Chambolle-pock and tseng’s methods: relationship and extension to the bilevel optimization.  06 2017.
  * [28] Yura Malitsky.  The primal-dual hybrid gradient method reduces to a primal method for linearly constrained optimization problems, 2019.
  * [29] Ilyes Mezghani, Sidhant Misra, and Deepjyoti Deka.  Stochastic ac optimal power flow: A data-driven approach.  Electric Power Systems Research, 189:106567, 2020.
  * [30] Daniel K. Molzahn and Ian A. Hiskens.  A survey of relaxations and approximations of the power flow equations.  Foundations and Trends® in Electric Energy Systems, 4(1-2):1–221, 2019.
  * [31] Ion Necoara and Dragos Clipici.  Efficient parallel coordinate descent algorithm for convex optimization problems with separable constraints: application to distributed mpc.  Journal of Process Control, 23(3):243–253, 2013.
  * [32] Yu. Nesterov.  A method of solving a convex programming problem with convergence rate $o(1/k^{2})$.  Soviet Mathematics Doklady, 27(2):372–376, 1983.
  * [33] Yu Nesterov.  Efficiency of coordinate descent methods on huge-scale optimization problems.  SIAM Journal on Optimization, 22(2):341–362, 2012.
  * [34] Yu. Nesterov.  Gradient methods for minimizing composite functions.  Mathematical Programming, 140(1):125–161, 2013.
  * [35] Yurii Nesterov.  Lectures on Convex Optimization, volume 137 of Springer Optimization and Its Applications.  Springer International Publishing, 2018.
  * [36] A. Papavasiliou.  Analysis of distribution locational marginal prices.  IEEE Transactions on Smart Grid, 9(5):4872–4882, 2018.
  * [37] Q. Peng and S. H. Low.  Distributed optimal power flow algorithm for radial networks, i: Balanced single phase case.  IEEE Transactions on Smart Grid, 9(1):111–121, 2018.
  * [38] Zheng Qu and Peter Richtárik.  Coordinate descent with arbitrary sampling i: algorithms and complexity.  Optimization Methods and Software, 31(5):829–857, 09 2016.
  * [39] Zheng Qu and Peter Richtárik.  Coordinate descent with arbitrary sampling ii: expected separable overapproximation.  Optimization Methods and Software, 31(5):858–884, 09 2016.
  * [40] Peter Richtárik and Martin Takáč.  Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function.  Mathematical Programming, 144(1):1–38, 2014.
  * [41] Peter Richtárik and Martin Takáč.  Parallel coordinate descent methods for big data optimization.  Mathematical Programming, 156(1436-4646):433–484, 2016.
  * [42] H. Robbins and D. Siegmund.  A convergence theorem for non negative almost supermartingales and some applications, pages 233–257.  Academic Press, 1971.
  * [43] Quoc Tran-Dinh and Deyi Liu.  Faster randomized primal-dual algorithms for nonsmooth composite convex minimization.  arXiv preprint arXiv:2003.01322, 2020.
  * [44] Quoc Tran-Dinh and Deyi Liu.  A new randomized primal-dual algorithm for convex optimization with fast last iterate convergence rates.  Optimization Methods and Software, pages 1–34, 09 2022.
  * [45] P. Tseng.  On accelerated proximal gradient methods for convex-concave optimization.  2008\.
  * [46] Allen J Wood, Bruce F Wollenberg, and Gerald B Sheblé.  Power generation, operation, and control.  John Wiley & Sons, 2013.
  * [47] Yangyang Xu and Shuzhong Zhang.  Accelerated primal-dual proximal block coordinate updating methods for constrained convex optimization.  Computational Optimization and Applications, 70(1):91–128, 2018\.

  *[LA]: Load Aggregator
  *[DERs]: Distributed Energy Resource
  *[DSO]: Distribution System Operator
  *[PPDLMP]: privacy-preserving DLMP solver
  *[DLMPs]: Distribution Locational Marginal Price
  *[LAs]: Load Aggregator
  *[ DSO]: Distribution System Operator
  *[DLMP]: Distribution Locational Marginal Price