Proximal Splitting Meets Variance Reduction

Fabian Pedregosa Kilian Fatras Mattia Casotto

ETH Zürich and UC Berkeley111Currently at Google AI, Canada USA and
Switzerland Univ. Bretagne-Sud, CNRS, IRISA Inria Rennes and Obelix Akur8
France

###### Abstract

Despite the raise to fame of stochastic variance reduced methods like Saga and
ProxSvrg, their use in non-smooth optimization is still limited to a few
simple cases. Existing methods require to compute the proximal operator of the
non-smooth term at each iteration, which, for complex penalties like the total
variation, overlapping group lasso or trend filtering, is an iterative process
that becomes unfeasible for moderately large problems. In this work we propose
and analyze Vr-Tos, a variance-reduced method to solve problems with an
arbitrary number of non-smooth terms. Like other variance reduced methods, it
only requires to evaluate one gradient per iteration and converges with a
constant step size, and so is ideally suited for large scale applications.
Unlike existing variance reduced methods, it admits multiple non-smooth terms
whose proximal operator only needs to be evaluated once per iteration. We
provide a convergence rate analysis for the proposed methods that achieves the
same asymptotic rate as their full gradient variants and illustrate its
computational advantage on 4 different large scale datasets.

##  1 Introduction

Stochastic variance reduced methods (Le Roux et al.,, 2012; Johnson and
Zhang,, 2013; Shalev-Shwartz and Zhang,, 2013) have been recently proposed as
an improved alternative to the venerable stochastic gradient descent (Sgd)
method (Robbins and Monro,, 1951). As Sgd, these methods only require to visit
a small batch of random examples per iteration. This makes them ideally suited
for large scale machine learning problems. Unlike Sgd, the variance of the
updates decreases to zero –hence the name– and converge with non-decreasing
step sizes.

While initial stochastic variance reduced methods only considered smooth
objectives, variants with support for a non-smooth term like ProxSvrg (Xiao
and Zhang,, 2014) and Saga (Defazio et al.,, 2014) were soon developed. These
methods are highly efficient whenever the nonsmooth part is _proximal_ , that
is, its proximal operator is available in closed form or at least fast to
compute. This includes penalties such as the $\ell_{1}$ or group lasso norm,
but not more complex ones like the overlapping group lasso (Jacob et al.,,
2009), multidimensional total variation (Barbero and Sra,, 2014) or trend
filtering (Kim et al.,, 2009), to name a few.

A key observation is that many of these complex penalties can be decomposed as
a sum of proximal terms. Proximal splitting methods like the three operator
splitting (Davis and Yin,, 2017) or the Condat-Vũ algorithm (Condat, 2013b, ;
Vũ,, 2013) then provide a principled approach to incorporate these penalties
into the optimization. However, these methods require to compute the full
gradient of the smooth term at each iteration, which can become costly in the
context of large scale machine learning problems as it involves a full pass
over the dataset. A question of key practical interest is whether these
proximal splitting methods can be accelerated through the use of stochastic
variance reduction techniques.

Our main contribution is the development and analysis of Vr-Tos, a stochastic
variance reduced method that can solve problems with a sum of proximal terms.

The proposed method bridges two previously distant families of algorithms and
inherit the best of both: like the three operator splitting of Davis and Yin,
(2017), it can solve problems with multiple proximal terms, and like variance
reduced stochastic methods its cost is independent on the number of smooth
terms, converges with a fixed step size, and reaches the same asymptotic
convergence rate than full gradient methods. Furthermore, we also develop a
sparse variant of the proposed algorithm which can take advantage of the
sparsity in the input data. The paper is organized as follows:

  * •

_Method._ §2 describes the Vr-Tos algorithm, and extends it in §2.1 to
leverage sparsity in the input data. §2.2 extends these methods to the case of
an arbitrary number of proximal terms.

  * •

_Analysis._ In §4 we provide a non-asymptotic convergence analysis of the
proposed method. We show that, like other variance reduced methods, it
converges with a fixed step size and can achieve the same asymptotic rate as
the full gradient variants.

  * •

_Experiments._ In §5 we compare the proposed method and related algorithms on
a logistic regression problem with overlapping group lasso penalty on 4
datasets.

###  1.1 Definitions and notation

By convention, we denote vectors and vector-valued functions in lowercase
boldface (e.g. ${\boldsymbol{x}}$) and matrices in uppercase boldface letters
(e.g. $\boldsymbol{D}$). The proximal operator of a convex lower
semicontinuous function $h$ is defined as
$\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname*{{arg\,min}}_{{\boldsymbol{z}}\in{\mathbb{R}}^{p}}\\{h({\boldsymbol{z}})+\frac{1}{2\gamma}\|{\boldsymbol{x}}-{\boldsymbol{z}}\|^{2}\\}$.
We say a function $f$ is $L$-smooth if it is differentiable and its gradient
is $L$-Lipschitz, while it is $\mu$-strongly convex if
$f-\frac{\mu}{2}\|\cdot\|^{2}$ is convex.

We denote by $[\,{\boldsymbol{x}}\,]_{b}$ the $b$-th coordinate in
${\boldsymbol{x}}$. This notation is overloaded so that for a collection of
blocks $T=\\{B_{1},B_{2},\ldots\\}$, $[{\boldsymbol{x}}]_{T}$ denotes the
vector ${\boldsymbol{x}}$ restricted to the coordinates in the blocks of $T$.
For convenience, when $T$ consists of a single block $B$ we use
$[{\boldsymbol{x}}]_{B}$ as a shortcut of $[{\boldsymbol{x}}]_{\\{B\\}}$.
Finally, we distinguish ${\mathbb{E}}$, the full expectation taken with
respect to all the randomness in the system, from $\mathbf{E}$, the
conditional expectation with respect to the random index sampled at iteration
$t$, conditioned on all randomness up to iteration $t$.

##  2 Methods

In this section we present our main contribution, the variance reduced three
operator splitting method. We will first consider problems with only two non-
smooth terms, and generalize this formulation to an arbitrary number in §2.2.

We consider the following optimization problem

$\displaystyle\vphantom{\sum^{i}}\operatorname*{minimize}_{{\boldsymbol{x}}\in{\mathbb{R}}^{p}}\,f({\boldsymbol{x}})+g({\boldsymbol{x}})+h({\boldsymbol{x}})\,,$
(OPT) $\displaystyle\vphantom{\sum_{i}^{n}}\text{ with
}f({\boldsymbol{x}})=\textstyle\frac{1}{n}\sum_{i=1}^{n}\psi_{i}({\boldsymbol{x}})+\omega({\boldsymbol{x}})~{}$

where each $\psi_{i}$ is convex and $L_{\psi}$-smooth, $\omega$ is convex and
$L_{\omega}$-smooth and $g,h$ are _proximal_ , i.e., convex and we have access
to their proximal operator.

This formulation allows to express a broad range of problems arising in
machine learning and signal processing: the finite-sum includes common loss
functions such as least squares or logistic loss; the two proximal terms $g,h$
can be extended to an arbitrary number and include penalties such as the group
lasso with overlap, total variation, $\ell_{1}$ trend filtering, etc.
Furthermore, the proximal terms can be extended-valued, thus allowing for
convex constraints through the use of the indicator function. With respect to
previous work, this significantly enlarges the class of functions stochastic
variance reduced methods can solve efficiently.

We allow the terms inside the finite sum to be an addition of two terms:
$\psi_{i}$ and $\omega$. This might seem superfluous since it is not more
general than the standard formulation with a single term. However, in practice
$\psi_{i}$ (e.g., a least squares or logistic loss, see Appendix F.1) can be
highly structured and allow for reduced storage schemes and/or have sparse
gradients (see §2.1), properties which might not be shared by $\omega$, (e.g.,
an $\ell_{2}$ regularization term).

Input: ${\boldsymbol{y}}_{0}\in{\mathbb{R}}^{p}$,
${\boldsymbol{\alpha}}_{0}\in{\mathbb{R}}^{n\times p}$, $\gamma>0$

1

2Temporary storage: ${\boldsymbol{z}}_{t}$, ${\boldsymbol{v}}_{t}$ and
${\boldsymbol{x}}_{t}$, all in ${\mathbb{R}}^{p}$

Result: approximate solution to (OPT)

3

4for  _$t=0,1,\ldots$ _ do

5

6 ${\boldsymbol{z}}_{t}=\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{y}}_{t})$

7 Sample $i\in\\{1,\ldots,n\\}$ uniformly at random

8
${\boldsymbol{v}}_{t}=\nabla\psi_{i}({\boldsymbol{z}}_{t})-{\boldsymbol{\alpha}}_{i,t}+\overline{{\boldsymbol{\alpha}}}_{t}+\nabla\omega({\boldsymbol{z}}_{t})$

9 ${\boldsymbol{x}}_{t}=\operatorname{\mathbf{prox}}_{\gamma
g}(2{\boldsymbol{z}}_{t}-{\boldsymbol{y}}_{t}-\gamma{\boldsymbol{v}}_{t})$

10
${\boldsymbol{y}}_{t+1}={\boldsymbol{y}}_{t}+{\boldsymbol{x}}_{t}-{\boldsymbol{z}}_{t}$

11 Update ${\boldsymbol{\alpha}}_{t+1}$ according to (1)

12

13return ${\boldsymbol{z}}_{t}$

Algorithm 1 Variance Reduced Tos (Vr-Tos)

Central to our algorithm is the concept of $q$-memorization (Hofmann et al.,,
2015), which we recall below. It provides a convenient abstraction over common
gradient memorization techniques like the ones in Saga and Svrg.

######  Definition 1.

A _uniform $q$-memorization_ algorithm selects at each iteration $t$ a random
index set $J_{t}$ of memory terms to update according to

${\boldsymbol{\alpha}}_{j,t+1}=\begin{cases}\nabla
f_{j}({\boldsymbol{z}}_{t})&\text{ if $j\in J_{t}$}\\\
{\boldsymbol{\alpha}}_{j,t}&\text{ otherwise\,,}\end{cases}$ (1)

such that any $j$ has the same probability $q/n$ of being updated.

We now introduce the variance-reduced three operator splitting (Vr-Tos), a
method to solve problems of the form (OPT). It is specified in Algorithm 1 and
takes as input a vector of coefficients
${\boldsymbol{y}}_{0}\in{\mathbb{R}}^{p}$, a table
${\boldsymbol{\alpha}}_{0}\in{\mathbb{R}}^{n\times p}$ to store previous
gradients and a step size $\gamma>0$. Although in the general case this table
is required to be of size $n\times p$, for linearly-parametrized loss
functions like the logistic or least squares loss this can be reduced to size
$n$ (Appendix F.1). Furthermore, the Svrg-like update detailed below avoids
the need for this storage at the expense of a lightly increased per iteration
cost.

The proposed method performs one evaluation of each of the proximal terms and
builds the gradient estimator ${\boldsymbol{v}}_{t}$ from the table of
previous gradients ${\boldsymbol{\alpha}}_{t}$ and the index $i$ sampled
uniformly at random. It is easy to see that ${\boldsymbol{v}}_{t}$ is an
unbiased estimate of the gradient, that is,
$\mathbf{E}\,{\boldsymbol{v}}_{t}=\nabla f({\boldsymbol{z}}_{t})$.

This method allows the memory terms to be updated using any scheme that
verifies the $q$-memorization framework (line 1). Some common schemes are:

  * •

_Saga -like update_. At each iteration, the algorithm updates the same
coefficient that has been sampled, i.e. $J_{t}=\\{i\\}$. In this scheme each
memory term has probability $1/n$ of being updated, and so $q=1$.

  * •

_Svrg -like update_. Fix parameter $q>0$ and draw at each iteration $r$ from a
uniform distribution in the $[0,1]$ interval. If $r<q/n$, the algorithm
performs a complete update
${\boldsymbol{\alpha}}_{j,t+1}=\nabla\psi_{j}({\boldsymbol{z}}_{t})$ for all
$j$, otherwise they are left unchanged.

Like in the Svrg algorithm (Johnson and Zhang,, 2013), it is possible to avoid
storing the memory terms since the $\overline{\boldsymbol{\alpha}}_{t}$ is
constant unless a full refresh is triggered. In this setting, only the
$p$-dimensional vectors $\overline{\boldsymbol{\alpha}}_{t}$ and
$\widetilde{{\boldsymbol{z}}}_{t}$ needs to be stored, where
$\widetilde{{\boldsymbol{z}}}_{t}$ is the value of ${\boldsymbol{z}}_{t}$ last
time a full refresh was triggered. This variant avoids the need to store
${\boldsymbol{\alpha}}_{t}$, at the cost of a slight per iteration cost, as
${\boldsymbol{\alpha}}_{i}=\nabla f_{i}(\widetilde{{\boldsymbol{z}}}_{t})$
needs to be computed at each iteration.

This memory update scheme was proposed by Hofmann et al., (2015), and unlike
the original Svrg algorithm the number of iterates between two full regresh is
a random variable instead of a fixed number of iterations.

###  2.1 Sparse Vr-Tos

#### Need for a sparse variant.

Modern web-scale optimization problems that arise in machine learning are not
only large, they are also often _sparse_. For example, in the LibSVM datasets
suite222https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/, 8 out of the
11 datasets with more than a million samples have a density below $0.01\%$,
and the largest one in number of samples has a density below 1 per million.
Linearly-parametrized loss functions of the form
$\psi_{i}({\boldsymbol{x}})=l_{i}({\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}})$
have a gradient of the form
$\nabla\psi_{i}({\boldsymbol{x}})={\boldsymbol{a}}_{i}l_{i}^{\prime}({\boldsymbol{a}}_{i}^{T}{\boldsymbol{x}})$,
which inherits the same sparsity pattern as the data ${\boldsymbol{a}}_{i}$.
Since the data might be extremely sparse, it is hence of great practical
interest to leverage sparsity in the partial gradients. This is the case in
generalized linear models such as least squares or logistic regression, where
${\boldsymbol{a}}_{i}$ are the rows of a data matrix.

In this subsection we assume that $g$ and $\omega$ are block separable, i.e.,
can be decomposed block coordinate-wise as
$g({\boldsymbol{x}})=\textstyle{\sum_{B\in\mathcal{B}}}g_{B}([{\boldsymbol{x}}]_{B})$
and
$\omega({\boldsymbol{x}})=\textstyle{\sum_{B\in\mathcal{B}}}\omega_{B}([{\boldsymbol{x}}]_{B})$,
where $\mathcal{B}$ is a partition of the coefficients into subsets which will
call _blocks_ and $g_{B},\omega_{B}$ only depends on coordinates in block $B$.
Furthermore, we will make use of the following notation:

  * •

_Extended support_. We define the extended support of $\nabla\psi_{i}$,
denoted $T_{i}$ as the set of blocks of $\mathcal{B}$ that intersect with its
support, formally defined as
$T_{i}\stackrel{{\scriptstyle\text{def}}}{{=}}\\{B:\text{supp}(\nabla
f_{i})\cap B\neq\varnothing,\,B\in\mathcal{B}\\}$. For totally separable
penalties such as the $\ell_{1}$ norm, the blocks are individual coordinates
and so the extended support covers the same coordinates as the support.

  * •

_Reweighting constants_. Let $\boldsymbol{P}_{i}$ be the projection onto the
extended support, i.e., the diagonal matrix where
$[{\boldsymbol{P}}_{i}]_{B,B}$ is the identity if $B\in T_{i}$ and zero
otherwise. For simplicity we assume that each block appears in at least one
$T_{i}$, as otherwise the problem can be reformulated without it. For each
block $B\in\mathcal{B}$ we define $d_{B}$ as the inverse frequency of that
block in the extended support, i.e.
$d_{B}=1/(\frac{1}{n}\sum_{i=1}^{n}\mathds{1}\\{B\in T_{i}\\})^{-1}$. For
notational convenience we define the block-diagonal matrix ${\boldsymbol{D}}$
as $[{\boldsymbol{D}}]_{B,B}=d_{B}\boldsymbol{I}$ for each block
$B\in\mathcal{B}$. Note that by definition
$\frac{1}{n}\sum_{i=1}^{n}{\boldsymbol{P}}_{i}={\boldsymbol{D}}^{-1}$.
Computation of this diagonal matrix should be done as a preprocessing step of
the algorithm.

Input: ${\boldsymbol{y}}_{0}\in{\mathbb{R}}^{p}$,
${\boldsymbol{\alpha}}_{0}\in{\mathbb{R}}^{n\times p}$, $\gamma>0$

1

2Temporary storage: ${\boldsymbol{z}}_{t}$, ${\boldsymbol{v}}_{t}$ and
${\boldsymbol{x}}_{t}$, all in ${\mathbb{R}}^{p}$

Result: approximate solution to (OPT)

3

4for  _$t=0,1,\ldots$ _ do

5

6 Sample $i\in\\{1,\ldots,n\\}$ uniformly at random

7 $T_{i}=\text{extended support of }\nabla\psi_{i}$

8
$[{\boldsymbol{z}}_{t}]_{T_{i}}=[\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}}_{t})]_{T_{i}}$

9
$[{\boldsymbol{v}}_{t}]_{T_{i}}=[\nabla\psi_{i}({\boldsymbol{z}}_{t})\\!-\\!{\boldsymbol{\alpha}}_{i,t}+{\boldsymbol{D}}(\overline{{\boldsymbol{\alpha}}}_{t}+\nabla\omega({\boldsymbol{z}}_{t}))]_{T_{i}}$

10
$[{\boldsymbol{x}}_{t}]_{T_{i}}=[\operatorname{\mathbf{prox}}_{\gamma\varphi_{i}}(2{\boldsymbol{z}}_{t}-{\boldsymbol{y}}_{t}-\gamma{\boldsymbol{v}}_{t})]_{T_{i}}$

11
$[{\boldsymbol{y}}_{t+1}]_{T_{i}}=[{\boldsymbol{y}}_{t}+{\boldsymbol{x}}_{t}-{\boldsymbol{z}}_{t}]_{T_{i}}$

12 Update ${\boldsymbol{\alpha}}_{t+1}$ according to (1)

13

14return $\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}}_{t})$

Algorithm 2 Sparse Vr-Tos

  * •

The _scaled proximal operator_ is defined for a function $\varphi$, step size
$\gamma>0$, positive definite matrix $\boldsymbol{H}$ and norm
$\|\cdot\|_{{\boldsymbol{H}}}^{2}\stackrel{{\scriptstyle\text{def}}}{{=}}\langle\cdot,{\boldsymbol{H}}\cdot\rangle$
as

$\operatorname{\mathbf{prox}}^{\boldsymbol{H}}_{\gamma\varphi}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname*{{arg\,min}}_{{\boldsymbol{z}}\in{\mathbb{R}}^{p}}\big{\\{}\,\varphi({\boldsymbol{z}})+\frac{1}{2\gamma}\|{\boldsymbol{x}}-{\boldsymbol{z}}\|_{\boldsymbol{H}}^{2}\,\big{\\}}$
(2)

We now have all necessary ingredients to present the sparse variant of Vr-Tos.
This is specified in Algorithm 2. In this variant, all operations are
restricted to the extended support.

The algorithm requires to compute the scaled proximal operators of $g$ and
$h$. By block separability of $g$ its scaled proximal operator can be computed
in block-wise as
$[\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
g}({\boldsymbol{x}})]_{B}=[\operatorname{\mathbf{prox}}_{(d_{B}\gamma)h}({\boldsymbol{x}})]_{B}$
for all $B\in\mathcal{B}$. Hence the cost of computing
$[{\boldsymbol{x}}_{t}]_{T_{i}}$ will depend on the extended support size and
not on the dimensionality.

We can unfortunately not guarantee the same complexity for
$[{\boldsymbol{z}}_{t}]_{T_{i}}$ since we do not have a closed form for the
scaled proximal operator of $h$ in general. We review some specific cases in
which it is possible to compute this scaled proximal operator in Appendix D.
Alternatively, in the next subsection we propose a reformulation that avoids
the need to compute this scaled proximal operator at the expense of higher
memory usage.

In the case that one proximal term is zero, the proposed algorithm with Saga-
like update of the memory terms defaults to the Sparse Saga variant of
Pedregosa et al., (2017). With Svrg-like update of the memory terms it instead
yields a novel sparse variant of ProxSvrg (Xiao and Zhang,, 2014). For both of
the proposed algorithms, when input is dense,
${\boldsymbol{P}}_{i}={\boldsymbol{D}}=\boldsymbol{I}$ and we recover
Algorithm 1.

###  2.2 Extension to an arbitrary number of proximal terms

The proposed method can be easily extended to the more general setting of an
objective function with an arbitrary number of proximal terms of the form

$\displaystyle\operatorname*{minimize}_{{\boldsymbol{x}}\in{\mathbb{R}}^{p}}\,f({\boldsymbol{x}})+\textstyle\sum_{j=1}^{k}g_{j}({\boldsymbol{x}})\,,$
$\displaystyle\text{ with
}f({\boldsymbol{x}})=\textstyle\frac{1}{n}\sum_{i=1}^{n}\psi_{i}({\boldsymbol{x}})+\omega({\boldsymbol{x}})~{},$

where $\psi_{i}$ and $\omega$ are as in (OPT) and $g_{1},\ldots,g_{k}$ are
proximal. This is done by expressing the above as a problem of the form (OPT)
in an enlarged space and then applying the proposed algorithm to this
reformulation. For this, we will introduce $k$ new variables which we will
constrain to be equal via an indicator function. The above problem can be
written equivalently as follows,

$\min_{{\boldsymbol{X}}\in{\mathbb{R}}^{k\times
p}}\,f(\overline{\boldsymbol{X}})+\underbrace{\textstyle\sum_{j=1}^{k}g_{j}({\boldsymbol{X}}_{j})}_{\stackrel{{\scriptstyle\text{def}}}{{=}}g({\boldsymbol{X}})}+\underbrace{\imath\\{{\boldsymbol{X}}_{1}\\!=\\!\cdots\\!=\\!{\boldsymbol{X}}_{k}\\}}_{\stackrel{{\scriptstyle\text{def}}}{{=}}h({\boldsymbol{X}})}\,,$

where we have split the original variable into $k$ variables
${\boldsymbol{X}}_{1},\ldots,{\boldsymbol{X}}_{k}$ and constrained them to be
equal using an indicator function in the last term. In this formulation the
first term is smooth, and the other two terms are proximal. The second term is
proximal since the variables in $g_{i}$ are decoupled, each $g_{i}$ is
proximal by assumption and the last term is an indicator function over a
linear subspace, and hence its scaled proximal operator can be computed in
closed form as follows (Lemma 15):

$\displaystyle[\operatorname{\mathbf{prox}}^{\boldsymbol{D}^{-1}}_{\gamma
h}({\boldsymbol{X}})]_{i,j}=\left({\textstyle\sum_{i=1}^{n}}a_{i,j}{\boldsymbol{X}}_{i,j}\right)/\left({\textstyle\sum_{i=1}^{n}}a_{i,j}\right)$
$\displaystyle\text{ with }a_{i,j}={\boldsymbol{D}}^{-1}_{ip+j,ip+j}~{},$ (3)

Hence, the problem with multiple proximal terms (2.2) can be formulated as a
problem with two proximal terms (OPT) and so it is possible to apply the
proposed method defined in the previous subsections. This gives a variance
reduced method for problems with an arbitrary number of proximal term. It is
worth noting that for the sparse variants this formulation avoids the
potentially difficult computation of the scaled proximal operator of $h$.

##  3 Related work

| Methods | incremental | non-decreasing | multiple non-smooth | sparse updates
---|---|---|---|---|---
| updates | step size | terms
| Vr-Tos |  |  |  |
| (_this work_) | ✓ | ✓ | ✓ | ✓
Saga |  |  |  |
(Defazio et al.,, 2014) | ✓ | ✓ | ✗ | ✓(Pedregosa et al.,, 2017)
ProxSvrg |  |  |  |
(Xiao and Zhang,, 2014) | ✓ | ✓ | ✗ | ✗ ${\dagger}$
Tos |  |  |  |
(Davis and Yin,, 2017) | ✗ | ✓ | ✓ | N/A
Stochastic Tos |  |  |  |
(Yurtsever et al.,, 2016) |  ✓ | ✗ | ✓ | ✗
Table 1: Comparison with related work. The proposed method is unique in that
it combines the advantages of variance-reduced methods (incremental updates,
non-decreasing step sizes and sparse updates) with the advantages of proximal
splitting (support for multiple non-smooth terms). $\dagger$: a sparse variant
of ProxSvrg follows as a special case of Algorithm 2 with $h=0$ and the SVRG-
like update of the memory terms.

We comment on the most closely related ideas, summarized in Table 1.

Methods that support objective functions of the form (OPT) with two or more
proximal terms and a smooth term accessed via its gradient have recently been
proposed. Examples are the the primal-dual hybrid gradient method (also known
as the Condat-Vũ) (Condat, 2013a, ; Vũ,, 2013),333We note that this method can
optimize the more general objective function
$f({\boldsymbol{x}})+g({\boldsymbol{x}})+h(\boldsymbol{L}{\boldsymbol{x}})$,
for an arbitrary linear operator $\boldsymbol{L}$ that is fixed to the
identity in our setting. the generalized forward-backward splitting (Raguet et
al.,, 2013) or the three operator splitting (Davis and Yin,, 2017). Due to its
excellent empirical performance and amenability to sparse updates we have
chosen this last method as the basis for the proposed method. The proposed Vr-
Tos method can be seen as a generalization of this last method, as both method
are identical when $n=1$.

A different stochastic variant of the three operator splitting was proposed by
Yurtsever et al., (2016) for the slightly more general case in which $f$ is
given by an expectation. Like the proposed algorithms, this method only needs
to evaluate the gradient of one element in the finite sum per iteration.
Unlike the proposed methods, the variance of the updates does not decrease to
zero and requires –as other non-variance reduced method– a decreasing step
size. Furthermore, all updates are dense even in the presence of sparse
gradients so the method performs poorly on large sparse problems.

(Balamurugan and Bach,, 2016) proposed a variance-reduced method to solve
problems a general class of saddle point problems including
$\min_{\boldsymbol{x}}\max_{\boldsymbol{u}}\frac{1}{n}\sum_{i=1}^{n}f_{i}({\boldsymbol{x}})+M({\boldsymbol{x}},{\boldsymbol{u}})$,
where $M(\cdot)$ is proximal. With
$M({\boldsymbol{x}},{\boldsymbol{u}})=g({\boldsymbol{x}})+\langle{\boldsymbol{x}},{\boldsymbol{u}}\rangle-h^{*}({\boldsymbol{u}})$,
this is equivalent to the problem in (OPT). However, the method requires $M$
to be strongly concave in ${\boldsymbol{u}}$, which is equivalent to $h$ being
smooth, and so is not applicable to the same class of problems as the proposed
method. We note that this requirement is not merely an artifact of the theory,
as the algorithm requires knowledge of this smoothness parameter.

Stochastic variance-reduced variants of ADMM have also been recently proposed,
see e.g. (Zheng and Kwok,, 2016; Yu and Huang,, 2017). Compared to the
proposed methods, none of the existing variants support sparse updates and
require tuning more than one step-size parameter.

##  4 Analysis

In this section we provide a non-asymptotic convergence rate analysis for the
proposed method:

  * •

All the proposed variants converge with a step size $1/(3L_{f})$, with
$L_{f}\stackrel{{\scriptstyle\text{def}}}{{=}}L_{\psi}+d_{\max}L_{\omega}$,
where $d_{\max}$ is the maximum element in the diagonal matrix
${\boldsymbol{D}}$ ($d_{\max}=1$ for non-sparse variants).

  * •

For Vr-Tos (Algorithm 1) we obtain convergence rates that asymptotically match
those of the full-gradient variant, i.e., $\mathcal{O}(1/t)$ convergence rate
for convex problems (Theorem 1) and a linear convergence rate under strong
convexity of $f$ and smoothness of $h$ (Theorem 3).

  * •

For the sparse variant, Sparse Vr-Tos (Algorithm 2), we obtain a linear
convergence rate under the same assumptions (Theorem 3). However, for general
convex objectives we could only obtain a worse $\mathcal{O}(1/\sqrt{t})$
convergence rate (Theorem 2).

In this section we will use the following extra notation. We define the
following primal ($\mathcal{P}$), and dual function ($\mathcal{D}$) as:

$\displaystyle\mathcal{P}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}f({\boldsymbol{x}})+g({\boldsymbol{x}})+h({\boldsymbol{x}})~{},$
$\displaystyle\mathcal{D}({\boldsymbol{u}})\stackrel{{\scriptstyle\text{def}}}{{=}}(f+g)^{*}(-{\boldsymbol{u}})+h^{*}({\boldsymbol{u}})~{},$
(4)

where ∗ denotes the Fenchel conjugate. We denote by ${\boldsymbol{x}}^{\star}$
an arbitrary minimizer of the primal objective and define the “dual iterate”
${\boldsymbol{u}}_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}{\boldsymbol{D}}^{-1}({\boldsymbol{y}}_{t}-{\boldsymbol{z}}_{t})/\gamma$
(${\boldsymbol{D}}=\boldsymbol{I}$ for the dense variants). We also define the
following generalized three operator splitting operator:

$\displaystyle\boldsymbol{G}_{\gamma}({\boldsymbol{y}})\stackrel{{\scriptstyle\text{def}}}{{=}}{\boldsymbol{y}}-{\boldsymbol{z}}_{\boldsymbol{y}}+\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
g}(2{\boldsymbol{z}}_{\boldsymbol{y}}-{\boldsymbol{y}}-\gamma{\boldsymbol{D}}\nabla
f({\boldsymbol{z}}_{\boldsymbol{y}}))\,,$ $\displaystyle~{}\text{with
${\boldsymbol{z}}_{\boldsymbol{y}}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}})$}~{},$ (5)

and its set of fixed points, which we denote
$\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$. Another quantity that will
appear often in the analysis is
$H_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}{1}/{(2nL_{f})}\sum_{i=1}^{n}\|{\boldsymbol{\alpha}}_{i,0}-\psi_{i}({\boldsymbol{x}}^{\star})\|^{2}$.

Throughout this section we make the following two technical assumptions:

#### Assumption 1: Regularity.

We assume each $\psi_{i}$ is $L_{\psi}$-smooth, $\omega$ is
$L_{\omega}$-smooth, $g$ and $h$ are proper (i.e., have nonempty domain),
lower semicontinuous (i.e., its sublevel sets are closed) convex functions. We
recall that lower semicontinuity is a weak form of continuity that allows
extended-valued functions with domain over a closed set.

Assumption 2: Qualification conditions. We assume the relative interior of
$\operatorname*{dom}g$ and $\operatorname*{dom}h$ have a non-empty
intersection. This is a very weak and standard assumption, which allows to
rule out pathological cases such as disjoint domains and allows to relate the
primal and dual optimal objective (see e.g.(Bauschke and Combettes,, 2017,
Proposition 15.13) or (Bertsekas,, 2015, Proposition 5.3.8)), a property
sometimes referred to as strong or total duality.

Sublinear convergence. The following theorem shows a $\mathcal{O}(1/t)$
convergence rate for Vr-Tos on arbitrary convex objectives.

One of the issues when analyzing the convergence of the three operator
splitting is that the objective function might be $+\infty$, for example when
both proximal terms are an indicator function. Following Chambolle and Pock,
(2015); Pedregosa and Gidel, (2018), we will state the convergence rate for
general functions in terms of the _saddle point suboptimality_ , defined as

$\displaystyle\mathcal{L}(\widetilde{\boldsymbol{x}},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},\widetilde{\boldsymbol{u}})\,,~{}\text{
with }~{}$
$\displaystyle\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}})\stackrel{{\scriptstyle\text{def}}}{{=}}f({\boldsymbol{x}})+g({\boldsymbol{x}})+\langle{\boldsymbol{x}},{\boldsymbol{u}}\rangle-h^{*}({\boldsymbol{u}})~{},$
(6)

where $\mathcal{L}$ is the Lagrangian associated with $\mathcal{P}$ and
$\mathcal{D}$. As Davis and Yin, (2017), we will also state convergence rates
in terms of the objective suboptimality under a Lipschitz assumption on $h$ in
(8).

######  Theorem 1.

Let $\overline{{\boldsymbol{x}}}_{t}$ denote the averaged (also known as
ergodic) iterate, i.e.,
$\overline{{\boldsymbol{x}}}_{t}={(\sum_{k=0}^{t}{\boldsymbol{x}}_{k})/(t+1)}$
and
$\overline{{\boldsymbol{u}}}_{t}=(\sum_{k=0}^{t}{\boldsymbol{u}}_{k})/(t+1)$.
Then the Vr-Tos method (Algorithm 1) converges for any step size $\gamma\leq
1/(3L_{f})$, and for $\gamma=1/(3L_{f})$ we have the following bound for all
$({\boldsymbol{x}},{\boldsymbol{u}})\in\operatorname*{dom}g\times\operatorname*{dom}h^{*}$:

${\mathbb{E}}\left[\mathcal{L}(\overline{\boldsymbol{x}}_{t},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},\overline{\boldsymbol{u}}_{t})\right]\leq\frac{10n}{q(t+1)}C_{0},$
(7)

with ${\boldsymbol{y}}={\boldsymbol{x}}+\gamma{\boldsymbol{u}}$,
${\boldsymbol{y}}^{\star}\in\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$, and
$C_{0}=\left[\frac{3L_{f}q}{20n}\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}\|^{2}+\frac{3L_{f}q}{2n}\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}^{\star}\|^{2}+H_{0}\right]~{}$,
where we recall
$H_{0}={1}/{(2nL_{f})}\sum_{i=1}^{n}\|{\boldsymbol{\alpha}}_{i,0}-\psi_{i}({\boldsymbol{x}}^{\star})\|^{2}$.

Furthermore, if $h$ is $\beta_{h}$-Lipschitz we have the following rate in
terms of the primal objective:

$\mathcal{P}(\overline{{\boldsymbol{x}}}_{t})-\mathcal{P}({\boldsymbol{x}}^{\star})\leq\frac{10n}{q(t+1)}\widetilde{C}_{0}~{},$
(8)

with
$\widetilde{C}_{0}=\frac{6L_{f}q}{20n}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\frac{3L_{f}q}{2n}\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}^{\star}\|^{2}+\frac{q}{15nL_{f}}\beta_{h}^{2}+H_{0}$.

The previous theorem gives a $\mathcal{O}(1/t)$ convergence rate in terms of
the saddle point suboptimality for arbitrary convex functions and
$\mathcal{O}(1/t)$ rate in function suboptimality under a Lipschitz assumption
on $h$, matching the strongest bounds of Saga (Defazio et al.,, 2014).

For their sparse variants, however, we have only been able to prove a slower
$\mathcal{O}(1/\sqrt{t})$ rate on the operator residual, despite the fact that
in practice the algorithm exhibits a much faster empirical convergence (see
§5). Appendix B contains a characterization of the fixed points of this
operator that justifies why this is a meaningful suboptimality criterion for
(OPT). Although there is no direct correspondence between rates on the
gradient and on objective values, lower bounds are asymptotically equivalent
(Nesterov,, 2012).

######  Theorem 2.

Sparse Vr-Tos (Algorithm 2) converges for every step size ${\gamma\leq
1/(3L_{f})}$. In particular, for $\gamma={1}/{(3L_{f})}$ and
${\boldsymbol{y}}_{t}$ obtained after $t\geq 1$ updates we have the bound

$\min_{k=0,\ldots,t}\left\\{{\mathbb{E}}\|{\boldsymbol{y}}_{k}-\boldsymbol{G}_{\gamma}({\boldsymbol{y}}_{k})\|\right\\}\leq\sqrt{\frac{C_{0}}{Lq(t+1)}}=\mathcal{O}\left(\frac{1}{\sqrt{t}}\right),$
(9)

with
$C_{0}=\frac{5d_{\max}n}{{Lq(t+1)}}\left[({2Lq}/{n})\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}^{\star}\|^{2}+H_{0}\right]$.

Linear convergence. The three operator splitting has been shown to have a
linear convergence rate under the assumption of strong convexity of the smooth
term and smoothness of one of the proximal terms (Davis and Yin,, 2015, §4.4).
Although this last condition is rarely verified in practice since its main
application is on non-smooth proximal terms, it is instructive to see that the
proposed method –despite the reduced cost per iteration– also enjoys a linear
convergence rate under the same assumptions.

######  Theorem 3 (Linear convergence).

Let $\psi_{i}$ be $\mu_{\psi}$-strongly convex and $\omega$ be
$\mu_{\omega}$-strongly convex, where $\mu_{\psi}+\mu_{\omega}>0$.
Furthermore, let $h$ be $L_{h}$-smooth. Then for any step size
$\gamma\leq{1}/{(3L_{f})}$, all the proposed methods converge geometrically in
expectation. For $\gamma={1}/{(3L_{f})}$, we have the following bound for
Algorithm 1 ($d_{\max{}}=1$ in this case) and Algorithm 2:

${\mathbb{E}}\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}\leq\left(1-\min\Big{\\{}\frac{q}{4n},\frac{1}{3d^{3}_{\max}\delta^{2}\kappa}\Big{\\}}\right)^{t}D_{0}\quad,$
(10)

with
$D_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}{d_{\max}}\left[\frac{q}{2\gamma(1-\gamma\mu)n}\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}^{\star}\|^{2}+H_{0}\right]$,
$\delta=(1+L_{h}/(3L_{f}))$, $\kappa=L_{f}/\mu$ and
${\boldsymbol{y}}^{\star}\in\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$.

| Method | step size | Proximal oracle | Convergence rate | Extra assumptions
---|---|---|---|---|---
Geometric | Saga |  |  |  |
(Defazio et al.,, 2014) | $\nicefrac{{$1$\,}}{{\,$3L_{f}$}}$ | $\mathbf{prox}_{\gamma(g+h)}$ | $\displaystyle\Big{(}1-\min\big{\\{}\textstyle\frac{1}{4n},\frac{1}{3\kappa}\big{\\}}\Big{)}^{t}C_{0}$ | Each $\psi_{i}$ is $\mu$-cvx
ProxSvrg |  |  |  | $f$ is $\mu$-cvx
(Xiao and Zhang,, 2014) | $\nicefrac{{$1$\,}}{{\,$10L_{f}$}}$ | $\mathbf{prox}_{\gamma(g+h)}$ | $\Big{(}\frac{1}{\kappa 0.6m}+\frac{2}{3}\Big{)}^{t}C_{0}$
Vr-Tos |  |  |  | Each $\psi_{i}$ is $\mu$-cvx
(this work) | $\nicefrac{{$1$\,}}{{\,$3L_{f}$}}$ | $\mathbf{prox}_{\gamma g}~{},~{}\mathbf{prox}_{\gamma h}$ | $\displaystyle\Big{(}1-\min\big{\\{}\textstyle\frac{q}{4n},\frac{1}{3d_{\max{}}^{2}\delta^{2}\kappa}\big{\\}}\Big{)}^{t}C_{0}$ | and $h$ is $L_{h}$-smooth
|  |  |  |
Saga |  |  |  |
(Defazio et al.,, 2014) | $\nicefrac{{$1$\,}}{{\,$3L_{f}$}}$ | $\mathbf{prox}_{\gamma(g+h)}$ |  $\mathcal{O}(1/t)$ | None
Sublinear | Stochastic Tos |  |  |  | $f$ is $\mu$-cvx +
| (Yurtsever et al.,, 2016) | $\mathcal{O}$($\nicefrac{{1\,}}{{\,$t$}}$) | $\mathbf{prox}_{\gamma g}~{},~{}\mathbf{prox}_{\gamma h}$ | $\mathcal{O}(1/t)$ | bound on gradients
Vr-Tos |  |  |  |
(this work, dense/sparse variant) | $\nicefrac{{$1$\,}}{{\,$3L_{f}$}}$ | $\mathbf{prox}_{\gamma g}~{},~{}\mathbf{prox}_{\gamma h}$ | $\mathcal{O}(1/{t})$ / $\mathcal{O}(1/\sqrt{t})$ | None
Table 2: Assumptions and properties of related incremental methods. In every
case, we take the step size recommended by the theory, where we assume
$\omega=0$ to make them comparable. Proximal oracle is the proximal operators
that are needed by the algorithm. Extra assumptions refer to those other than
Assumptions 1 and 2. The linear rates use the quantities $\delta=(1+\gamma
L_{h})$, $\kappa=L_{f}/\mu$. For ProxSvrg, $m$ denotes the epoch size and the
convergence rate is relative to the number of epochs and not iterations like
the rest.

###  4.1 Discussion

Comparison of convergence rates. We summarize the obtained convergence rates
for the proposed methods and compare them against the best known rates for
related stochastic methods in Table 2. In the linearly-convergent regime, we
obtain rates that are similar to Saga but with the rate factor multiplied by
${1}/{(\delta^{2}d^{3}_{\max})}$, quantity that depends on the smoothness of
$g$ and the sparsity of the gradients.

An improved ProxSvrg variant. The analysis of ProxSvrg (Xiao and Zhang,, 2014)
requires that the step size verifies an implicit equation that depends among
other things on the strong convexity parameter. For typical choices of the
parameters this is $1/(10L_{f})$ (Xiao and Zhang,, 2014, Theorem 1). In
contrast, Sparse Vr-Tos with Svrg-like sampling with $h=0$ yields a variant of
ProxSvrg with more favorable properties. First, none of its parameters depend
on the strong convexity constant (while still obtaining a linear convergence
rate since $L_{h}=0$ in this case), which is most often unknown. Second, it
admits the much larger step size $1/(3L_{f})$, which is, to the best of our
knowledge, the largest step size of any Svrg variant. Third, it can leverage
sparsity in the input data through sparse updates.

Linear convergence without smoothness of the proximal term. Theorem 3 requires
smoothness of one of the proximal terms to guarantee linear convergence.
Despite this, linear convergence is observed in practice without this
assumption (Figure 1). This has also been observed in the case of the original
(non-variance reduced) three operator splitting (Davis and Yin,, 2017;
Pedregosa and Gidel,, 2018), although an explanation for this is still an open
problem. Furthermore, the lack of linear convergence when both proximal terms
are non-smooth does not seem to be a limitation of the proof, as a
counterexample was provided in (Davis and Yin,, 2015, Appendix D.6). In this
work, the authors constructed a strongly monotone operator with a sublinear
convergence.

Step size adaptivity to linear convergence. A practical consequence of the
above theorems is that using the same step size $\gamma={1}/{(3L_{f})}$ we
obtain a sublinear convergence by Theorem 1 and a linear rate (under
additional assumptions) by Theorem 3. That is, one can use the “universal”
step size ${1}/{(3L_{f})}$ and automatically obtain linear convergence
whenever the assumptions of Theorem 3 are verified.

Limitations. The following are some scenarios under which the proposed method
is expected to perform poorly. The cost in computation and storage scales
linearly with the number of proximal terms, hence it cannot cope with other
scenarios with many nonsmooth terms such as empirical risk minimization with
the hinge loss or group lasso with overlap with a large number of overlaps
(for instance $>100$). Also, there are still penalties that cannot be reduced
to a sum of proximal terms, such as the nuclear norm. Algorithms based on
Frank-Wolfe (Jaggi,, 2013) or with approximate proximal operators (Schmidt et
al.,, 2011) might be better suited in such regimes.

##  5 Experiments

Dataset | #samples | #dimensions | density | $L_{f}/\mu$
---|---|---|---|---
RCV1 (full) (Lewis et al.,, 2004) | 697,641 | 47,236 | $1.5\times 10^{-3}$ | 2.50 $\times 10^{4}$
URL (Ma et al.,, 2009) | 2,396,130 | 3,231,961 | $3.5\times 10^{-5}$ | 1.28 $\times 10^{7}$
KDD10 (Yu et al.,, 2010) | 19,264,097 | 29,890,095 | 9.8 $\times 10^{-7}$ | $5.2\times 10^{8}$
Criteo (Juan et al.,, 2016) | 45,840,617 | 1,000,000 | $3.8\times 10^{-5}$ | $1.1\times 10^{7}$

Figure 1: Top: Description of considered datasets. Bottom: Suboptimality vs
time of different algorithms on a logistic regression with overlapping group
lasso penalty problem.

Although the proposed methods can be applied more broadly, we consider for the
experiments a logistic regression problem with squared $\ell_{2}$
regulrization and an overlapping group lasso penalty (Jacob et al.,, 2009).
Following Jacob et al., (2009) we choose groups of 10 variables with 2
variables of overlap between two successive groups:
$\\{\\{1,\ldots,10\\},\\{8,\ldots,18\\},\\{16,\ldots,26\\},\ldots\\}$. The
amount of group regularization was chosen such that the solution has roughly
$10\%$ of non-zero coefficients and the of $\ell_{2}$ regularization was fixed
to ${1}/{n}$. We consider the following methods:

  * •

The proposed method Sparse Vr-Tos (Algorithm 2), where the overlapping group
lasso penalty is split as a sum of two non-overlapping group lasso penalties,
for which the proximal operator is available in closed form. We used the
formulation with 3 proximal terms of §2.2 to better leverage sparsity in the
dataset and consider Saga and Svrg-like updates, denoted Vr-Tos (Saga variant)
and Vr-Tos (Svrg variant) respectively. This implementation is publicly
available in the C-OPT package.444http://openopt.github.io/copt/

It is worth noting that while original penalty is _not_ block separable, each
of the terms in the splitting as two group lasso penalties _is_ block
separable. This will allow us to make a much more efficient use of sparsity
than what is possible on on methods like Saga and ProxSvrg.

  * •

The three operator splitting (denoted Tos), in its recently proposed variant
with adaptive step size (Pedregosa and Gidel,, 2018).

  * •

The stochastic three operator spitting of (Yurtsever et al.,, 2016) with the
same splitting as Vr-Tos, denoted Stos.

  * •

Saga and ProxSvrg, where the proximal operator is evaluated approximately
using 10 iterations of the Douglas-Rachford method.

The above methods were compared on 4 large-scale datasets described in the
table of Figure 1. Further details and implementation aspects are discussed in
Appendix F.1.

The best performing algorithms overall are the proposed Vr-Tos variants, which
are over an order of magnitude faster than the second best method, the
adaptive three operator splitting. The stochastic three operator splitting,
not being able to take advantage of the sparsity in the gradients, performs
poorly in this benchmark, appearing as a straight line. Saga and ProxSvrg were
the slowest since they require to compute a costly proximal operator at each
iteration and are unable to leverage the sparsity of the dataset due to the
non-block-separability of the non-smooth term.

It is worth noting from Figure 1 that the two variants of Sparse Vr-Tos
exhibit an empirical linear convergence, despite the fact that the theory only
predicts in this regime a much slower $\mathcal{O}(1/\sqrt{t})$ convergence
rate (Theorem 1).

We provide extra experiments in Appendix F.2.

##  6 Future work

This work can be extended in several ways. As highlighted in §4.1, a
theoretical explanation for the empirical linear convergence without smoothess
of any proximal term, even for the full gradient algorithm, is lacking. We
conjecture _partly smooth_ is a sufficient condition on the penalties to
ensure local linear convergence, as recently proven for related methods (Liang
et al.,, 2018). Second, we conjecture that the convergence rate of the sparse
variant can be improved up to to $\mathcal{O}(1/t)$. A third direction for
future work would be the development an extension that allow for a linear
operator inside one of the proximal terms, as in (Condat, 2013b, ; Zhao and
Cevher,, 2018; Yan,, 2018).

## Acknowledgements

The authors warmly thank Vincent Roulet, Vlad Niculae, Rémi Leblond and
Federico Vaggi for their feedback on the manuscript, as well as Adrien Taylor,
Alexandre D’Aspremont, Gabriel Peyré, Guillaume Obozinski, P. Balamurugan,
Francis Bach and Marwa El Halabi for fruitful discussions.

This work has been done while FP was under funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodorowska-
Curie grant agreement 748900. KF is funded through the project OATMIL
ANR-17-CE23-0012 of the French National Research Agency (ANR). Computing time
on was donated by Amazon through the program “AWS Cloud Credits for Research”.

## References

  * Balamurugan and Bach, (2016) Balamurugan, P. and Bach, F. (2016).  Stochastic Variance Reduction Methods for Saddle-Point Problems.  Advances in Neural Information Processing Systems.
  * Barbero and Sra, (2014) Barbero, Á. and Sra, S. (2014).  Modular proximal optimization for multidimensional total-variation regularization.  arXiv preprint arXiv:1411.0589.
  * Bauschke et al., (2012) Bauschke, H. H., Boţ, R. I., Hare, W. L., and Moursi, W. M. (2012).  Attouch–Théra duality revisited: paramonotonicity and operator splitting.  Journal of Approximation Theory.
  * Bauschke and Combettes, (2017) Bauschke, H. H. and Combettes, P. L. (2017).  Convex analysis and monotone operator theory in Hilbert spaces.  Springer Science & Business Media.
  * Bertsekas, (2015) Bertsekas, D. P. (2015).  Convex optimization algorithms.  Athena Scientific Belmont.
  * Chambolle and Pock, (2015) Chambolle, A. and Pock, T. (2015).  On the ergodic convergence rates of a first-order primal–dual algorithm.  Mathematical Programming.
  * (7) Condat, L. (2013a).  A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms.  Journal of Optimization Theory and Applications.
  * (8) Condat, L. (2013b).  A direct algorithm for 1D total variation denoising.  IEEE Signal Processing Letters.
  * Davis and Yin, (2015) Davis, D. and Yin, W. (2015).  A three-operator splitting scheme and its optimization applications.  preprint arXiv:1504.01032v1.
  * Davis and Yin, (2017) Davis, D. and Yin, W. (2017).  A three-operator splitting scheme and its optimization applications.  Set-Valued and Variational Analysis.
  * Defazio et al., (2014) Defazio, A., Bach, F., and Lacoste-Julien, S. (2014).  SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives.  In Advances in Neural Information Processing Systems.
  * Giselsson and Boyd, (2016) Giselsson, P. and Boyd, S. (2016).  Linear Convergence and Metric Selection in Douglas-Rachford Splitting and ADMM.  IEEE Transactions on Automatic Control.
  * Hofmann et al., (2015) Hofmann, T., Lucchi, A., Lacoste-Julien, S., and McWilliams, B. (2015).  Variance Reduced Stochastic Gradient Descent with Neighbors.  In Advances in Neural Information Processing Systems.
  * Iusem, (1998) Iusem, A. N. (1998).  On Some Properties of Generalized Proximal Point Methods for Variational Inequalities.  Journal of Optimization Theory and Applications.
  * Jacob et al., (2009) Jacob, L., Obozinski, G., and Vert, J.-P. (2009).  Group lasso with overlap and graph lasso.  In Proceedings of the 26th annual international conference on machine learning. ACM.
  * Jaggi, (2013) Jaggi, M. (2013).  Revisiting Frank-Wolfe: projection-free sparse convex optimization.  In International Conference on Machine Learning.
  * Johnson, (2013) Johnson, N. (2013).  A dynamic programming algorithm for the fused lasso and $L_{0}$-segmentation.  Journal of Computational and Graphical Statistics.
  * Johnson and Zhang, (2013) Johnson, R. and Zhang, T. (2013).  Accelerating stochastic gradient descent using predictive variance reduction.  In Advances in Neural Information Processing Systems.
  * Juan et al., (2016) Juan, Y., Zhuang, Y., Chin, W.-S., and Lin, C.-J. (2016).  Field-aware factorization machines for CTR prediction.  In Proceedings of the 10th ACM Conference on Recommender Systems. ACM.
  * Kim et al., (2009) Kim, S.-J., Koh, K., Boyd, S., et al. (2009).  $\ell_{1}$ trend filtering.  SIAM review.
  * Le Roux et al., (2012) Le Roux, N., Schmidt, M., and Bach, F. (2012).  A stochastic gradient method with an exponential convergence rate for finite training sets.
  * Lewis et al., (2004) Lewis, D. D., Yang, Y., Rose, T. G., and Li, F. (2004).  RCV1: A new benchmark collection for text categorization research.  Journal of machine learning research, 5(Apr):361–397.
  * Liang et al., (2018) Liang, J., Fadili, J., and Peyré, G. (2018).  Local linear convergence analysis of primal–dual splitting methods.  Optimization.
  * Ma et al., (2009) Ma, J., Saul, L. K., et al. (2009).  Identifying suspicious URLs: an application of large-scale online learning.  In Proceedings 26th ACM international conference on machine learning.
  * Nesterov, (2004) Nesterov, Y. (2004).  Introductory lectures on convex optimization.  Springer.
  * Nesterov, (2012) Nesterov, Y. (2012).  How to make the gradients small.  Optima.
  * Pedregosa and Gidel, (2018) Pedregosa, F. and Gidel, G. (2018).  Adaptive Three Operator Splitting.  Proceedings of the 35th International Conference on Machine Learning.
  * Pedregosa et al., (2017) Pedregosa, F., Leblond, R., and Lacoste-Julien, S. (2017).  Breaking the Nonsmooth Barrier: A Scalable Parallel Method for Composite Optimization.  Advances in Neural Information Processing System 30 (NIPS).
  * Raguet et al., (2013) Raguet, H., Fadili, J., and Peyré, G. (2013).  A generalized forward-backward splitting.  SIAM Journal on Imaging Sciences.
  * Robbins and Monro, (1951) Robbins, H. and Monro, S. (1951).  A Stochastic Approximation Method.  Ann. Math. Statist.
  * Rockafellar, (1997) Rockafellar, R. T. (1997).  Convex analysis.
  * Rockafellar and Wets, (1998) Rockafellar, R. T. and Wets, R. J.-B. (1998).  Variational analysis.  Springer.
  * Schmidt et al., (2011) Schmidt, M., Le Roux, N., and Bach, F. (2011).  Convergence rates of inexact proximal-gradient methods for convex optimization.  In Advances in neural information processing systems 24.
  * Shalev-Shwartz and Zhang, (2013) Shalev-Shwartz, S. and Zhang, T. (2013).  Stochastic dual coordinate ascent methods for regularized loss minimization.  Journal of Machine Learning Research.
  * Vũ, (2013) Vũ, B. C. (2013).  A splitting algorithm for dual monotone inclusions involving cocoercive operators.  Advances in Computational Mathematics.
  * Xiao and Zhang, (2014) Xiao, L. and Zhang, T. (2014).  A proximal stochastic gradient method with progressive variance reduction.  SIAM Journal on Optimization.
  * Yan, (2018) Yan, M. (2018).  A New Primal–Dual Algorithm for Minimizing the Sum of Three Functions with a Linear Operator.  Journal of Scientific Computing.
  * Yu et al., (2010) Yu, H.-F., Lo, H.-Y., Hsieh, H.-P., Lou, J.-K., McKenzie, T. G., Chou, J.-W., Chung, P.-H., Ho, C.-H., Chang, C.-F., Wei, Y.-H., et al. (2010).  Feature engineering and classifier ensemble for KDD cup 2010.  In KDD Cup.
  * Yu and Huang, (2017) Yu, Y. and Huang, L. (2017).  Fast stochastic variance reduced admm for stochastic composition optimization.  In Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence (IJCAI-17).
  * Yurtsever et al., (2016) Yurtsever, A., Vu, C. B., and Cevher, V. (2016).  Stochastic Three-Composite Convex Minimization.  In Advances in Neural Information Processing Systems.
  * Zhao and Cevher, (2018) Zhao, R. and Cevher, V. (2018).  Stochastic Three-Composite Convex Minimization with a Linear Operator.  In Proceedings of the 21st International Conference on Artificial Intelligence and Statistics.
  * Zheng and Kwok, (2016) Zheng, S. and Kwok, J. T. (2016).  Stochastic variance-reduced ADMM.  arXiv preprint arXiv:1604.07070.

Proximal Splitting Meets Variance Reduction

Supplementary material

#### Outline.

The supplementary material of this paper is organized as follows.

  * •

Appendix A presents basic definitions and properties that will be used
throughout the proofs but which are not specific to our methods. Most of these
can be found in convex optimization textbooks, such as (Bauschke and
Combettes,, 2017; Nesterov,, 2004).

  * •

Appendix B give a characterization of the fixed points of the three operator
splitting, relating the set of fixed points of the three operator splitting to
the solutions of primal and dual objectives. This is a stronger result than
the one stated in (Davis and Yin,, 2017) and used in some of our proofs.

  * •

Appendix C gives the proofs of those results in the Analysis section of the
paper.

  * •

Appendix D discusses splitting strategies for different penalties and examines
some cases in which the scaled proximal operator can be computed in closed
form.

  * •

Appendix F discusses implementation aspects of the proposed algorithms.

##  Appendix Appendix A Basic definitions and properties

######  Definition 2 (proper function).

A function $f:\mathcal{X}\subseteq{\mathbb{R}}^{p}\to]-\infty,\infty]$ is said
to be proper if its domain is not empty.

######  Definition 3 (Fenchel conjugate).

The Fenchel conjugate of a function
$f:\mathcal{X}\subseteq{\mathbb{R}}^{p}\to]-\infty,\infty]$ is defined as

$f^{*}({\boldsymbol{x}}^{\star})=\sup\left\\{\left.\left\langle{\boldsymbol{x}}^{\star},{\boldsymbol{x}}\right\rangle-f\left({\boldsymbol{x}}\right)\right|{\boldsymbol{x}}\in\mathcal{X}\right\\}.$
(11)

######  Definition 4 (lower semicontinuity).

We say that a proper convex function $f$ is lower-semicontinuous if all of its
levelsets
${\\{{\boldsymbol{x}}\in\operatorname*{dom}(f)\,|\,f({\boldsymbol{x}})\leq\alpha\\}}$
are closed.

######  Definition 5 (relative interior).

The relative interior of a convex set $C\subseteq{\mathbb{R}}^{p}$ is defined
as

$\operatorname{relint}(C)\stackrel{{\scriptstyle\text{def}}}{{=}}\\{{\boldsymbol{x}}\in
C:\forall{{\boldsymbol{y}}\in
C}\;\exists{\lambda>1}:\lambda{\boldsymbol{x}}+(1-\lambda){\boldsymbol{y}}\in
C\\}$ (12)

######  Definition 6 (Bregman divergence).

The Bregman divergence associated with a convex function $f$ for points
${\boldsymbol{x}},{\boldsymbol{y}}$ in its domain is defined as:

$B_{f}({\boldsymbol{x}},{\boldsymbol{y}})\stackrel{{\scriptstyle\text{def}}}{{=}}f({\boldsymbol{x}})-f({\boldsymbol{y}})-\langle\nabla
f({\boldsymbol{y}}),{\boldsymbol{x}}-{\boldsymbol{y}}\rangle$

Note that this is always positive due to the convexity of $f$.

######  Definition 7 (Proximal operators).

Here, we redefine 2 variants of a critical notion. The proximal operator is
defined for a function $\varphi$, step size $\gamma>0$ as:

${\operatorname{\mathbf{prox}}_{\gamma\varphi}}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}{\operatorname*{{arg\,min}}_{{\boldsymbol{z}}\in{\mathbb{R}}^{p}}}\big{\\{}\,\varphi({\boldsymbol{z}})+\frac{1}{2\gamma}\|{\boldsymbol{x}}-{\boldsymbol{z}}\|\big{\\}}\,$
(13)

The _scaled proximal opeartor_ is defined for a function $\varphi$, step size
$\gamma>0$ and positive definite matrix $\boldsymbol{H}$ as the solution of
the following optimization problem

${\operatorname{\mathbf{prox}}^{\boldsymbol{H}}_{\gamma\varphi}}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}{\operatorname*{{arg\,min}}_{{\boldsymbol{z}}\in{\mathbb{R}}^{p}}}\big{\\{}\,\varphi({\boldsymbol{z}})+\frac{1}{2\gamma}\|{\boldsymbol{x}}-{\boldsymbol{z}}\|_{\boldsymbol{H}}^{2}\,\big{\\}}\,~{}\text{
with
$\|\cdot\|_{{\boldsymbol{H}}}^{2}\stackrel{{\scriptstyle\text{def}}}{{=}}\langle\cdot,{\boldsymbol{H}}\cdot\rangle$}.$
(14)

######  Lemma 1 (subgradient characterization of proximal operator).

Let $g$ be a convex proper lower semicontinuous function. Then for any
${\boldsymbol{x}}$, ${\boldsymbol{H}}$ positive definite matrix and any
$\gamma>0$ we have the following characterization of proximal operator:

${\boldsymbol{z}}=\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}}_{\gamma
g}({\boldsymbol{x}})\iff{\boldsymbol{H}}({\boldsymbol{x}}-{\boldsymbol{z}})/{\gamma}\in\partial
g({\boldsymbol{z}})$ (15)

###### Proof.

By the definition of proximal operator we have that
${\boldsymbol{z}}=\operatorname{\mathbf{prox}}_{\gamma
g}^{{\boldsymbol{H}}}({\boldsymbol{x}})$ is equivalent to

$\displaystyle{\boldsymbol{z}}\in\operatorname*{{arg\,min}}_{{\boldsymbol{z}}^{\prime}\in{\mathbb{R}}^{p}}g({\boldsymbol{z}}^{\prime})+\frac{1}{2\gamma}\|{\boldsymbol{z}}^{\prime}-{\boldsymbol{x}}\|_{{\boldsymbol{H}}}^{2}$
(16) $\displaystyle\iff 0\in\partial
g({\boldsymbol{z}})+\frac{{\boldsymbol{H}}}{\gamma}({\boldsymbol{z}}-{\boldsymbol{x}})$
(17)
$\displaystyle\iff\frac{{\boldsymbol{H}}}{\gamma}({\boldsymbol{x}}-{\boldsymbol{z}})\in\partial
g({\boldsymbol{z}})$ (18)

where the first equivalence is a consequence of the first order optimality
conditions. ∎

######  Lemma 2 (Conjugate-inverse identity).

Let $h$ be a convex, proper lower semicontinuous function. Then

${\boldsymbol{u}}\in\partial
h({\boldsymbol{z}})\iff{\boldsymbol{z}}\in\partial
h^{*}({\boldsymbol{u}})~{}.$ (19)

In other words, $(\partial h)^{-1}=\partial h^{*}$.

###### Proof.

See e.g. (Bauschke and Combettes,, 2017, Corollary 16.30) or (Rockafellar and
Wets,, 1998, Proposition 11.3). ∎

######  Lemma 3 (Generalized variance decomposition).

Let $\zeta_{i}$ be a random variable and let $\mathbf{E}$ the expectation with
respect to this random variable. Furthermore, let ${\boldsymbol{Q}}_{i}$ be an
orthogonal projection such that
${\boldsymbol{Q}}_{i}\boldsymbol{\zeta}_{i}=\boldsymbol{\zeta}_{i}$,
$\mathbf{E}{\boldsymbol{Q}}_{i}$ is invertible and
$\boldsymbol{A}=(\mathbf{E}{\boldsymbol{Q}}_{i})^{-1}$. Then we have

$\mathbf{E}\|\boldsymbol{\zeta}_{i}-{\boldsymbol{Q}}_{i}\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\|^{2}=\mathbf{E}\|\boldsymbol{\zeta}_{i}\|^{2}-\|\mathbf{E}\boldsymbol{\zeta}_{i}\|^{2}_{\boldsymbol{A}}~{}.$
(20)

###### Proof.

The assumption of ${\boldsymbol{Q}}_{i}$ being an orthogonal projection
implies that it is symmetric and idempotent. Developing the square we have

$\displaystyle\mathbf{E}\|\boldsymbol{\zeta}_{i}-{\boldsymbol{Q}}_{i}\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\|^{2}$
$\displaystyle=\mathbf{E}\|\boldsymbol{\zeta}_{i}\|^{2}+\mathbf{E}\langle{\boldsymbol{Q}}_{i}\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i},{\boldsymbol{Q}}_{i}\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\rangle-2\mathbf{E}\langle\boldsymbol{\zeta}_{i},{\boldsymbol{Q}}_{i}\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\rangle$
(21)
$\displaystyle=\mathbf{E}\|\boldsymbol{\zeta}_{i}\|^{2}+\mathbf{E}\langle{\boldsymbol{Q}}_{i}{\boldsymbol{Q}}_{i}\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i},\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\rangle-2\mathbf{E}\langle\boldsymbol{{\boldsymbol{Q}}}_{i}\zeta_{i},\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\rangle$
(22) (by symmetry of ${\boldsymbol{Q}}_{i}$)
$\displaystyle=\mathbf{E}\|\boldsymbol{\zeta}_{i}\|^{2}+\mathbf{E}\langle{\boldsymbol{Q}}_{i}\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i},\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\rangle-2\mathbf{E}\langle\boldsymbol{\zeta}_{i},\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\rangle$
(23) (idempotence of ${\boldsymbol{Q}}_{i}$ and assumption
${\boldsymbol{Q}}_{i}\boldsymbol{\zeta}_{i}=\boldsymbol{\zeta}_{i}$
respectively)
$\displaystyle=\mathbf{E}\|\boldsymbol{\zeta}_{i}\|^{2}+\langle\mathbf{E}\boldsymbol{\zeta}_{i},\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\rangle-2\langle\mathbf{E}\boldsymbol{\zeta}_{i},\boldsymbol{A}\mathbf{E}\boldsymbol{\zeta}_{i}\rangle$
(24) (taking expectations )
$\displaystyle=\mathbf{E}\|\boldsymbol{\zeta}_{i}\|^{2}-\|\mathbf{E}\boldsymbol{\zeta}_{i}\|^{2}_{\boldsymbol{A}}$
(25)

∎

######  Lemma 4 (Smooth inequality 1).

Let $f_{i}$ be $L_{f}$-smooth and convex for $i=1,\ldots,n$. Then it is
verified that

$\frac{1}{n}\sum_{i=1}^{n}\|\nabla f_{i}({\boldsymbol{x}})-\nabla
f_{i}({\boldsymbol{y}})\|^{2}\leq
2L_{f}B_{f}({\boldsymbol{x}},{\boldsymbol{y}})~{}.$ (26)

###### Proof.

Since each $f_{i}$ is $L_{f}$-smooth, it is verified (Nesterov,, 2004, Theorem
2.1.5) that

$\|\nabla f_{i}({\boldsymbol{x}})-\nabla f_{i}({\boldsymbol{y}})\|^{2}\leq
2L_{f}(f_{i}({\boldsymbol{x}})-f_{i}({\boldsymbol{y}})-\langle\nabla
f_{i}({\boldsymbol{y}}),{\boldsymbol{x}}-{\boldsymbol{y}}\rangle)$ (27)

The result is obtained by averaging over $i$. ∎

######  Lemma 5 (Bound on matrix norm).

Let ${\boldsymbol{D}}$ be a diagonal matrix with strictly positive diagonal
elements, let $d_{\max}$ and $d_{\min}$ denote its maximum and minimum
diagonal entry respectively. Then for any
${\boldsymbol{x}}\in{\mathbb{R}}^{p}$ we have the following inequalities

$d_{\min{}}\|{\boldsymbol{x}}\|^{2}\leq\|{\boldsymbol{x}}\|^{2}_{{\boldsymbol{D}}}\leq
d_{\max}\|{\boldsymbol{x}}\|^{2}$ (28)

###### Proof.

By definition of the ${\boldsymbol{D}}$-norm we have

$\displaystyle\|{\boldsymbol{x}}\|^{2}_{{\boldsymbol{D}}}=\sum_{i=1}^{n}{\boldsymbol{D}}_{i,i}{{\boldsymbol{x}}_{i}^{2}}\leq\sum_{i=1}^{n}d_{\text{max}}{{\boldsymbol{x}}_{i}^{2}}=d_{\text{max}}\|{\boldsymbol{x}}\|^{2}$
(29)
$\displaystyle\|{\boldsymbol{x}}\|^{2}_{{\boldsymbol{D}}}=\sum_{i=1}^{n}{\boldsymbol{D}}_{i,i}{{\boldsymbol{x}}_{i}^{2}}\geq\sum_{i=1}^{n}d_{\text{min}}{{\boldsymbol{x}}_{i}^{2}}=d_{\text{min}}\|{\boldsymbol{x}}\|^{2}$
(30)

The result follows from chaining both inequalities ∎

######  Lemma 6 (Properties of proximal operator).

Let $g$ be a convex lower semicontinuous function and ${\boldsymbol{H}}$ a
symmetric positive definite matrix. Then for all
${\boldsymbol{y}},\widetilde{{\boldsymbol{y}}}$ we have the following
inequality, often referred to as firm nonexpansiveness:

$\|\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}({\boldsymbol{y}})-\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}(\widetilde{{\boldsymbol{y}}})\|_{{\boldsymbol{H}}}^{2}\leq\langle\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}({\boldsymbol{y}})-\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}(\widetilde{{\boldsymbol{y}}}),{\boldsymbol{y}}-\widetilde{{\boldsymbol{y}}}\rangle_{{\boldsymbol{H}}}~{}$
(31)

Furthermore, if $g$ is $L_{g}$-smooth and ${\boldsymbol{H}}$ has smallest
singular value $\sigma_{\min}$ and largest singular value $\sigma_{\max}$,
then we also have the following bound:

$\|{\boldsymbol{y}}-\widetilde{{\boldsymbol{y}}}\|\leq(\sigma_{\max}/\sigma_{\min}+\gamma
L\sigma_{\max})\|\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}({\boldsymbol{y}})-\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}(\widetilde{{\boldsymbol{y}}})\|$ (32)

###### Proof.

_First inequality_. Let
${\boldsymbol{z}}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}({\boldsymbol{y}})$,
$\widetilde{\boldsymbol{z}}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}(\widetilde{{\boldsymbol{y}}})$. By the subgradient characterization of
Lemma 1 we have

${\boldsymbol{H}}^{-1}({\boldsymbol{y}}-{\boldsymbol{z}})/\gamma\in\partial
g({\boldsymbol{z}})\quad,\quad{\boldsymbol{H}}^{-1}(\widetilde{\boldsymbol{y}}-\widetilde{\boldsymbol{z}})/\gamma\in\partial
g(\widetilde{\boldsymbol{z}})$ (33)

Since the subdifferential of a convex function is monotonous, in particular
$\partial g$ is monotonous, and so we have

$\displaystyle\langle{\boldsymbol{H}}^{-1}({\boldsymbol{y}}-{\boldsymbol{z}})/\gamma-{\boldsymbol{H}}^{-1}(\widetilde{\boldsymbol{y}}-\widetilde{\boldsymbol{z}})/\gamma,{\boldsymbol{z}}-\widetilde{\boldsymbol{z}}\rangle\geq
0$
$\displaystyle\iff\langle({\boldsymbol{y}}-{\boldsymbol{z}})-(\widetilde{\boldsymbol{y}}-\widetilde{\boldsymbol{z}}),{\boldsymbol{z}}-\widetilde{\boldsymbol{z}}\rangle_{{\boldsymbol{H}}^{-1}}\\!\geq
0$ (34)
$\displaystyle\iff\langle{\boldsymbol{y}}-\widetilde{\boldsymbol{y}},{\boldsymbol{z}}-\widetilde{\boldsymbol{z}}\rangle_{{\boldsymbol{H}}^{-1}}\geq\|{\boldsymbol{z}}-\widetilde{\boldsymbol{z}}\|_{{\boldsymbol{H}}^{-1}}$
(35)

which proves the first part of the lemma (firm nonexpansive).

_Second inequality_. To prove the second inequality we will use a
generalization of the argument from Giselsson and Boyd, (2016, Proposition 1).
Let $g_{\gamma}$ be defined as

$g_{\gamma}\stackrel{{\scriptstyle\text{def}}}{{=}}\gamma
g+\frac{1}{2}\|\cdot\|_{{\boldsymbol{H}}^{-1}}^{2}\quad.$ (36)

By the subgradient characterization of the proximal operator (Lemma 1) and the
conjugate-inverse identity (Lemma 2), we have

$\displaystyle{\boldsymbol{z}}=\operatorname{\mathbf{prox}}^{{\boldsymbol{H}}^{-1}}_{\gamma
g}({\boldsymbol{y}})$ (37)
$\displaystyle\Longleftrightarrow{\boldsymbol{z}}\in\\{{\boldsymbol{x}}\,|\,{\boldsymbol{H}}^{-1}({\boldsymbol{y}}-{\boldsymbol{x}})\in\gamma\nabla
g({\boldsymbol{x}})\\}\qquad\text{(Lemma~{}\ref{lemma:prox_characterization})}$
(38)
$\displaystyle{\Longleftrightarrow}{\boldsymbol{z}}\in\\{{\boldsymbol{x}}\,|\,{\boldsymbol{H}}^{-1}{\boldsymbol{y}}\in\nabla
g_{\gamma}({\boldsymbol{x}})\\}$ (39)
$\displaystyle\Longleftrightarrow{\boldsymbol{z}}\in(\nabla
g_{\gamma})^{-1}({\boldsymbol{H}}^{-1}{\boldsymbol{y}})$ (40)
$\displaystyle\Longleftrightarrow{\boldsymbol{z}}=\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}{\boldsymbol{y}})\quad\text{(Lemma~{}\ref{lemma:conjugate_inverse})}$
(41)

where $g^{*}_{\gamma}$ denotes the convex conjugate of $g_{\gamma}$. Note that
we can write $\nabla$ in the last term instead of $\partial$ for
$g^{*}_{\gamma}$ because this function is $1$-smooth with respect to the
${\boldsymbol{H}}^{-1}$-norm by the strong convexity of $g_{\gamma}$.

The term $\frac{1}{2}\|\cdot\|_{{\boldsymbol{H}}^{-1}}^{2}$ is
$\sigma^{-1}_{\min}$-smooth and so $g_{\gamma}$ is $(\sigma^{-1}_{\min}+\gamma
L_{g})$-smooth. By the duality between Lipschitz gradient and strong convexity
(see e.g., Rockafellar and Wets, (1998)), $g_{\gamma}^{*}$ is
$1/(\sigma^{-1}_{\min}+\gamma L)$-strongly convex. Then for arbitrary
${\boldsymbol{y}}$ and $\widetilde{\boldsymbol{y}}$ we have

$\displaystyle\|\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}{\boldsymbol{y}})-\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}})\|\|{\boldsymbol{H}}^{-1}{\boldsymbol{y}}-{\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}}\|$
(42) $\displaystyle\qquad\geq\langle\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}{\boldsymbol{y}})-\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}}),{\boldsymbol{H}}^{-1}{\boldsymbol{y}}-{\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}}\rangle\qquad\text{
(by Cauchy-Schwarz)}$ (43)
$\displaystyle\qquad\geq\frac{1}{1/\sigma_{\min}+\gamma
L_{g}}\|{\boldsymbol{H}}^{-1}{\boldsymbol{y}}-{\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}}\|^{2}\qquad\text{(by
strong convexity)}$ (44)

The result is trivial when ${\boldsymbol{y}}=\widetilde{{\boldsymbol{y}}}$. We
can then assume ${\boldsymbol{y}}\neq\widetilde{{\boldsymbol{y}}}$, and
dividing both sides by
${\|{\boldsymbol{H}}^{-1}{\boldsymbol{y}}-{\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}}\|}$
(non-zero by assumption) we obtain

$\displaystyle\|\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}{\boldsymbol{y}})-\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}})\|$
$\displaystyle\geq\frac{1}{1/\sigma_{\min}+\gamma
L_{g}}\|{\boldsymbol{H}}^{-1}{\boldsymbol{y}}-{\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}}\|$
(45) $\displaystyle\geq\frac{\sigma^{-1}_{\max}}{1/\sigma_{\min}+\gamma
L_{g}}\|{\boldsymbol{y}}-\widetilde{{\boldsymbol{y}}}\|~{}.$ (46)

where the last inequality we have used that
${\boldsymbol{y}}\to\frac{1}{2}{\boldsymbol{y}}^{T}{\boldsymbol{H}}^{-1}{\boldsymbol{y}}$
is strongly convex with strong convexity parameter $\sigma_{\max}^{-1}$ and
${\boldsymbol{H}}^{-1}{\boldsymbol{y}}$ is the gradient of this function.

Using now the equivalence between $\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}\cdot)$ of Eq. (41) and the proximal
operator we finally have the claimed bound:

$\displaystyle\|{\boldsymbol{y}}-\widetilde{{\boldsymbol{y}}}\|$
$\displaystyle\leq{(\sigma_{\max}/\sigma_{\min}+\gamma
L_{g}\sigma_{\max}})\|\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}{\boldsymbol{y}})-\nabla
g^{*}_{\gamma}({\boldsymbol{H}}^{-1}\widetilde{{\boldsymbol{y}}})\|$ (47)
$\displaystyle={(\sigma_{\max}/\sigma_{\min}+\gamma
L_{g}\sigma_{\max})}\|\operatorname{\mathbf{prox}}_{\gamma
g}({\boldsymbol{y}})-\operatorname{\mathbf{prox}}_{\gamma
g}(\widetilde{{\boldsymbol{y}}})\|~{}.$ (48)

∎

######  Lemma 7 (Block firm non-expansiveness).

Let ${\boldsymbol{x}},\widetilde{{\boldsymbol{x}}}$ be two arbitrary vectors
in ${\mathbb{R}}^{p}$, $g$ be a block-separable convex lower semicontinuous
function with blocks $\mathcal{B}$. Let
${\boldsymbol{z}}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname{\mathbf{prox}}_{\gamma
g}({\boldsymbol{x}})$,
$\widetilde{{\boldsymbol{z}}}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname{\mathbf{prox}}_{\gamma
g}(\widetilde{{\boldsymbol{x}}})$. Then for any subset
$\mathcal{A}\subseteq\mathcal{B}$ it is verified that:

$\langle[{\boldsymbol{z}}-\widetilde{{\boldsymbol{z}}}]_{\mathcal{A}},[{\boldsymbol{x}}-\widetilde{{\boldsymbol{x}}}]_{\mathcal{A}}\rangle\geq\|[{\boldsymbol{z}}-\widetilde{{\boldsymbol{z}}}]_{\mathcal{A}}\|^{2}\,.$
(49)

###### Proof.

By the block-separability of $g$, the proximal operator is the concatenation
of the proximal operators of the blocks. In other words, for any block
$A\in\mathcal{A}$ we have:

$[{\boldsymbol{z}}]_{A}=\operatorname{\mathbf{prox}}_{\gamma
g_{A}}([{\boldsymbol{x}}]_{A})~{},\quad[\widetilde{{\boldsymbol{z}}}]_{A}=\operatorname{\mathbf{prox}}_{\gamma
g_{A}}([\widetilde{{\boldsymbol{x}}}]_{A})\,,$ (50)

where $g_{A}$ is the restriction of $g$ to $A$. By firm non-expansiveness of
the proximal operator (see e.g. Bauschke and Combettes, (2017, Proposition
4.2)) we have that:

$\langle[{\boldsymbol{z}}]_{A}-[\widetilde{{\boldsymbol{z}}}]_{A},[{\boldsymbol{x}}]_{A}-[\widetilde{{\boldsymbol{x}}}]_{A}\rangle\geq\|[{\boldsymbol{z}}]_{A}-[\widetilde{{\boldsymbol{z}}}]_{A}\|^{2}\,.$

Summing over the blocks in $\mathcal{A}$ yields the desired result. ∎

##  Appendix Appendix B Fixed point characterization

In this subsection we provide a characterization of the fixed points of the
three operator splitting.

The following theorem characterizes the set of fixed points
$\boldsymbol{G}_{\gamma}$ (defined in (4)) as the weighted Minkowski sum of
primal and dual solutions. We will denote by
$\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$ the set of fixed points of
$\boldsymbol{G}_{\gamma}$. This characterization seems to be new, and it will
be used in some of the later proofs.

######  Theorem 4 (Fixed point for operator splitting).

Let $\mathcal{P}^{\star}$ denote the set of minimizers of the primal objective
and $\mathcal{D}^{\star}$ the set of minimizers of the dual objective. Then
the set of fixed points of the three splitting is

$\operatorname*{Fix}(\boldsymbol{G}_{\gamma})=\mathcal{P}^{\star}+\gamma\,{\boldsymbol{D}}\mathcal{D}^{\star}=\left\\{{\boldsymbol{x}}+\gamma{\boldsymbol{D}}{\boldsymbol{u}}\,|\,{\boldsymbol{x}}\in\mathcal{P}^{\star},{\boldsymbol{u}}\in\mathcal{D}^{\star}\right\\}~{}.$
(51)

###### Proof.

We first characterize the fixed points of $\boldsymbol{G}_{\gamma}$ by a
subdifferential inclusion. Given ${\boldsymbol{y}}\in{\mathbb{R}}^{p}$, let
${\boldsymbol{z}}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}})$ and
${\boldsymbol{u}}\stackrel{{\scriptstyle\text{def}}}{{=}}{\boldsymbol{D}}^{-1}({\boldsymbol{y}}-{\boldsymbol{z}})/\gamma$.
Consider the following sequence of equivalences:

$\displaystyle{\boldsymbol{y}}=\boldsymbol{G}_{\gamma}({\boldsymbol{y}})$
$\displaystyle\iff\begin{cases}{\boldsymbol{z}}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
g}(2{\boldsymbol{z}}-{\boldsymbol{y}}-\gamma{\boldsymbol{D}}\nabla
f({\boldsymbol{z}}))\\\
{\boldsymbol{z}}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}})\end{cases}$ (by definition of $\boldsymbol{G}_{\gamma}$)
(52)
$\displaystyle\iff\begin{cases}{\boldsymbol{D}}^{-1}(-\frac{\gamma}{\gamma}({\boldsymbol{y}}-{\boldsymbol{z}})-\gamma{\boldsymbol{D}}\nabla
f({\boldsymbol{z}}))/{\gamma}\in\partial g({\boldsymbol{z}})\\\
{\boldsymbol{D}}^{-1}({\boldsymbol{y}}-{\boldsymbol{z}})/{\gamma}\in\partial
h({\boldsymbol{z}})\end{cases}$ (by Lemma 1) (53)
$\displaystyle\iff\begin{cases}-{\boldsymbol{u}}\in\partial(f+g)({\boldsymbol{z}})\\\
{\boldsymbol{u}}\in\partial h({\boldsymbol{z}})\end{cases}$ (54)
$\displaystyle\iff\begin{cases}{\boldsymbol{z}}\in\partial(f+g)^{*}(-{\boldsymbol{u}})\\\
{\boldsymbol{z}}\in\partial h^{*}({\boldsymbol{u}})\end{cases}$ (by Lemma 2)
(55)

The rest of the proof is divided in two parts, proving in the first part that
$\operatorname*{Fix}(\boldsymbol{G}_{\gamma})\subseteq{\mathcal{P}^{\star}+\gamma{\boldsymbol{D}}\mathcal{D}^{\star}}\,$,
and the reverse inclusion in the second part.

_Part 1_. Our goal is to prove
$\operatorname*{Fix}(\boldsymbol{G}_{\gamma})\subseteq{\mathcal{P}^{\star}+\gamma\mathcal{D}^{\star}}\,$.
Let ${\boldsymbol{y}}\in\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$ and
${\boldsymbol{z}}$, ${\boldsymbol{u}}$ be as defined above. From their
definition we immediately have
${\boldsymbol{z}}+\gamma{\boldsymbol{D}}{\boldsymbol{u}}={\boldsymbol{y}}$,
and so we only need to prove that ${\boldsymbol{z}},{\boldsymbol{u}}$ are
minimizers of the primal and dual objective respectively. By definition of
${\boldsymbol{z}}$ we have the following subdifferential inclusions

$\displaystyle{\boldsymbol{z}}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}})$
$\displaystyle\iff\frac{{\boldsymbol{D}}^{-1}}{\gamma}({\boldsymbol{y}}-{\boldsymbol{z}})={\boldsymbol{u}}\in\partial
h({\boldsymbol{z}})$ (56) $\displaystyle\iff{\boldsymbol{z}}\in\partial
h^{*}({\boldsymbol{u}})~{},$ (57)

where we have used Lemma 1 for the first equivalence and Lemma 2 for the
second one. Adding together (56) with the first line of (54), and (57) minus
the first line of (55) gives

$\displaystyle\qquad\,0\in\partial h({\boldsymbol{z}})+\partial
g({\boldsymbol{z}})+\nabla f({\boldsymbol{z}})$ (58) $\displaystyle\text{ and
}~{}0\in\partial
h^{*}({\boldsymbol{u}})-\partial(f+g)^{*}(-{\boldsymbol{u}})~{},$ (59)

and so by the first-order optimality conditions ${\boldsymbol{z}}$ and
${\boldsymbol{u}}$ are minimizers of the primal and dual objectives
respectively. We have proved
$\operatorname*{Fix}(\boldsymbol{G}_{\gamma})\subseteq\mathcal{P}^{\star}+\gamma{\boldsymbol{D}}\mathcal{D}^{\star}\,$.

_Part 2_. Our goal now is to prove the inverse inclusion,
${\mathcal{P}^{\star}+\gamma{\boldsymbol{D}}\mathcal{D}^{\star}}\subseteq\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$.
Let
$({\boldsymbol{x}},{\boldsymbol{u}})\in\mathcal{P}^{\star}\times\mathcal{D}^{\star}$,
we will prove that
${\boldsymbol{y}}\stackrel{{\scriptstyle\text{def}}}{{=}}{\boldsymbol{x}}+\gamma{\boldsymbol{D}}{\boldsymbol{u}}$
is a fixed point of $\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$.

We start by recalling the notion of _paramonotinicity_ , which will play a key
role in this part of the proof. This notion was introduced by Iusem, (1998)
and is key to characterizing the set of fixed points of related methods, such
as the Douglas-Rachford splitting (Bauschke et al.,, 2012). An operator
$\boldsymbol{C}$ is said to be paramonotonic if the following implication is
verified

$\begin{rcases*}\boldsymbol{a}^{\star}\in\boldsymbol{C}\boldsymbol{a}\\\
\boldsymbol{b}^{\star}\in\boldsymbol{C}\boldsymbol{b}\\\
\langle\boldsymbol{a}^{\star}-\boldsymbol{b}^{\star},\boldsymbol{a}-\boldsymbol{b}\rangle=0\end{rcases*}\implies\boldsymbol{a}^{\star}\in\boldsymbol{C}\boldsymbol{b}\text{
and }\boldsymbol{b}^{\star}\in\boldsymbol{C}\boldsymbol{a}~{}.$ (60)

The usefulness of this notion in this case comes from the fact that the
subdifferential of a convex proper lower semicontinuous function is
paramonotonic (Iusem,, 1998, Proposition 2.2). Hence we have that $\partial h$
and $\partial(f+g)$ are paramonotonic.

By the first-order optimality conditions on the primal and dual loss we have
that there exists elements ${\boldsymbol{u}}_{\boldsymbol{z}}$ and
${\boldsymbol{z}}_{\boldsymbol{u}}$ such that

$\displaystyle{\boldsymbol{u}}_{\boldsymbol{z}}\in\partial
h({\boldsymbol{z}})\cap(-\partial(f+g)({\boldsymbol{z}}))$ (61)
$\displaystyle{\boldsymbol{z}}_{\boldsymbol{u}}\in\partial
h^{*}({\boldsymbol{u}})\cap(\partial(f+g)^{*}(-{\boldsymbol{u}}))~{},$ (62)

where the second inclusion can be written equivalently using the conjugate-
inverse identity (Lemma 2) as

${\boldsymbol{u}}\in\partial
h({\boldsymbol{z}}_{\boldsymbol{u}})\cap(-\partial(f+g)({\boldsymbol{z}}_{\boldsymbol{u}}))~{}.$
(63)

Using Eq. (61) and (63) we have by monotony of $\partial h$ and
$\partial(f+g)$

$\langle{\boldsymbol{u}}_{\boldsymbol{z}}-{\boldsymbol{u}},{\boldsymbol{z}}-{\boldsymbol{z}}_{\boldsymbol{u}}\rangle\geq
0~{}\text{ and
}\langle{\boldsymbol{u}}_{\boldsymbol{z}}-{\boldsymbol{u}},{\boldsymbol{z}}-{\boldsymbol{z}}_{\boldsymbol{u}}\rangle\leq
0$ (64)

from where we necessarily have
$\langle{\boldsymbol{u}}_{\boldsymbol{z}}-{\boldsymbol{u}},{\boldsymbol{z}}-{\boldsymbol{z}}_{\boldsymbol{u}}\rangle=0$.
We hence have by paramonotonicity of $\partial h$

$\begin{rcases}{\boldsymbol{u}}_{\boldsymbol{z}}\in\partial
h({\boldsymbol{z}})\\\ {\boldsymbol{u}}\in\partial
h({\boldsymbol{z}}_{\boldsymbol{u}})\\\
\langle{\boldsymbol{u}}_{\boldsymbol{z}}-{\boldsymbol{u}},{\boldsymbol{z}}-{\boldsymbol{z}}_{\boldsymbol{u}}\rangle=0\end{rcases}\implies{\boldsymbol{u}}\in\partial
h({\boldsymbol{z}})$ (65)

Similarly, by paramonotonicity of $\partial(f+g)$ we have

$\begin{rcases}-{\boldsymbol{u}}_{\boldsymbol{z}}\in\partial(f+g)({\boldsymbol{z}})\\\
-{\boldsymbol{u}}\in\partial(f+g)({\boldsymbol{z}}_{\boldsymbol{u}})\\\
\langle{\boldsymbol{u}}_{\boldsymbol{z}}-{\boldsymbol{u}},{\boldsymbol{z}}-{\boldsymbol{z}}_{\boldsymbol{u}}\rangle=0\end{rcases}\implies-{\boldsymbol{u}}\in\partial(f+g)({\boldsymbol{z}})$
(66)

Combining the last two equations we have by the definition of
${\boldsymbol{y}}$ the following inclusions

$\begin{cases}-{\boldsymbol{u}}\in\partial(f+g)({\boldsymbol{z}})\\\
{\boldsymbol{u}}\in\partial h({\boldsymbol{z}})\end{cases}$ (67)

which by Eq. (54) implies that
${\boldsymbol{y}}\in\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$ (note that
these are all equivalences from (52) to (54)). This concludes the proof. ∎

######  Corollary 1 (Minimizer of our objective).

Let ${\boldsymbol{y}}\in\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$. Then we
have that
${\boldsymbol{z}}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}})$ is a minimizer of the primal objective $\mathcal{P}$,
$\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}})=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}(2{\boldsymbol{z}}-{\boldsymbol{y}}-\gamma{\boldsymbol{D}}\nabla
f({\boldsymbol{z}}))$, and
${\boldsymbol{u}}={\boldsymbol{D}}^{-1}({\boldsymbol{y}}-\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}}))/\gamma$ is a minimizer of the dual objective.

###### Proof.

This follows from the first part of the proof of Theorem 4. In particular, Eq.
(58) shows that
${\boldsymbol{z}}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}})$ is a minimizer of the primal objective, while Eq. (59)
shows that
${\boldsymbol{u}}={\boldsymbol{D}}^{-1}({\boldsymbol{y}}-\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}}))/\gamma$ is a minimizer of the dual objective.

The identity $\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}})=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}(2{\boldsymbol{z}}-{\boldsymbol{y}}-\gamma{\boldsymbol{D}}\nabla
f({\boldsymbol{z}}))$ comes from the definition of fixed point (52). ∎

##  Appendix Appendix C Iteration complexity analysis

In this section we provide a proof for the convergence rate analysis of the
proposed methods of §4. We will start by with the proof of linear convergence
(Theorem 3) and then prove the sublinear convergence rate (Theorem 1), as this
last theorem reuses many elements from the first.

Unless explicitly stated (e.g., in Theorem 1), the results are only proven for
the sparse variants. Since the dense variants are a special case of the sparse
variants with ${\boldsymbol{P}}_{i}=\boldsymbol{I}$,
${\boldsymbol{D}}=\boldsymbol{I}$, the results for the dense variants follow
as a special case.

#### Structure of this appendix.

  * •

Appendix C.1 provides technical lemmas that will be used in later proofs.

  * •

Appendix C.2 provides a proof for the linear convergence (under assumptions)
of the proposed methods (Theorem 3).

  * •

Appendix C.3 provides a sublinear convergence rate for the dense variants of
the proposed methods (Theorem 1).

  * •

Appendix C.4 provides a (weaker) sublinear convergence rate for the sparse
variants of the proposed methods (Theorem 2).

#### Extra notation for this section.

  * •

We define
$f_{i}({\boldsymbol{x}})=\psi_{i}({\boldsymbol{x}})+\omega({\boldsymbol{x}})$

  * •

To provide a unified analysis of the dense and sparse algorithm, we define the
following auxiliary function:

$\xi_{i}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\psi_{i}({\boldsymbol{x}})+\sum_{B\in
T_{i}}d_{B}\omega_{B}([{\boldsymbol{x}}]_{B})~{}.$ (68)

Note that $\frac{1}{n}\sum_{i=1}^{n}\xi_{i}=f$. Since $\psi_{i}$ is
$L_{\psi}$-smooth and $\omega$ is $L_{\omega}$-smooth we have that $\xi_{i}$
is $L_{f}$-smooth, with $L_{f}=L_{\psi}+d_{\max}L_{\omega}$ (as defined in
§4).

  * •

Contrary to full gradient algorithms, in stochastic variance reduced methods
the objective function is not guaranteed to decrease at each iteration. To
compensate for this, a common approach is to add a positive term that
decreases throughout the iterations. The resulting function is often called a
_Lyapunov_ function. Throughout this paper, the positive term that we will add
is the following:

$H_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{n}\sum_{i=1}^{n}m_{i,t}$
(69)

where $m_{i,t}$ are positive constants initialized as

$m_{i,0}=\frac{1}{2L_{f}}\|{\boldsymbol{\alpha}}_{i,0}-\nabla\psi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
(70)

and updated at each iteration as

$m_{i,t+1}=\begin{cases}B_{f_{i}}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})&\text{
if ${\boldsymbol{\alpha}}_{i}$ has been updated}\\\ m_{i,t}&\text{
otherwise}\end{cases}~{},$ (71)

for all $i\in\\{1,\ldots,n\\}$. This term is a hybrid between those used by
Defazio et al., (2014) and Hofmann et al., (2015). Like Defazio et al.,
(2014), it will allow us to obtain a large $1/(3L_{f})$ step size, contrary to
the $<1/4L_{f}$ step size of Hofmann et al., (2015). Like Hofmann et al.,
(2015) (and unlike Defazio et al., (2014)), it will allow to initialize
${\boldsymbol{\alpha}}_{0}$ arbitrarily.

  * •

For convenience, we denote by $\langle\cdot,\cdot\rangle_{(i)}$ (resp.
$\|\cdot\|_{(i)}$) the scalar product (resp. norm) restricted to blocks in the
extended support, i.e.,
$\langle{\boldsymbol{x}},{\boldsymbol{y}}\rangle_{(i)}\stackrel{{\scriptstyle\text{def}}}{{=}}\langle{\boldsymbol{x}},{\boldsymbol{P}}_{i}{\boldsymbol{y}}\rangle$
and
$\|{\boldsymbol{x}}\|_{(i)}\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{\langle{\boldsymbol{x}},{\boldsymbol{x}}\rangle_{(i)}}$.

  * •

We denote by $d_{\max}$ the maximum entry in the diagonal matrix
${\boldsymbol{D}}$, with ${\boldsymbol{D}}$ as defined in §2.1.

### Appendix C.1 Preliminaries

In this subsection we state some key lemmas that are used in both the proof of
linear and sublinear convergence.

######  Lemma 8 (Strong convexity inequality).

Let $\psi_{i}$ be $\mu_{\psi}$-strongly convex. Let $\omega$ be
$\mu_{\omega}$-strongly convex (where we allow $\mu_{\psi}=\mu_{\omega}=0$).
Then with $\mu=\mu_{\psi}+\mu_{\omega}$ we have the following inequality for
arbitrary ${\boldsymbol{x}}$ and ${\boldsymbol{y}}$ in the domain:

$\displaystyle f({\boldsymbol{x}})$ $\displaystyle\geq
f({\boldsymbol{y}})+\langle\nabla
f({\boldsymbol{y}}),{\boldsymbol{x}}-{\boldsymbol{y}}\rangle+\frac{1}{2n(L_{f}-\mu)}\sum_{i=1}^{n}\|\nabla\xi_{i}({\boldsymbol{x}})-\nabla\xi_{i}({\boldsymbol{y}})\|^{2}$
(72) $\displaystyle\qquad+\frac{\mu
L_{f}}{2d_{\text{max}}(L_{f}-\mu)}\|{\boldsymbol{x}}-{\boldsymbol{y}}\|^{2}+\frac{\mu}{L_{f}-\mu}\langle\nabla
f({\boldsymbol{x}})-\nabla
f({\boldsymbol{y}}),{\boldsymbol{y}}-{\boldsymbol{x}}\rangle$

###### Proof.

We start by proving that $\xi_{i}$ is $\mu$-strongly convex when restricted to
his support. Let $\boldsymbol{a},\boldsymbol{b}$ be arbitrary vectors in
${\mathbb{R}}^{p}$. Then we have the following sequence of inequalities:

$\displaystyle\langle\nabla$
$\displaystyle\xi_{i}(\boldsymbol{a})-\nabla\xi_{i}(\boldsymbol{b}),\boldsymbol{a}-\boldsymbol{b}\rangle$
(73)
$\displaystyle=\langle\nabla\psi_{i}(\boldsymbol{a})+{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla\omega(\boldsymbol{a})-\nabla\psi_{i}(\boldsymbol{b})-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla\omega(\boldsymbol{b}),\boldsymbol{a}-\boldsymbol{b}\rangle$
(74) (by definition of $\nabla\xi_{i}$)
$\displaystyle=\langle\nabla\psi_{i}(\boldsymbol{a})-\nabla\psi_{i}(\boldsymbol{b}),\boldsymbol{a}-\boldsymbol{b}\rangle+\langle{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla\omega(\boldsymbol{a})-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla\omega(\boldsymbol{b}),\boldsymbol{a}-\boldsymbol{b}\rangle$
(75)
$\displaystyle\geq\mu_{\psi}\|\boldsymbol{a}-\boldsymbol{b}\|^{2}+\langle{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla\omega(\boldsymbol{a})-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla\omega(\boldsymbol{b}),\boldsymbol{a}-\boldsymbol{b}\rangle$
(76) (by strong convexity of $\nabla\psi_{i}$)
$\displaystyle=\mu_{\psi}\|\boldsymbol{a}-\boldsymbol{b}\|^{2}+\textstyle\sum_{B\in\mathcal{B}}d_{B}\langle\nabla\omega_{B}([\boldsymbol{a}]_{B})-\nabla\omega_{B}([\boldsymbol{b}]_{B}),[\boldsymbol{a}]_{B}-[\boldsymbol{b}]_{B}\rangle$
(77) (by block separability of $\omega$ and definition of
${\boldsymbol{P}}_{i}{\boldsymbol{D}}$)
$\displaystyle\geq\mu_{\psi}\|\boldsymbol{a}-\boldsymbol{b}\|+\textstyle\sum_{B\in\mathcal{B}}d_{B}\mu_{\omega}\|[\boldsymbol{a}]_{B}-[\boldsymbol{b}]_{B}\|^{2}$
(78) (strong convexity of $\omega_{B}$, consequence of strong cvx of $\omega$)
$\displaystyle\geq\mu_{\psi}\|\boldsymbol{a}-\boldsymbol{b}\|^{2}+\mu_{\omega}\|\boldsymbol{a}-\boldsymbol{b}\|_{(i)}^{2}$
(79) (using $d_{B}\geq 1$ by definition)
$\displaystyle\geq\underbrace{\left(\mu_{\psi}+\mu_{\omega}\right)}_{=\mu}\|\boldsymbol{a}-\boldsymbol{b}\|_{(i)}^{2}~{}.$
(80)

We have proved that $\xi_{i}$ is $\mu$-strongly convex on the subspace
generated by the extended support (i.e., with respect to the norm
$\|\cdot\|_{(i)}$). Since it is also $L_{f}$-smooth by (68), we can apply
(Defazio et al.,, 2014, Lemma 4) to obtain the following inequality, valid for
all $\boldsymbol{a}$ and $\boldsymbol{b}$ in its domain:

$\displaystyle\xi_{i}(\boldsymbol{a})$
$\displaystyle\geq\xi_{i}(\boldsymbol{b})+\langle\nabla\xi_{i}(\boldsymbol{b}),\boldsymbol{a}-\boldsymbol{b}\rangle+\frac{1}{2(L_{f}-\mu)}\|\nabla\xi_{i}(\boldsymbol{a})-\nabla\xi_{i}(\boldsymbol{b})\|^{2}$
(81) $\displaystyle\qquad+\frac{\mu
L_{f}}{2(L_{f}-\mu)}\|\boldsymbol{a}-\boldsymbol{b}\|_{(i)}^{2}+\frac{\mu}{L_{f}-\mu}\langle\nabla\xi_{i}(\boldsymbol{a})-\nabla\xi_{i}(\boldsymbol{b}),\boldsymbol{b}-\boldsymbol{a}\rangle$

We will apply the previous inequality at $\boldsymbol{a}={\boldsymbol{x}}$,
$\boldsymbol{b}={\boldsymbol{y}}$ and average over all $i$. Note that
$\frac{1}{n}\sum_{i=1}^{n}\xi_{i}({\boldsymbol{x}})=\psi({\boldsymbol{x}})+\omega({\boldsymbol{x}})=f({\boldsymbol{x}})$
by definition of $\xi_{i}$ and so we can write

$\displaystyle f({\boldsymbol{x}})$ $\displaystyle\geq
f({\boldsymbol{y}})+\frac{1}{n}\sum_{i=1}^{n}\langle\nabla\xi_{i}({\boldsymbol{y}}),{\boldsymbol{x}}-{\boldsymbol{y}}\rangle+\frac{1}{2n(L_{f}-\mu)}\sum_{i=1}^{n}\|\nabla\xi_{i}({\boldsymbol{x}})-\nabla\xi_{i}({\boldsymbol{y}})\|^{2}$
$\displaystyle\qquad+\frac{\mu
L_{f}}{2(L_{f}-\mu)}\frac{1}{n}\sum_{i=1}^{n}\|{\boldsymbol{y}}-{\boldsymbol{x}}\|_{(i)}^{2}+\frac{\mu}{L_{f}-\mu}\frac{1}{n}\sum_{i=1}^{n}\langle\nabla\xi_{i}({\boldsymbol{x}})-\nabla\xi_{i}({\boldsymbol{y}}),{\boldsymbol{y}}-{\boldsymbol{x}}\rangle$
(82)

We can simplify the terms in this inequality as follows:

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\langle\nabla\xi_{i}({\boldsymbol{y}}),{\boldsymbol{x}}-{\boldsymbol{y}}\rangle=\langle\frac{1}{n}\sum_{i=1}^{n}\nabla\xi_{i}({\boldsymbol{y}}),{\boldsymbol{x}}-{\boldsymbol{y}}\rangle=\langle\nabla
f({\boldsymbol{y}}),{\boldsymbol{x}}-{\boldsymbol{y}}\rangle$ (83)
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\|{\boldsymbol{x}}-{\boldsymbol{y}}\|_{(i)}^{2}=\|{\boldsymbol{x}}-{\boldsymbol{y}}\|^{2}_{{\boldsymbol{D}}^{-1}}\geq\frac{1}{d_{\text{max}}}\|{\boldsymbol{x}}-{\boldsymbol{y}}\|^{2}$
(84)
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\langle\nabla\xi_{i}({\boldsymbol{x}})-\nabla\xi_{i}({\boldsymbol{y}}),{\boldsymbol{y}}-{\boldsymbol{x}}\rangle=\langle\nabla
f({\boldsymbol{x}})-\nabla
f({\boldsymbol{y}}),{\boldsymbol{y}}-{\boldsymbol{x}}\rangle$ (85)

The second equality results by the definition of ${\boldsymbol{D}}$ which
gives: $\mathbf{E}[\boldsymbol{P}_{i}]={\boldsymbol{D}}^{-1}$.Using the
previous identities (and inequality) into (82) we finally obtain the desired
bound:

$\displaystyle f({\boldsymbol{x}})$ $\displaystyle\geq
f({\boldsymbol{y}})+\langle\nabla
f({\boldsymbol{y}}),{\boldsymbol{x}}-{\boldsymbol{y}}\rangle+\frac{1}{2n(L_{f}-\mu)}\sum_{i=1}^{n}\|\nabla\xi_{i}({\boldsymbol{x}})-\nabla\xi_{i}({\boldsymbol{y}})\|^{2}$
(86) $\displaystyle\qquad+\frac{\mu
L_{f}}{2d_{\text{max}}(L_{f}-\mu)}\|{\boldsymbol{x}}-{\boldsymbol{y}}\|^{2}+\frac{\mu}{L_{f}-\mu}\langle\nabla
f({\boldsymbol{x}})-\nabla
f({\boldsymbol{y}}),{\boldsymbol{y}}-{\boldsymbol{x}}\rangle$

∎

######  Lemma 9 (Bound on gradient estimate variance).

Let $\mathbf{E}$ denote the conditional expectation with respect to the random
index $i$ selected at the $t$-th iteration. Then we have the following
inequality:

$\displaystyle\mathbf{E}\|{\boldsymbol{v}}_{t}-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\leq(1+\beta^{-1})2L_{f}{H}_{t}+(1+\beta)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
(87) $\displaystyle\qquad-\beta\|\nabla f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star})\|^{2}~{},$

valid for any $\beta>0$.

###### Proof.

Let
$\psi\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{n}\sum_{j=1}^{n}\psi_{j}$.
Then we have the following sequence of inequalities:

$\displaystyle\mathbf{E}\|{\boldsymbol{v}}_{t}-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star})\|^{2}=\mathbf{E}\|\underbrace{\nabla\xi_{i}({\boldsymbol{z}}_{t})-{\boldsymbol{\alpha}}_{i,t}+\boldsymbol{D}_{i}\overline{{\boldsymbol{\alpha}}}_{t}-\boldsymbol{D}_{i}\nabla
f({\boldsymbol{x}}^{\star})}_{=\boldsymbol{\zeta}_{i}}\|^{2}$ (88)
$\displaystyle=\mathbf{E}\|\overbrace{\left[-{\boldsymbol{\alpha}}_{i,t}+\boldsymbol{D}_{i}\overline{{\boldsymbol{\alpha}}}_{t}+\nabla\xi_{i}({\boldsymbol{x}}^{\star})-\boldsymbol{D}_{i}\nabla
f({\boldsymbol{x}}^{\star})\right]+\big{[}\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})}^{\boldsymbol{\zeta}_{i}}-\overbrace{(\nabla
f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star}))}^{\mathbf{E}\boldsymbol{\zeta}_{i}}\big{]}\|^{2}$
$\displaystyle\qquad+\|\overbrace{\nabla f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star})}^{\mathbf{E}\boldsymbol{\zeta}_{i}}\|^{2}$ (89)
(by Lemma 3 with ${\boldsymbol{Q}}_{i}=\boldsymbol{I}$,
$\boldsymbol{A}=\boldsymbol{I}$, and where we have also added and substracted
$\nabla\xi_{i}({\boldsymbol{x}}^{\star})$)
$\displaystyle\leq(1+\beta^{-1})\mathbf{E}\|{\boldsymbol{\alpha}}_{i,t}-\psi_{i}({\boldsymbol{x}}^{\star})-\boldsymbol{D}_{i}\overline{{\boldsymbol{\alpha}}}_{t}+\boldsymbol{D}_{i}\nabla\psi({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad+(1+\beta)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})-\nabla
f({\boldsymbol{z}}_{t})+\nabla f({\boldsymbol{x}}^{\star})\|^{2}+\|\nabla
f({\boldsymbol{z}}_{t})-\nabla f({\boldsymbol{x}}^{\star})\|^{2}$ (90) (by
Young’s inequality
$\|\boldsymbol{a}+\boldsymbol{b}\|^{2}\leq(1+\beta^{-1})\|\boldsymbol{a}\|^{2}+(1+\beta)\|\boldsymbol{b}\|^{2}$
and definition of $\xi_{i}$)
$\displaystyle=(1+\beta^{-1})\mathbf{E}\|{\boldsymbol{\alpha}}_{i,t}-\psi_{i}({\boldsymbol{x}}^{\star})\|^{2}-(1+\beta^{-1})\mathbf{E}\|\overline{{\boldsymbol{\alpha}}}_{t}-\nabla\psi({\boldsymbol{x}}^{\star})\|_{{\boldsymbol{D}}}^{2}$
$\displaystyle\qquad+(1+\beta)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}-\beta\|\nabla
f({\boldsymbol{z}}_{t})-\nabla f({\boldsymbol{x}}^{\star})\|^{2}~{},$ (91)

where in the last equivalence we have applied Lemma 3 both to the first term
(with ${\boldsymbol{Q}}_{i}={\boldsymbol{P}}_{i}$,
$\boldsymbol{A}={\boldsymbol{D}}$) and to the second term (this time with
$\boldsymbol{A}=\boldsymbol{I}$, ${\boldsymbol{Q}}_{i}=\boldsymbol{I}$). In
all, and dropping the negative second term we have the inequality

$\displaystyle\mathbf{E}\|{\boldsymbol{v}}_{t}-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\leq(1+\beta^{-1})\mathbf{E}\|{\boldsymbol{\alpha}}_{i,t}-\psi_{i}({\boldsymbol{x}}^{\star})\|^{2}+(1+\beta)\mathbf{E}\|\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
(92) $\displaystyle\qquad-\beta\|\nabla f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star})\|^{2}$

We will now bound the first term of the above inequality. Let $\mathcal{J}$
denote the set of indices for which the memory terms have been updated at
least once and $\mathcal{J}^{c}$ its complement. For $j\in\mathcal{J}$, we
denote by ${\boldsymbol{\phi}}_{j,t}$ the iterate at which
${\boldsymbol{\alpha}}_{j}$ was last updated, i.e.,
${\boldsymbol{\alpha}}_{j,t}=\nabla\psi({\boldsymbol{\phi}}_{j,t})$ for all
$j$ and $t$. Then we have

$\displaystyle\mathbf{E}\|{\boldsymbol{\alpha}}_{i,t}-\nabla\psi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle=\frac{1}{n}\left(\sum_{j\in\mathcal{J}}\|{\boldsymbol{\alpha}}_{i,t}-\nabla\psi_{i}({\boldsymbol{x}}^{\star})\|^{2}+2L_{f}\sum_{j\in\mathcal{J}^{c}}\xi_{t,j}\right)$
(93)
$\displaystyle\leq\frac{2L_{f}}{n}\left(\sum_{j\in\mathcal{J}}B_{\psi_{i}}({\boldsymbol{\phi}}_{i,k},{\boldsymbol{x}}^{\star})+\sum_{j\in\mathcal{J}^{c}}\xi_{t,j}\right)$
(94) (by Lemma 4 and also using $L_{f}\geq L_{\psi}$)
$\displaystyle\leq\frac{2L_{f}}{n}\left(\sum_{j\in\mathcal{J}}B_{f_{i}}({\boldsymbol{\phi}}_{i,k},{\boldsymbol{x}}^{\star})+\sum_{j\in\mathcal{J}^{c}}\xi_{t,j}\right)$
(95) (adding $B_{\omega}$, which is positive by convexity of $\omega$)
$\displaystyle=2L_{f}H_{t}~{}.$ (96)

Finally, plugging this bound back in (92) we obtain the desired inequality:

$\displaystyle\mathbf{E}\|{\boldsymbol{v}}_{t}-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\leq(1+\beta^{-1})2L_{f}H_{t}+(1+\beta)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
(97) $\displaystyle\qquad-\beta\|\nabla f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star})\|^{2}$

∎

######  Lemma 10 (Evolution of $H_{t}$).

Let $\mathbf{E}$ denote the conditional expectation with respect to the random
index $i$ selected at the $t$-th iteration. Then for every iteration $t\geq 0$
we have (with $q=1$ for Saga variants):

$\mathbf{E}H_{t+1}=\frac{q}{n}B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})+\left(1-\frac{q}{n}\right)H_{t}~{}.$
(98)

###### Proof.

By definition of $m_{j,t+1}$ in Eq. (68), for a fixed index $j$ we have:

$\mathbf{E}[m_{j,t+1}]=\frac{q}{n}B_{f_{j}}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{*})+(1-\frac{q}{n})m_{j,t}~{}.$
(99)

Hence averaging over all indices we get

$\displaystyle\mathbf{E}[H_{t+1}]$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{n}\mathbf{E}[m_{j,t+1}]$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{n}\left(\frac{q}{n}B_{f_{j}}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{*})+(1-\frac{q}{n})m_{j,t}\right)$
$\displaystyle=\frac{q}{n}B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{*})+(1-\frac{q}{n})H_{t}$
(100)

∎

### Appendix C.2 Linear convergence: proof of Theorem 3

The proof is structured as follows:

  * •

We start by proving an inequality that relates
$\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}^{\star}\|^{2}$ with
$\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}$, where
${\boldsymbol{y}}^{\star}$ is a fixed point of $\boldsymbol{G}_{\gamma}$. This
inequality will be central in both proofs of linear and sublinear convergence.
We call this the “master recurrence inequality” (Lemma 11),

  * •

As is often the case in variance reduced methods, a recurrence purely in terms
of the iterates as the one in Lemma 11 does not provide the monotonic decrease
required to prove a linear convergence rate. To overcome this, we will make
use of an auxiliary function which is always larger than the suboptimality
criterion and which _does_ verify a monotonic decrease in expectation. This is
often referred to as a _Lyapunov_ function. The Lyapunov function that we will
use is the following:

$V_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}c\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+{H}_{t}\quad,$
(101)

with $H_{t}$ as defined in (69) and ${\boldsymbol{y}}^{\star}$ an arbitrary
fixed point of $\boldsymbol{G}_{\gamma}$.

  * •

Finally, in Theorem 3 we use the decrease of the Lyapunov function prove the
desired rates of convergence.

######  Lemma 11 (Master recurrence inequality).

Let $\\{{\boldsymbol{y}}_{t},{\boldsymbol{x}}_{t},{\boldsymbol{z}}_{t}\\}$ be
the iterates produced by any of the proposed algorithms,
${\boldsymbol{y}}^{\star}\in\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$ and
${\boldsymbol{x}}^{\star}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}}^{\star})$ (with ${\boldsymbol{D}}=\boldsymbol{I}$ for the
dense variants). Then we have the following inequality, valid for all
$\beta>0$ and $s>0$:

$\displaystyle\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}^{\star}\|^{2}$
$\displaystyle\leq\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+(s-1)\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}$
$\displaystyle\qquad+\frac{\gamma^{2}}{s}(1+\beta^{-1})2L_{f}H_{t}$
$\displaystyle\qquad+(\frac{\gamma^{2}}{s}(1+\beta)-\frac{\gamma}{L_{f}})\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad+(-2\frac{\gamma^{2}\beta}{s}\mu-2\gamma\frac{L_{f}-\mu}{L_{f}})B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})$
$\displaystyle\qquad-\frac{\gamma\mu}{d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$
(102)

###### Proof.

Developing the square we have

$\displaystyle\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}^{\star}\|^{2}$
$\displaystyle=\mathbf{E}\|{\boldsymbol{y}}_{t}+({\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t})-{\boldsymbol{y}}^{\star}\|^{2}$
(103)
$\displaystyle=\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}+2\mathbf{E}\langle{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t},{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\rangle~{}.$
(104)

We will now work towards bounding the last term of this expression.

Let $i$ denote the random index selected at iteration $t$. Note that by
definition of ${\boldsymbol{D}}$ we have
$\mathbf{E}[\boldsymbol{P}_{i}]={\boldsymbol{D}}^{-1}$ and so we can write:

$\displaystyle\mathbf{E}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star},\boldsymbol{P}_{i}({\boldsymbol{y}}_{t}-{{\boldsymbol{y}}^{\star}})\rangle$
$\displaystyle=\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star},{\boldsymbol{y}}_{t}-{{\boldsymbol{y}}^{\star}}\rangle_{{\boldsymbol{D}}^{-1}}$
(105)
$\displaystyle\stackrel{{\scriptstyle}}{{\geq}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{D}}^{-1}}^{2}~{},$
(106)

where in the last inequality we have used Eq. (31) with
${\boldsymbol{z}}_{t}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}}_{t})$ and
${\boldsymbol{x}}^{\star}=\operatorname{\mathbf{prox}}_{\gamma
g}^{{\boldsymbol{D}}}({\boldsymbol{y}}^{\star})$. Using once again the
identity $\mathbf{E}[\boldsymbol{P}_{i}]={\boldsymbol{D}}^{-1}$ and noting
that ${\boldsymbol{z}}_{t}$ does not depend on $i$ we have
$\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{D}}^{-1}}^{2}=\mathbf{E}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}_{{\boldsymbol{P}}_{i}}$
and so in all, we have

$\mathbf{E}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star},{\boldsymbol{P}}_{i}({\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star})\rangle\geq\mathbf{E}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{P}}_{i}}^{2}~{}.$
(107)

Furthermore, by the blockwise version of the firm non-expansiveness of the
prox (Lemma 7), from the definition of ${\boldsymbol{x}}_{t}$ in Vr-Tos we
also have the following inequality, with
${\boldsymbol{x}}_{t}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
g}(2{\boldsymbol{z}}_{t}-{\boldsymbol{y}}_{t}-\gamma{\boldsymbol{v}}_{t})$ and
${\boldsymbol{x}}^{\star}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}(2{\boldsymbol{x}}^{\star}-{\boldsymbol{y}}^{\star}-\gamma{\boldsymbol{D}}\nabla
f({{\boldsymbol{x}}^{\star}}))$, where this last equality is a consequence of
Colollary 1:

$\langle
2{\boldsymbol{z}}_{t}-{\boldsymbol{y}}_{t}-\gamma{\boldsymbol{v}}_{t}-2{\boldsymbol{x}}^{\star}+{\boldsymbol{y}}^{\star}+\gamma{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({{\boldsymbol{x}}^{\star}}),{\boldsymbol{P}}_{i}({\boldsymbol{x}}_{t}-{{\boldsymbol{x}}^{\star}})\rangle-\|{\boldsymbol{x}}_{t}-{{\boldsymbol{x}}^{\star}}\|_{{\boldsymbol{P}}_{i}}^{2}\geq
0~{},$ (108)

which taking conditional expectation gives

$\mathbf{E}\langle
2{\boldsymbol{z}}_{t}-{\boldsymbol{y}}_{t}-\gamma{\boldsymbol{v}}_{t}-2{\boldsymbol{x}}^{\star}+{\boldsymbol{y}}^{\star}+\gamma{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({{\boldsymbol{x}}^{\star}}),{\boldsymbol{P}}_{i}({\boldsymbol{x}}_{t}-{{\boldsymbol{x}}^{\star}})\rangle-\mathbf{E}\|{\boldsymbol{x}}_{t}-{{\boldsymbol{x}}^{\star}}\|_{{\boldsymbol{P}}_{i}}^{2}\geq
0~{}.$ (109)

We now have the following sequence of inequalities:

$\displaystyle\mathbf{E}\langle{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t},{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\rangle$
(110)
$\displaystyle\quad=\mathbf{E}\langle{\boldsymbol{x}}_{t}-{\boldsymbol{z}}_{t},{\boldsymbol{P}}_{i}({\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star})\rangle$
(111) (definition of ${\boldsymbol{y}}_{t+1}$)
$\displaystyle\quad=\mathbf{E}\langle{\boldsymbol{x}}_{t}-{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}+{\boldsymbol{x}}^{\star},{\boldsymbol{P}}_{i}({\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star})\rangle$
(112) (adding and substracting ${\boldsymbol{x}}^{\star}$)
$\displaystyle\quad=\langle{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star},{\boldsymbol{P}}_{i}({\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star})\rangle-\mathbf{E}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star},{\boldsymbol{P}}_{i}({\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star})\rangle$
(113)
$\displaystyle\quad\leq\mathbf{E}\langle{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star},{\boldsymbol{P}}_{i}({\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star})\rangle-\mathbf{E}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{P}}_{i}}^{2}\qquad\text{
(by Eq~{}\eqref{eq:strict_nonexpansiveness_h})}$ (114)
$\displaystyle\quad\leq\mathbf{E}\langle
2{\boldsymbol{z}}_{t}-\gamma{\boldsymbol{v}}_{t}-2{\boldsymbol{x}}^{\star}+\gamma{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{P}}_{i}({\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star})\rangle\\!-\\!\mathbf{E}\|{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{P}}_{i}}^{2}\\!-\\!\mathbf{E}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{P}}_{i}}^{2}$
(115) (adding Eq. (109)) $\displaystyle\quad\leq\mathbf{E}\Big{[}\langle
2{\boldsymbol{z}}_{t}-\gamma{\boldsymbol{v}}_{t}-2{\boldsymbol{x}}^{\star}+\gamma{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{P}}_{i}({\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star})\rangle-\|{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{P}}_{i}}^{2}$
(116)
$\displaystyle\qquad\qquad\qquad-\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{P}}_{i}}^{2}\Big{]}\quad\qquad\text{
(by linearity of expectation)}$
$\displaystyle\quad\leq\mathbf{E}\Big{[}-\left(\|{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{P}}_{i}}^{2}-2\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star},{\boldsymbol{P}}_{i}({\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star})\rangle+\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{P}}_{i}}^{2}\right)$
(117)
$\displaystyle\qquad\qquad\qquad+\langle-\gamma{\boldsymbol{v}}_{t}+\gamma{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{P}}_{i}({\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star})\rangle\Big{]}\quad\qquad\text{
(reordering terms)}$
$\displaystyle\quad=\mathbf{E}\left[-\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}-({\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star})\|_{{\boldsymbol{P}}_{i}}^{2}-\gamma\langle{\boldsymbol{v}}_{t}-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{x}}_{t}-{\boldsymbol{x}}^{\star}\rangle\right]$
(118) (completing the square)
$\displaystyle\quad=\mathbf{E}\Big{[}-\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t}\|_{{\boldsymbol{P}}_{i}}^{2}-\langle\gamma{\boldsymbol{v}}_{t}-\gamma{{\boldsymbol{P}}_{i}}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{x}}_{t}-{\boldsymbol{z}}_{t}\rangle$
$\displaystyle\qquad\qquad-\langle\gamma{\boldsymbol{v}}_{t}-\gamma{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\rangle\Big{]}\qquad\text{
(adding and substracting ${\boldsymbol{z}}_{t}$)}$
$\displaystyle\quad\leq\Big{(}\frac{s}{2}-1\Big{)}\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}+\frac{\gamma^{2}}{2s}\mathbf{E}\|{\boldsymbol{v}}_{t}-{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad\qquad-\gamma\langle\nabla f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\rangle~{},$
(119)

where in the last inequality we have used Young’s inequality:
$|\langle{\boldsymbol{a}},\boldsymbol{b}]\rangle|\leq\frac{s}{2}\|{\boldsymbol{a}}\|^{2}+\frac{1}{2s}\|\boldsymbol{b}\|^{2}$
and definition of ${\boldsymbol{y}}_{t+1}$ for the first term and computed the
expectation in the last term.

Replacing this last inequality into (104) we obtain

$\displaystyle\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}^{\star}\|^{2}$
$\displaystyle\leq\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+(s-1)\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}$
(120)
$\displaystyle+\frac{\gamma^{2}}{s}\mathbf{E}\|{\boldsymbol{v}}_{t}-{{\boldsymbol{P}}_{i}}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star})\|^{2}-2\gamma\langle\nabla
f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\rangle$

We will proceed to further bound the second and last terms using previous
results. For the second term, we can use the bound
$\mathbf{E}\|{\boldsymbol{v}}_{t}-{\boldsymbol{P}}_{i}{\boldsymbol{P}}_{i}{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star})\|^{2}\leq(1+\beta^{-1})2L_{f}H_{t}+{(1+\beta)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}}-\beta\|\nabla
f({\boldsymbol{z}}_{t})-\nabla f({\boldsymbol{x}}^{\star})\|^{2}$ from Lemma
9, giving:

$\displaystyle\frac{\gamma^{2}}{2s}\mathbf{E}\|{\boldsymbol{v}}_{t}-{\boldsymbol{D}}\nabla
f({\boldsymbol{x}}^{\star})\|_{(i)}^{2}$
$\displaystyle\leq\frac{\gamma^{2}}{s}(1+\beta^{-1})L_{f}H_{t}+\frac{\gamma^{2}}{2s}(1+\beta)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
(121) $\displaystyle\qquad-\frac{\gamma^{2}\beta}{2s}\|\nabla
f({\boldsymbol{z}}_{t})-\nabla f({\boldsymbol{x}}^{\star})\|^{2}~{},$

The third term can be bounded using the strong convexity inequality of Lemma 8
with
${\boldsymbol{y}}={\boldsymbol{z}}_{t},{\boldsymbol{x}}={\boldsymbol{x}}^{\star}$,
to obtain

$\displaystyle-$ $\displaystyle\gamma\langle\nabla
f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\rangle=\gamma\langle\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\rangle+\gamma\langle\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}^{\star}-{\boldsymbol{z}}_{t}\rangle$
(122) $\displaystyle\leq\gamma\langle\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\rangle+\gamma(\frac{L_{f}-\mu}{L_{f}}(f({\boldsymbol{x}}^{\star})-f({\boldsymbol{z}}_{t}))$
$\displaystyle\qquad-\frac{1}{2L_{f}}\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}-\frac{\mu}{2d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}-\frac{\mu}{L_{f}}\langle\nabla
f({\boldsymbol{x}}^{\star}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\rangle)$
(123)
$\displaystyle\leq-\gamma\frac{L_{f}-\mu}{L_{f}}B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})+\gamma\left(-\frac{1}{2L_{f}}\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}-\frac{\mu}{2d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}\right)$
(124)

Using the bound for these two terms in (119) we have

$\displaystyle\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}^{\star}\|^{2}$
$\displaystyle\leq\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+(s-1)\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}$
(${\boldsymbol{y}}_{t}$ and ${\boldsymbol{y}}^{\star}$ do not depend on $i$)
$\displaystyle\qquad+\frac{\gamma^{2}}{s}(1+\beta^{-1})2L_{f}H_{t}$
$\displaystyle\qquad+(\frac{\gamma^{2}}{s}(1+\beta)-\frac{\gamma}{L})\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad-\frac{\gamma^{2}\beta}{s}\|\nabla
f({\boldsymbol{z}}_{t})-\nabla f({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad-2\gamma\frac{L-\mu}{L}B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})-\frac{\gamma\mu}{d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$
(125)

We now use the bound $-\|\nabla f({\boldsymbol{z}}_{t})-\nabla
f({\boldsymbol{x}}^{\star})\|^{2}\leq-2\mu
B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})$ (Nesterov,, 2004,
Theorem 2.1.10) to obtain:

$\displaystyle\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}^{\star}\|^{2}$
$\displaystyle\leq\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+(s-1)\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}$
$\displaystyle\qquad+\frac{\gamma^{2}}{s}(1+\beta^{-1})2L_{f}H_{t}$
$\displaystyle\qquad+(\frac{\gamma^{2}}{s}(1+\beta)-\frac{\gamma}{L})\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad+(-2\frac{\gamma^{2}\beta}{s}\mu-2\gamma\frac{L-\mu}{L})B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})$
$\displaystyle\qquad-\frac{\gamma\mu}{d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$
(126)

∎

######  Lemma 12 (Lyapunov inequality).

Let $\\{{\boldsymbol{y}}_{t},{\boldsymbol{x}}_{t},{\boldsymbol{z}}_{t}\\}$ be
the iterates produced by any of the proposed algorithms for $t\geq 0$,
${\boldsymbol{y}}^{\star}\in\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$ and
${\boldsymbol{x}}^{\star}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}}^{\star})$. Let the Lyapunov function $V_{t}$ be as defined
in (101). Then we have the following inequality:

$\displaystyle\mathbf{E}V_{t+1}$ $\displaystyle\leq
V_{t}+c(s-1)\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}+\left(\frac{2L_{f}c\gamma^{2}}{s}(1+\beta^{-1})-\frac{q}{n}\right)H_{t}$
(127)
$\displaystyle\qquad+c\left(\frac{\gamma^{2}}{s}(1+\beta)-\frac{\gamma}{L_{f}}\right)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad+\left(-2\frac{c\gamma^{2}\beta}{s}\mu-2c\gamma\frac{L_{f}-\mu}{L_{f}}+\frac{q}{n}\right)B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})-c\frac{\gamma\mu}{d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$

###### Proof.

We will first compute the conditional expectation of the Lyapunov term,
$\mathbf{E}V_{t+1}$. For the first term we can use the bound in Lemma 11 to
obtain

$\displaystyle
c\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}^{\star}\|^{2}$
$\displaystyle\leq
c\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+c(s-1)\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}$
$\displaystyle\qquad+\frac{\gamma^{2}}{s}(1+\beta^{-1})2cL_{f}H_{t}$
$\displaystyle\qquad+c(\frac{\gamma^{2}}{s}(1+\beta)-\frac{\gamma}{L_{f}})\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad+c(-2\frac{\gamma^{2}\beta}{s}\mu-2\gamma\frac{L_{f}-\mu}{L_{f}})B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})$
$\displaystyle\qquad-c\frac{\gamma\mu}{d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$
(128)

Using Lemma 10, the second term of the Lyapunov function gives :

$\displaystyle\mathbf{E}\,{H}_{t+1}=(1-\frac{q}{n})H_{t}+\frac{q}{n}B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})~{}.$
(129)

and so adding both inequalities we have

$\displaystyle\mathbf{E}V_{t+1}$
$\displaystyle\leq\overbrace{c\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+H_{t}}^{V_{t}}+c(s-1)\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}\|^{2}$
$\displaystyle\qquad+(\frac{2L_{f}c\gamma^{2}}{s}(1+\beta^{-1})-\frac{q}{n})H_{t}$
$\displaystyle\qquad+c(\frac{\gamma^{2}}{s}(1+\beta)-\frac{\gamma}{L_{f}})\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad+(-2c\frac{\gamma^{2}\beta}{s}\mu-2c\gamma\frac{L_{f}-\mu}{L_{f}}+\frac{q}{n})B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})$
$\displaystyle\qquad-c\frac{\gamma\mu}{d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$
(130)

which completes the proof. ∎

######  Theorem 3.

Let $\psi_{i}$ be $\mu_{\psi}$-strongly convex and $\omega$ be
$\mu_{\omega}$-strongly convex, where $\mu_{\psi}+\mu_{\omega}>0$.
Furthermore, let $h$ be $L_{h}$-smooth. Then for any step size
$\gamma\leq{1}/{(3L_{f})}$, all the proposed methods converge geometrically in
expectation. For $\gamma={1}/{(3L_{f})}$, we have the following bound for
Algorithm 1 ($d_{\max{}}=1$ in this case) and Algorithm 2:
${\mathbb{E}}\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}\leq\left(1-\min\Big{\\{}\frac{q}{4n},\frac{1}{3d^{3}_{\max}\delta^{2}\kappa}\Big{\\}}\right)^{t}D_{0}\quad,$
(131) with
$D_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}{d_{\max}}\left[\frac{q}{2\gamma(1-\gamma\mu)n}\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}^{\star}\|^{2}+H_{0}\right]$,
$\delta=(1+L_{h}/(3L_{f}))$, $\kappa=L_{f}/\mu$ and
${\boldsymbol{y}}^{\star}\in\operatorname*{Fix}(\boldsymbol{G}_{\gamma})$.

###### Proof.

From the Lyapunov inequality of Lemma 12 with $s=1$ we have the following
sequence of inequalities

$\displaystyle\mathbf{E}V_{t+1}-(1-\rho)V_{t}\leq$ $\displaystyle\qquad\rho
V_{t}+\left({2Lc\gamma^{2}}(1+\beta^{-1})-\frac{q}{n}\right)H_{t}+c\left({\gamma^{2}}(1+\beta)-\frac{\gamma}{L_{f}}\right)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\quad+\left(-2{c\gamma^{2}\beta}\mu-2c\gamma\frac{L_{f}-\mu}{L_{f}}+\frac{q}{n}\right)B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})-c\frac{\gamma\mu}{d_{\max}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}$
(132)
$\displaystyle\leq\rho\left(c\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+H_{t}\right)+\left(\frac{2L_{f}c\gamma^{2}}{s}(1+\beta^{-1})-\frac{q}{n}\right)H_{t}$
$\displaystyle\qquad+c(\frac{\gamma^{2}}{s}(1+\beta)-\frac{\gamma}{L_{f}})\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}+(-2\frac{c\gamma^{2}\beta}{s}\mu-2c\gamma\frac{L_{f}-\mu}{L_{f}}+\frac{q}{n})B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})$
$\displaystyle\qquad\qquad-c\frac{\gamma\mu}{d_{\max}^{3}(1+\gamma
L_{h})^{2}}\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}$ (133) (using
Lemma 6 on the last term, where we have bounded $d_{\min}\geq 1$)
$\displaystyle\leq c\left[\rho-\frac{\gamma\mu}{d_{\max}^{3}(1+\gamma
L_{h})^{2}}\right]\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}^{\star}\|^{2}+\left[\rho+2L_{f}c\gamma^{2}(1+\beta^{-1})-\frac{q}{n}\right]H_{t}$
$\displaystyle\qquad+c\gamma\left[\gamma(1+\beta)-\frac{1}{L_{f}}\right]\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}$
$\displaystyle\qquad\qquad+\left[-2c\gamma^{2}\beta\mu-2c\gamma\frac{L_{f}-\mu}{L_{f}}+\frac{q}{n}\right]B_{f}({\boldsymbol{z}}_{t},{\boldsymbol{x}}^{\star})$
(134)

It is worth noting that Eq. (133) is the only part of the proof in which we
use the smoothness of $h$.

Taking the coefficients

$c=\frac{q}{2\gamma(1-\gamma\mu)n}~{},\quad\beta=2~{},\quad\rho=\min\left\\{\frac{q}{4n},\frac{1}{3d_{\max}\delta^{2}\kappa}\right\\}~{},$
(135)

With $\delta=d_{\max}(1+\frac{L_{h}}{3L_{f}})$ and $\kappa={L_{f}}/{\mu}$. One
can verify that all square brackets are non-positive for $\gamma\leq
1/(3L_{f})$ (the coefficients are the same, except for the first square
bracket, than those that appear in (Defazio et al.,, 2014, Theorem 1). We
hence have

$\mathbf{E}V_{t+1}\leq(1-\rho)V_{t}~{},$ (136)

which chaining expectations gives

${\mathbb{E}}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}^{\star}\|^{2}\leq{\mathbb{E}}V_{t+1}\leq(1-\rho)^{t+1}V_{0}~{}.$
(137)

This gives a geometric convergence on ${\boldsymbol{y}}_{t}$. However, we
would like to have a convergence rate in terms of the primal iterate
${\boldsymbol{x}}_{t}$.

By Theorem 4 we have that
${\boldsymbol{x}}^{\star}=\operatorname{\mathbf{prox}}^{{\boldsymbol{D}}^{-1}}_{\gamma
h}({\boldsymbol{y}}^{\star})$, and in this case the minimizer is unique
because of strong convexity. Then by firm nonexpansiveness of the prox (Lemma
6) we have

$\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{D}}}^{2}\leq\|{\boldsymbol{y}}_{t+1}-{{\boldsymbol{y}}^{\star}}\|^{2}_{{\boldsymbol{D}}}$
(138)

which combined with Lemma 5 and bounding $d_{\min}$ by $1$ (by definition all
diagonal entries in ${\boldsymbol{D}}$ are $\geq 1$) gives

$\displaystyle\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}$
$\displaystyle\leq\frac{1}{d_{\min}}\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|_{{\boldsymbol{D}}}^{2}\leq\frac{1}{d_{\min}}\|{\boldsymbol{y}}_{t+1}-{{\boldsymbol{y}}^{\star}}\|^{2}_{{\boldsymbol{D}}}$
(139)
$\displaystyle\leq\left(\frac{d_{\max}}{d_{\min}}\right)\|{\boldsymbol{y}}_{t+1}-{{\boldsymbol{y}}^{\star}}\|^{2}\leq{d_{\max}}\|{\boldsymbol{y}}_{t+1}-{{\boldsymbol{y}}^{\star}}\|^{2}$
(140)

Combining this with (137) gives the following bound in
${\boldsymbol{z}}_{t+1}$:

${\mathbb{E}}\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}\leq{d_{\max}}\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{{\boldsymbol{y}}^{\star}}\|^{2}\leq(1-\rho)^{t+1}{d_{\max}}V_{0}~{},$
(141)

and the claimed bound follows from definition of $V_{0}$.

∎

### Appendix C.3 Proof of sublinear convergence rate – dense algorithms

In this section we give a proof of convergence for the dense variants of the
proposed algorithms (Algorithm 1). Because we will not be considering the
sparse variants, we assume ${\boldsymbol{D}}=\boldsymbol{I}$ without explicit
mention.

######  Lemma 13 (Bound on gradient estimate variance, Variant 2).

$\mathbf{E}\|{\boldsymbol{v}}_{t}-\nabla
f({\boldsymbol{z}}_{t})\|^{2}\leq(1+\eta)\mathbf{E}\|\nabla
f_{i}({\boldsymbol{z}}_{t})-\nabla
f_{i}({\boldsymbol{x}}^{\star})\|^{2}+2(1+\eta^{-1})L_{f}H_{t}$ (142)

###### Proof.

$\displaystyle\mathbf{E}\|{\boldsymbol{v}}_{t}-\nabla
f({\boldsymbol{z}}_{t})\|^{2}$ (143)
$\displaystyle=\mathbf{E}\|\nabla\psi_{i}({\boldsymbol{z}}_{t})-{\boldsymbol{\alpha}}_{i,t}+(\overline{{\boldsymbol{\alpha}}}_{t}+\nabla\omega({\boldsymbol{z}}_{t}))-\nabla
f({\boldsymbol{z}}_{t})\|^{2}$ (144)
$\displaystyle=\mathbf{E}\|\underbrace{\nabla
f_{i}({\boldsymbol{z}}_{t})-{\boldsymbol{\alpha}}_{i,t}}_{\xi}+\underbrace{\overline{{\boldsymbol{\alpha}}}_{t}-\nabla
f({\boldsymbol{z}}_{t})}_{-\mathbf{E}\xi}\|^{2}$ (145) (By definition of
$f_{i}$) $\displaystyle\leq\mathbf{E}\|\nabla
f_{i}({\boldsymbol{\phi}}_{i}^{t})-\nabla f_{i}({\boldsymbol{z}}_{t})\|^{2}$
(146) (Applying Lemma 3 and by definition of $f_{i}$)
$\displaystyle\leq\mathbf{E}\|\nabla f_{i}({\boldsymbol{\phi}}_{i}^{t})-\nabla
f_{i}({\boldsymbol{x}}^{\star})+\nabla f_{i}({\boldsymbol{x}}^{\star})-\nabla
f_{i}({\boldsymbol{z}}_{t})\|^{2}$ (147) (Adding and substracting $\nabla
f_{i}({\boldsymbol{x}}^{\star})$)
$\displaystyle\leq(1+\eta^{-1})\mathbf{E}\|\nabla
f_{i}({\boldsymbol{\phi}}_{i}^{t})-\nabla
f_{i}({\boldsymbol{x}}^{\star})\|^{2}+(1+\eta)\mathbf{E}\|\nabla
f_{i}({\boldsymbol{z}}_{t})-\nabla f_{i}({\boldsymbol{x}}^{\star})\|^{2}$
(148) (Applying Young’s inequality) $\displaystyle\leq
2L_{f}(1+\eta^{-1})H_{t}+(1+\eta)\mathbf{E}\|\nabla
f_{i}({\boldsymbol{z}}_{t})-\nabla f_{i}({\boldsymbol{x}}^{\star})\|^{2}$
(149) (Applying Lemma 6 from (Defazio et al.,, 2014) on the first term)

∎

######  Lemma 14 (Saddle point recursive inequality).

Let $\gamma\leq 1/L$ and
${\boldsymbol{y}}_{t},{\boldsymbol{x}}_{t},{\boldsymbol{u}}_{t}$ be the
iterates generated by either Vr-Tos (Algorithm 1). Then we have the following
inequality for any
$({\boldsymbol{x}},{\boldsymbol{u}})\in\operatorname*{dom}(g)\times\operatorname*{dom}(h)$,
with ${\boldsymbol{y}}={\boldsymbol{x}}+\gamma{\boldsymbol{u}}$:

$\displaystyle
2\gamma\mathbf{E}(\mathcal{L}({\boldsymbol{x}}_{t},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t}))+\mathbf{E}\|{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}\|^{2}$
(150)
$\displaystyle\qquad\leq\mathbf{E}\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}\|^{2}+2\gamma^{2}(1+\eta)\mathbf{E}\|\nabla\xi_{i}({\boldsymbol{z}}_{t})-\nabla\xi_{i}({\boldsymbol{x}}^{\star})\|^{2}+4\gamma^{2}(1+\eta^{-1})L_{f}H_{t}$

###### Proof.

By the convexity and the $L$-smoothness inequality, $f$ verifies the following
inequalities for an arbitrary ${\boldsymbol{x}}$:

$\displaystyle f({\boldsymbol{z}}_{t})-f({\boldsymbol{x}})\leq\langle\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\rangle$ (151)
$\displaystyle
f({\boldsymbol{x}}_{t})-f({\boldsymbol{z}}_{t})\leq\langle\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t}-{\boldsymbol{z}}_{t}\rangle+\frac{L}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t}\|^{2}~{}$
(152) $\displaystyle
f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}})\leq\langle\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle+\frac{L}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t}\|^{2}~{},$
(153)

where the last equation is derived from adding the previous two. We now derive
inequalities for $g$ and $h^{*}$. From the subdifferential characterization of
the proximal operator (Lemma 1), the update
${\boldsymbol{z}}_{t}=\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{y}}_{t})$ implies the inclusion

${\boldsymbol{u}}_{t}=\frac{1}{\gamma}({\boldsymbol{y}}_{t}-{\boldsymbol{z}}_{t})\in\partial
h({\boldsymbol{z}}_{t})~{}\implies{\boldsymbol{z}}_{t}\in\partial
h^{*}({\boldsymbol{u}}_{t})$ (154)

where the implication is a consequence of the Fenchel-Young inequality, see
e.g. (Bauschke and Combettes,, 2017, Proposition 16.10) or (Rockafellar and
Wets,, 1998, Proposition 11.3). Similarly, the update
${\boldsymbol{x}}_{t}=\operatorname{\mathbf{prox}}_{\gamma
g}(2{\boldsymbol{z}}_{t}-{\boldsymbol{y}}_{t}-\gamma{\boldsymbol{v}}_{t})$ in
its turn gives the inclusion

$\displaystyle\frac{1}{\gamma}(2{\boldsymbol{z}}_{t}-{\boldsymbol{y}}_{t}-\gamma{\boldsymbol{v}}_{t}-{\boldsymbol{x}}_{t})\in\partial
g({\boldsymbol{x}}_{t})$ (155)

By convexity of $g$ and $h^{*}$ we then have the inequalities

$\displaystyle h^{*}({\boldsymbol{u}}_{t})-h^{*}({\boldsymbol{u}})$
$\displaystyle\leq\langle{\boldsymbol{z}}_{t},{\boldsymbol{u}}_{t}-{\boldsymbol{u}}\rangle~{}.$
(156) $\displaystyle g({\boldsymbol{x}}_{t})-g({\boldsymbol{x}})$
$\displaystyle\leq\frac{1}{\gamma}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t},{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle-\langle{\boldsymbol{u}}_{t}+{\boldsymbol{v}}_{t},{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle$
(157)

Adding (153) and (157) we obtain

$\displaystyle
f({\boldsymbol{x}}_{t})+g({\boldsymbol{x}}_{t})-f({\boldsymbol{x}})-g({\boldsymbol{x}})$
$\displaystyle\leq\frac{1}{\gamma}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t},{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle+\frac{L}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t}\|^{2}-\langle{\boldsymbol{u}}_{t},{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle$
$\displaystyle\qquad+\langle{\boldsymbol{v}}_{t}-\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle$ (158)

Using these, we can now write the following sequence of inequalities for the
Lagrangian suboptimality

$\displaystyle\mathcal{L}({\boldsymbol{x}}_{t},{\boldsymbol{u}}_{t})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t})$
$\displaystyle=f({\boldsymbol{x}}_{t})-f({\boldsymbol{x}})+g({\boldsymbol{x}}_{t})-g({\boldsymbol{x}})+\langle{\boldsymbol{x}}_{t}-{\boldsymbol{x}},{\boldsymbol{u}}_{t}\rangle$
(159)
$\displaystyle\stackrel{{\scriptstyle\eqref{eq:primal_subopt0}}}{{\leq}}\frac{1}{\gamma}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t},{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle+\frac{L}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t}\|^{2}+\langle{\boldsymbol{v}}_{t}-\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle$ (160)
$\displaystyle\mathcal{L}({\boldsymbol{x}}_{t},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}}_{t},{\boldsymbol{u}}_{t})$
$\displaystyle=h^{*}({\boldsymbol{u}}_{t})-h^{*}({\boldsymbol{u}})+\langle{\boldsymbol{x}}_{t},{\boldsymbol{u}}-{\boldsymbol{u}}_{t}\rangle\stackrel{{\scriptstyle\eqref{eq:ineq_g}}}{{\leq}}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t},{\boldsymbol{u}}_{t}-{\boldsymbol{u}}\rangle$
(161)

Adding these two last equations we have

$\displaystyle\mathcal{L}({\boldsymbol{x}}_{t},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t})\leq\frac{1}{\gamma}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t},{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle+\frac{L}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t}\|^{2}+\langle{\boldsymbol{v}}_{t}-\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle$
$\displaystyle\qquad+\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t},{\boldsymbol{u}}_{t}-{\boldsymbol{u}}\rangle$
(162)
$\displaystyle~{}=\frac{1}{\gamma}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t},{\boldsymbol{x}}_{t}+\gamma{\boldsymbol{u}}_{t}-{\boldsymbol{x}}-\gamma{\boldsymbol{u}}\rangle+\frac{L}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t}\|^{2}+\langle{\boldsymbol{v}}_{t}-\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle$ (163)
$\displaystyle~{}=\frac{1}{\gamma}\langle\boldsymbol{{\boldsymbol{y}}}_{t}-{\boldsymbol{y}}_{t+1},{\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}\rangle+\frac{L}{2}\|{\boldsymbol{y}}_{t}-{\boldsymbol{y}}_{t+1}\|^{2}+\langle{\boldsymbol{v}}_{t}-\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t}-{\boldsymbol{x}}\rangle$ (164)
(using
${\boldsymbol{y}}_{t+1}-{\boldsymbol{y}}_{t}={\boldsymbol{x}}_{t}-{\boldsymbol{z}}_{t}$)