# Adaptive Three Operator Splitting

Fabian Pedregosa  Gauthier Gidel

###### Abstract

We propose and analyze an adaptive step-size variant of the Davis-Yin three
operator splitting. This method can solve optimization problems composed by a
sum of a smooth term for which we have access to its gradient and an arbitrary
number of potentially non-smooth terms for which we have access to their
proximal operator. The proposed method sets the step-size based on local
information of the objective –hence allowing for larger step-sizes–, only
requires two extra function evaluations per iteration and does not depend on
any step-size hyperparameter besides an initial estimate. We provide an
iteration complexity analysis that matches the best known results for the non-
adaptive variant: sublinear convergence for general convex functions and
linear convergence under strong convexity of the smooth term and smoothness of
one of the proximal terms. Finally, an empirical comparison with related
methods on 6 different problems illustrates the computational advantage of the
proposed method.

optimization, proximal splitting, line search, three operator splitting

##  1 Introduction

Minimizing the sum of a smooth and a non-smooth term is at the core of many
optimization problems that arise in machine learning and signal processing
(Rudin et al., 1992; Candès et al., 2006; Chambolle & Pock, 2016). In a few
but important cases, such as $\ell_{1}$ or group lasso regularization, the
non-smooth term is simple enough so that its proximal operator is available in
closed form or at least fast to compute. In this case, highly scalable methods
such as proximal gradient descent (Beck & Teboulle, 2009; Nesterov et al.,
2013) or proximal coordinate descent (Richtárik & Takáč, 2014) have shown
state of the art performance. However, the desire to model increasingly
complex phenomena has led to the development of a flurry of penalties with
costly to compute proximal operator. Examples are the overlapping group lasso
(Jacob et al., 2009), multidimensional total variation (Barbero & Sra, 2014)
or trend filtering (Kim et al., 2009), to name a few.

A key observation is that, despite the difficulty in computing its proximal
operator, many of these penalties can be decomposed as a sum of terms for
which we have access to their proximal operator. Proximal splitting methods
like the three operator splitting (Davis & Yin, 2017) offer a principled way
to incorporate these penalties into the optimizer. In this work we will
describe a method to solve optimization problems of the form

$\operatorname*{minimize}_{{\boldsymbol{x}}\in{\mathbb{R}}^{p}}\,\vphantom{\sum_{i}^{n}}f({\boldsymbol{x}})+g({\boldsymbol{x}})+h({\boldsymbol{x}})~{},$
(OPT)

where $f$ is convex and $L_{f}$-smooth (i.e., differentiable with
$L_{f}$-Lipschitz gradient) and $g,h$ are both convex but potentially non-
smooth. We further assume $g$ and $h$ are _proximal_ , i.e., we have access to
the proximal operator.

This formulation allows to express a broad range of problems arising in
machine learning and signal processing: the smooth term includes the least
squares or logistic loss functions; the two proximal terms can be extended to
an arbitrary number via a product space formulation and as we will see in §4.1
include many important penalties such as the group lasso with overlap, total
variation, $\ell_{1}$ trend filtering, etc. Furthermore, the penalties can be
extended-valued, thus allowing an intersection for convex constraints through
the use of the indicator function.

The three operator splitting (TOS) method (Davis & Yin, 2017) is a recently
proposed method for problems of the form (OPT). At each iteration, it only
requires to evaluate once the gradient of $f$ and the proximal operator of $g$
and $h$. It also relies on one step-size parameter, and while it can be set
based on the Lipschitz constant of the gradient of $f$, this is not entirely
satisfactory for two reasons. First, this constant is often costly to compute.
Second, this constant is a global upper bound on the Lipschitz constant, while
locally the Lipschitz constant might be smaller, allowing for larger step-
sizes.

_Adaptive step-size_ methods, also known as inexact and backtracking line
search, instead choose the step-size by verifying a sufficient decrease
condition at each iteration. This allows to take larger step-sizes and has
proven to be an important ingredient in the practical implementation of first
and second-order methods (Nocedal & Wright, 2006).

#### Outline and main contributions.

Our main contribution is the development and analysis of an adaptive variant
of the TOS algorithm. The proposed algorithm does not depend on any step-size
hyperparameter (besides an initial estimate) and enjoys similar convergence
guarantees as the non adaptive variant. The paper is organized as follows:

  * •

_Methods_. §2 describes the proposed algorithm, extended in §2.1 to an
arbitrary number of proximal terms.

  * •

_Analysis_. §3 provides a convergence analysis based on an interpretation of
the algorithm as a saddle-point optimization method. This significantly
departs from the analysis of Davis & Yin (2017) for the non adaptive variant
and results in improved and more general rates.

  * •

_Applications_. §4 discusses the application to different penalties and
presents an empirical comparison on 6 different problems and 5 different
penalties.

#### Notation.

We denote vectors with boldface lower case letters (i.e., ${\boldsymbol{x}}$),
and matrices and vector-valued functions in boldface upper case (i.e.,
${\boldsymbol{X}}$, ${\boldsymbol{T}}(\cdot)$). $\|\cdot\|$ denotes the
euclidean vector norm. Given a matrix
${\boldsymbol{X}}\in{\mathbb{R}}^{n\times p}$, we denote by
$\overline{{\boldsymbol{X}}}$ the average along rows, that is,
$\overline{{\boldsymbol{X}}}=\nicefrac{{1}}{{n}}\sum_{i=1}^{n}{\boldsymbol{X}}_{i}$.
We make extensive use of the proximal operator, defined for a convex function
$\varphi$ and $\gamma>0$ as

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

The domain of a function
$f\mathrel{\mathop{\ordinarycolon}}{\mathbb{R}}^{p}\to]-\infty,\infty]$ is
$\operatorname*{dom}f\stackrel{{\scriptstyle\text{def}}}{{=}}\\{{\boldsymbol{x}}\in{\mathbb{R}}^{p}\,|\,f({\boldsymbol{x}})<\infty\\}$.
The indicator function is denoted $\imath\\{\text{condition}\\}$, which is $0$
if condition is verified and $+\infty$ otherwise. Basic properties and
definitions of convex functions are provided for convenience in Appendix A.

###  1.1 Related work

Proximal splitting methods that can solve problems involving a sum of terms by
accessing the proximal operators of their constituents can be traced back to
the 1970s in the works of Glowinski & Marroco (1975); Gabay & Mercier (1976);
Lions & Mercier (1979). There has been a surge in interest in these methods in
the last years due to their applicability in machine learning (Parikh & Boyd,
2013), signal processing (Combettes & Pesquet, 2011) and parallel optimization
(Boyd et al., 2011).

Algorithms to solve problems 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 generalized forward-backward splitting (Raguet et al., 2013),
the three operator splitting (TOS) (Davis & Yin, 2017), the primal-dual hybrid
gradient (PDHG) method, proposed in (Condat, 2013b; Vũ, 2013) and analyzed by
Chambolle & Pock (2015) and the very recent primal-dual three operator
splitting (Yan, 2018). We note that the last two methods can optimize a more
general objective function in which $h({\boldsymbol{x}})$ is replaced with
$h(\boldsymbol{K}{\boldsymbol{x}})$ for an arbitrary matrix $\boldsymbol{K}$.
The original formulation of these methods requires to set the step-size based
on criteria such as the Lipschitz constant of the gradient of the smooth term,
but variants with adaptive step-size have recently emerged.

An adaptive step-size variant of the PDHG algorithm has recently been proposed
by Malitsky & Pock (2018, §5). Compared to the proposed method, it requires
one less function evaluation per iteration but since the original algorithm
has two step-sizes, it still relies on one step-size hyperparameter.
Convergence rates are not derived.

A different adaptive step-size strategy was proposed by Giselsson et al.
(2016) as a general scheme for averaged operators. TOS is averaged for step-
sizes $<2/L_{f}$, and we denote the combination of both methods LSAO-TOS. An
$\mathcal{O}(1/\sqrt{t})$ convergence rate in terms of the operator residual
norm is derived. Unfortunately, this quantity is difficult to relate to the
more common objective function suboptimality used in the other contributions.

Another adaptive step-size variant of TOS was proposed without proof in the
technical report Davis & Yin (2015, Algorithm 3). It uses the same sufficient
decrease inequality as our method, although the iterates are defined
differently. We found the algorithm sometimes non-convergent and did not
consider it further.

In contrast, we provide a convergence analysis for our method that achieves a
$\mathcal{O}(1/t)$ convergence rate for the ergodic (i.e., averaged) iterate,
and linear convergence under stronger assumptions, matching and in some cases
even improving the best known rates of the non adaptive variant.

| Method | Adaptive | Sublinear rate | Linear rate
---|---|---|---|---
| Adaptive TOS |  |  |
|  (_this work_) | ✓ | ✓ | ✓
TOS |  |  |
(Davis & Yin, 2017) | ✗ | ✓ | ✓
LSAO-TOS |  |  |
(Giselsson et al., 2016) | ✓ |  ✓111Convergence rate in terms of operator residuals. |  ✗
PDHG | ✗ | ✓ | ✗
(Condat, 2013b; Vũ, 2013)
PDHG-LS |  |  |
(Malitsky & Pock, 2018) | ✓ | ✗ | ✗

##  2 Methods

In this section we present our main contribution, a three operator splitting
method with adaptive step-size. The method is detailed in Algorithm 1 and
requires at each iteration to evaluate once the gradient of $f$ and the
proximal operators of $g$ and $h$, and perform two function evaluations of
$f$. At iteration $t$ the candidate step-size $\gamma_{t}$ is chosen as to
verify the following _sufficient decrease_ condition between the iterates
${\boldsymbol{z}}_{t}$ and ${\boldsymbol{x}}_{t+1}$ (Line 1):

$\displaystyle f({\boldsymbol{x}}_{t+1})\leq
Q_{t}({\boldsymbol{x}}_{t+1},\gamma_{t})\,,\,\text{ with \,$Q_{t}$ defined
as}$ $\displaystyle
Q_{t}({\boldsymbol{x}},\gamma)\\!\stackrel{{\scriptstyle\text{def}}}{{=}}\\!f({\boldsymbol{z}}_{t})+\langle\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}\\!-\\!{\boldsymbol{z}}_{t}\rangle\\!+\\!\frac{1}{2\gamma}\|{\boldsymbol{x}}\\!-\\!{\boldsymbol{z}}_{t}\|^{2}.$

This inequality can be interpreted as a quadratic upper bound condition on $f$
at ${\boldsymbol{x}}_{t+1}$: the right-hand side is a quadratic $Q_{t}$ which
is tangent to $f$ at ${\boldsymbol{z}}_{t}$ with amplitude
$(2\gamma_{t})^{-1}$, and both sides are evaluated at
${\boldsymbol{x}}_{t+1}$, defined in Line 1. The under-

${\boldsymbol{x}}$${\boldsymbol{z}}_{t}$${\boldsymbol{x}}_{t+1}$pt$Q_{t}({\boldsymbol{x}})$$f({\boldsymbol{x}})$

lying principle of choosing the step-size based on the minimization of a
quadratic upper bound has already been successful for the proximal-gradient
method, where it is also referred to as backtracking (Beck & Teboulle, 2009)
or full relaxation (Nesterov et al., 2013). In fact, the proposed method
coincides with the aforementioned when one of the proximal terms is constant.

By the properties of $L_{f}$-smooth functions, the sufficient decrease
condition is verified for any $\gamma_{t}\leq 1/L_{f}$. Hence the step-size
search loop always has a finite terminationand the step-size is lower bounded
by $\gamma_{t}\geq\min\\{\tau/L_{f},\gamma_{0}\\}$. The practical advantage of
this strategy is that it allows to consider a step-size potentially larger
than $1/L_{f}$ and verify whether the above is verified at each iteration. If
it is, then the algorithm uses the current step-size, and if not, it decreases
the step-size by a factor which we denote $\tau$.

#### Growing step-size strategies.

We consider two different strategies to initialize next iterate step-size. The
first strategy (_Variant 1_) is the simplest and consists in initializing the
next step-size with the current one (Line 1). In this variant, the step-size
is only allowed to decrease.

The second strategy (_Variant 2_) allows the step-size to increase but in
exchange requires the proximal term $h$ to be Lipschitz continuous (note, not
smooth as $f$ but only Lipschitz). This is the case of most penalties (i.e.,
$\ell_{1}$, group lasso, total variation, etc.) but not of indicator functions
and so is less general than the first variant. As we will see in the
applications section, the ability to grow the step-size has an important
effect on its empirical performance.

Input: ${\boldsymbol{z}}_{0}\in{\mathbb{R}}^{p}$,
${\boldsymbol{u}}_{0}\in{\mathbb{R}}^{p}$, $\gamma_{0}>0,\tau\in(0,1)$

1

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

3

4 repeat $\triangleright$ step-size search loop

5

6
${\boldsymbol{x}}_{t+1}=\operatorname{\mathbf{prox}}_{\gamma_{t}g}({\boldsymbol{z}}_{t}-\gamma_{t}{\boldsymbol{u}}_{t}-\gamma_{t}\nabla
f({\boldsymbol{z}}_{t}))$

7 if  _$f({\boldsymbol{x}}_{t+1})\leq
Q_{t}({\boldsymbol{x}}_{t+1},\gamma_{t})~{}$_ then

8

9 break $\triangleright$ sufficient decrease verified

10 else

11 $\gamma_{t}=\tau\gamma_{t}$ $\triangleright$ decrease step-size

12

13

14
${\boldsymbol{z}}_{t+1}=\operatorname{\mathbf{prox}}_{\gamma_{t}h}({\boldsymbol{x}}_{t+1}+\gamma_{t}{\boldsymbol{u}}_{t})$

15
${\boldsymbol{u}}_{t+1}={\boldsymbol{u}}_{t}+({\boldsymbol{x}}_{t+1}-{\boldsymbol{z}}_{t+1})/\gamma_{t}$

16 $\triangleright$ choose step-size for next iteration, two variants

17 Variant 1

18

19 $\gamma_{t+1}=\gamma_{t}$

20

21 Variant 2 $\triangleright$ only if $h$ is $\beta_{h}$-Lipschitz

22

23
$\delta_{t}=Q_{t}({\boldsymbol{x}}_{t+1},\gamma_{t})-f({\boldsymbol{x}}_{t+1})$

24 Choose any
$\gamma_{t+1}\in[\gamma_{t},\sqrt{\gamma_{t}^{2}+\gamma_{t}\delta_{t}(2\beta_{h})^{-2}}]$

25

return _${\boldsymbol{x}}_{t+1}$ , ${\boldsymbol{u}}_{t+1}$_

Algorithm 1 Adaptive Three Operator Splitting

#### Initial and default values.

The proposed method takes as input 4 parameters, which we briefly discuss,
together with a growing step-size heuristic for Variant 2:

  * •

_Initial guess ${\boldsymbol{z}}_{0}$ and ${\boldsymbol{u}}_{0}$_.
${\boldsymbol{z}}_{0}$ is an initial guess of the primal problem (OPT), while
${\boldsymbol{u}}_{0}$ is an initial guess for a minimizer of a (yet to be
defined) dual function (11). In practice, we initialize both variables to
zero.

  * •

_Initial step-size $\gamma_{0}$_. To estimate a starting value for the step-
size, we start with $\varepsilon=10^{-3}$,
$\widetilde{\boldsymbol{z}}={\boldsymbol{z}}_{0}-\varepsilon\nabla
f({\boldsymbol{z}}_{0})$ and divide $\varepsilon$ by 10 until
$f(\widetilde{\boldsymbol{z}})\leq f({\boldsymbol{z}}_{0})$. Then we solve
$f(\widetilde{\boldsymbol{z}})=Q_{0}(\widetilde{\boldsymbol{z}})$ for
$\gamma_{0}$ and double that estimate, giving

$\gamma_{0}=4(f({\boldsymbol{z}}_{0})-f(\widetilde{\boldsymbol{z}}))\|\nabla
f({\boldsymbol{z}}_{0})\|^{-2}~{}.$ (2)

  * •

The _line search decrease parameter $\tau$_ regulates the factor by which the
step-size is decreased each time the line search condition is unsuccessful.
This is a parameter that is common to all line search methods and can be set
to any value $\tau\in(0,1)$. Following (Malitsky & Pock, 2018) we set it to
$\tau=0.7$.

  * •

_step-size growth._ Variant 2 allows the step-size to grow by an amount that
depends on $\beta_{h}^{-2}$. This quantity can be arbitrarily large (e.g.,
vanishing regularization), and so choosing the largest admissible step-size
might result in too many decrease corrections. This can be avoided e.g. by
limiting its growth to double every $20$ iterations. Line 1 then becomes:

$\gamma_{t+1}=\min\\{\gamma_{t}2^{0.05},\sqrt{\gamma_{t}^{2}+\gamma_{t}\delta_{t}(2\beta_{h})^{-2}}\\}~{}.$
(3)

Upon termination, the algorithm returns two vectors. The first vector is an
approximate solution to (OPT), while the second vector is an approximate
minimizer of a dual objective which we will detail in §3.

#### Special cases and related methods.

We mention two notable special cases of this algorithm. First, for any step-
size $\gamma_{t}\leq 1/L_{f}$, the line search condition will always succeed
by the properties of $L_{f}$-smooth functions and so the step-size in Variant
1 is constant. Defining
${\boldsymbol{y}}_{t}={\boldsymbol{x}}_{t}+\gamma_{t}{\boldsymbol{u}}_{t-1}$,
it is easy to verify that Algorithm 1 (Variant 1) can be written with a
constant step-size $\gamma=\gamma_{t}$ as an iteration of the form

$\displaystyle{\boldsymbol{z}}_{t}=\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{y}}_{t})$ (4)
$\displaystyle{\boldsymbol{x}}_{t+1}\\!=\\!\operatorname{\mathbf{prox}}_{\gamma
g}(2{\boldsymbol{z}}_{t}-{\boldsymbol{y}}_{t}-\gamma\nabla
f({\boldsymbol{z}}_{t}))$
$\displaystyle{\boldsymbol{y}}_{t+1}={\boldsymbol{y}}_{t}-{\boldsymbol{z}}_{t}+{\boldsymbol{x}}_{t+1}~{},$

which is the standard (non-overrelaxed) form of the three operator splitting
(Davis & Yin, 2017, Algorithm 1). The adaptive variant requires two extra
function evaluations $f({\boldsymbol{z}}_{t})$ and $f({\boldsymbol{x}}_{t+1})$
for the line search condition in Line 1, but as we will see in the
experimental section, most often the ability to take larger step outweighs
this extra cost.

Second, for $h=0$, we have from lines 1 and 1 that ${\boldsymbol{u}}_{t}=0$
and in this case (ignoring growing step-size strategies), this algorithm
simplifies to the proximal gradient descent with line search of (Beck &
Teboulle, 2009).

Algorithm 1 can be written equivalently in a way that highlights similarities
and differences with the PDHG method. Using Moreau’s decomposition
$\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{x}})={\boldsymbol{x}}-\gamma\operatorname{\mathbf{prox}}_{\gamma
h^{*}}({\boldsymbol{x}}/\gamma)$ yields the following recurrence

$\displaystyle{\boldsymbol{u}}_{t+1}$
$\displaystyle=\operatorname{\mathbf{prox}}_{h^{*}/\gamma}({\boldsymbol{u}}_{t}+{\boldsymbol{x}}_{t}/\gamma)~{},$
(5) $\displaystyle{\boldsymbol{x}}_{t+2}$
$\displaystyle=\operatorname{\mathbf{prox}}_{\gamma
g}({\boldsymbol{x}}_{t+1}\\!-\gamma(\nabla
f({\boldsymbol{z}}_{t+1})+2{\boldsymbol{u}}_{t+1}\\!-{\boldsymbol{u}}_{t}))\,.$

This form is almost identical to Algorithm 3.2 in (Condat, 2013b), but with a
different step-size and the gradient evaluated at the extrapolated
${\boldsymbol{z}}_{t+1}={\boldsymbol{x}}_{t+1}-\gamma({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{t})$
instead of the previous iterate ${\boldsymbol{x}}_{t+1}$ in PDHG.

###  2.1 Extension to multiple proximal terms

We now consider the problem of minimizing an objective of the form:

$\operatorname*{minimize}_{{\boldsymbol{x}}\in{\mathbb{R}}^{p}}\,\vphantom{\sum_{i}^{n}}\varphi({\boldsymbol{x}})+\textstyle\sum_{j=1}^{k}h_{j}({\boldsymbol{x}})~{},$
(OPT-$k$)

where $\varphi$ is $L_{\varphi}$-smooth and each $h_{j}$ is proximal. The
adaptive three operator splitting can be used to solve problems of this form
by reducing them to a problem of the form (OPT) in an enlarged space. Consider
consider the following problem in ${\mathbb{R}}^{k\times p}$,

$\operatorname*{minimize}_{{\boldsymbol{X}}\in{\mathbb{R}}^{k\times
p}}\underbrace{\varphi(\overline{\boldsymbol{X}})}_{=f({\boldsymbol{X}})}+\underbrace{\textstyle\sum_{j=1}^{k}h_{j}({\boldsymbol{X}}_{j})}_{=h({\boldsymbol{X}})}+\underbrace{\imath\\{{\boldsymbol{X}}_{1}\\!=\\!\cdots\\!=\\!{\boldsymbol{X}}_{k}\\}}_{=g({\boldsymbol{X}})}.$

It is easy to see that this problem shares the same set of solutions as
(OPT-$k$) with the correspondence
${\boldsymbol{x}}=\overline{{\boldsymbol{X}}}$, as the last term forces all
the ${\boldsymbol{X}}_{i}$ terms to be equal. In this formulation the first
term is smooth, the second term is proximal (variables in $h_{i}$ are
separated) and the proximal operator of the last term is given by
$\overline{{\boldsymbol{X}}}\,\mathbf{1}^{T}$. Hence Algorithm 1 can be
applied to this problem. Deriving the complete algorithm is now merely a
matter of replacing $f,g,h$ by its appropriate values in Algorithm 1 and is
specified in Appendix B. The resulting adaptive algorithm seems to be new also
in this extended formulation.

It is also possible to swap the definitions of $g$ and $h$, which results in a
different algorithm that can be seen as an adaptive variant of the Generalized
Forward-Backward splitting of Raguet et al. (2013). However, this formulation
is less convenient for our purpose, since in this case the $h$ term is always
an indicator function and so it would not be possible to apply variant 2 of
our algorithm.

##  3 Analysis

In this section we provide a convergence rate analysis of the proposed method.
We start by a characterization the set of fixed points of the algorithm,
followed by a discussion on the gap function used to measure suboptimality.
Finally, we present convergence rates for two different function classes. All
proofs can be found in Appendix C.

Assumption 1: Regularity. We assume that $f$ is convex and $L_{f}$-smooth in
${\mathbb{R}}^{p}$ and that $g$ and $h$ are proper (i.e., have nonempty
domain), lower semicontinuous (i.e., its sublevel sets are closed) convex
functions. We note that lower semicontinuity is a weak form of continuity that
allows extended-valued functions (such as the indicator function) over a
closed domain.

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 weak and standard assumption that, together with the
regularity assumption, guarantees strong (also known as total) duality
(Bertsekas, 2015, Prop. 5.3.8).

Using the definition of Fenchel conjugate, we can can reformulate (OPT) as a
saddle-point problem:

$\displaystyle\min_{{\boldsymbol{x}}\in{\mathbb{R}}^{d}}f({\boldsymbol{x}})+g({\boldsymbol{x}})+h({\boldsymbol{x}})$
(6)
$\displaystyle\quad=\min_{{\boldsymbol{x}}\in{\mathbb{R}}^{d}}f({\boldsymbol{x}})+g({\boldsymbol{x}})+\max_{{\boldsymbol{u}}\in{\mathbb{R}}^{d}}\left\\{\langle{\boldsymbol{x}},{\boldsymbol{u}}\rangle-h^{*}({\boldsymbol{u}})\right\\}$
(7)
$\displaystyle\quad=\min_{{\boldsymbol{x}}\in{\mathbb{R}}^{d}}\max_{{\boldsymbol{u}}\in{\mathbb{R}}^{d}}\underbrace{f({\boldsymbol{x}})+g({\boldsymbol{x}})+\langle{\boldsymbol{x}},{\boldsymbol{u}}\rangle-h^{*}({\boldsymbol{u}})}_{\stackrel{{\scriptstyle\text{def}}}{{=}}\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}})}~{}.$
(8)

We recall that a saddle point of $\mathcal{L}$ is a pair
$({\boldsymbol{x}}^{\star},{\boldsymbol{u}}^{\star})$ such that the following
is verified for any $({\boldsymbol{x}},{\boldsymbol{u}})$ in the domain
(Hiriart-Urruty & Lemaréchal, 1993, §4.1):

$\mathcal{L}({\boldsymbol{x}}^{\star}\\!,{\boldsymbol{u}})\leq\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}^{\star})~{}.$
(9)

A consequence of strong duality is the equivalence between the saddle points
of $\mathcal{L}$ and the minimizers of the primal and dual objectives. Let $P$
and $D$ denote these primal and dual objectives:

$\displaystyle P({\boldsymbol{x}})$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}f({\boldsymbol{x}})+g({\boldsymbol{x}})+h({\boldsymbol{x}})$
(10) $\displaystyle D({\boldsymbol{u}})$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}(f+g)^{*}(-{\boldsymbol{u}})+h^{*}({\boldsymbol{u}}).$
(11)

Then if $({\boldsymbol{x}}^{\star},{\boldsymbol{u}}^{\star}$) is a saddle
point of $\mathcal{L}$, ${\boldsymbol{x}}^{\star}$ is a minimizer of $P$ and
${\boldsymbol{u}}^{\star}$ is a minimizer of $D$. Likewise, a pair of
minimizers of $P$ and $D$ form a saddle point of $\mathcal{L}$.

###  3.1 Fixed point characterization

A common first step in the analysis of optimization methods is the study of
its set of fixed or stationary points. While this does not necessarily imply
convergence, knowing which elements will be left invariant by the method
improves our understanding and is a stepping stone for further analysis. We
will show that the set of fixed points of the algorithm has a particularly
simple and elegant structure: the Cartesian product of primal and dual
solutions.

For the purpose of analysis it will be useful to express Algorithm 1 as an
iteration of the form,
$({\boldsymbol{z}}_{t+1},{\boldsymbol{u}}_{t+1})=\boldsymbol{T}_{\gamma_{t}}({\boldsymbol{z}}_{t},{\boldsymbol{u}}_{t})$,
where the operator $\boldsymbol{T}_{\gamma}$ is defined as

$\displaystyle\boldsymbol{T}_{\gamma}$
$\displaystyle({\boldsymbol{z}},{\boldsymbol{u}})\stackrel{{\scriptstyle\text{def}}}{{=}}({\boldsymbol{z}}^{+},{\boldsymbol{u}}^{+}),\text{
with }$ (12)
$\displaystyle\begin{cases}{\boldsymbol{z}}^{+}=\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{x}}({\boldsymbol{z}},{\boldsymbol{u}})+\gamma{\boldsymbol{u}})\\\
{\boldsymbol{u}}^{+}={\boldsymbol{u}}+({\boldsymbol{x}}({\boldsymbol{z}},{\boldsymbol{u}})-{\boldsymbol{z}}^{+})/\gamma\\\
{\boldsymbol{x}}({\boldsymbol{z}},{\boldsymbol{u}})=\operatorname{\mathbf{prox}}_{\gamma
g}({\boldsymbol{z}}-\gamma({\boldsymbol{u}}+\nabla
f({\boldsymbol{z}})))~{}.\end{cases}$

The following theorem characterizes the set of fixed points of this operator,
denoted $\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})$.

######  Theorem 1.

Let $\mathcal{P}^{\star}$ denote the set of minimizers of the primal objective
$P$ and $\mathcal{D}^{\star}$ the set of minimizers of the dual objective $D$
. Then the fixed points of ${\boldsymbol{T}}_{\gamma}$ are given by

$\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})=\mathcal{P}^{\star}\times\mathcal{D}^{\star}~{}.$
(13)

###  3.2 Gap function

The progress of optimization methods is commonly measured in terms of a gap or
merit function that is zero at optimum and nonzero otherwise. An appropriate
gap function for many first-order methods is the suboptimality of the
objective function. However, this is not an appropriate suboptimality measure
for this algorithm, as the objective function might be $+\infty$ at an
iterate, for example when the two proximal terms are an indicator function.

Davis & Yin (2015) avoid the issue by either evaluating $h$ and $g$ at
different iterates (Davis & Yin, 2015, Corollary D.5.1) or assuming Lipschitz
continuity of one of the proximal terms (Davis & Yin, 2015, Corollary D.5.2).

In this work we take an alternative approach, and instead use the following
_saddle point suboptimality_ criterion to measure the progress of our
algorithm:

${\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1})}~{}.$
(14)

From the definition of saddle point in Eq. (9), this criterion is non-positive
for all $({\boldsymbol{x}},{\boldsymbol{u}})$ if and only if
$({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t+1})$ is a saddle point, and is so
an appropriate suboptimality criterion. Furthermore, contrary to the primal
objective function, this is defined for all iterates without further
assumptions. Finally, we mention that this criteria has been previously used
in the analysis of primal-dual methods, see e.g., Chambolle & Pock (2016,
2015) and Gidel et al. (2017) for a discussion of saddle point gap functions.

This suboptimality criteria can also be related to the primal and dual gap, as
minimizing (14) over ${\boldsymbol{x}}$ and maximizing over ${\boldsymbol{u}}$
one recovers the primal-dual gap
$P({\boldsymbol{x}}_{t})-D({\boldsymbol{u}}_{t})$ by definition of Fenchel
conjugate.

###  3.3 Sublinear convergence

The following theorem gives a sublinear convergence rate for Algorithm 1. This
convergence will be given in terms of the weighted ergodic (i.e., averaged)
sequence. Denoting by $s_{t}$ the sum of all step-sizes up to iteration $t$,
i.e.,
$s_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{i=0}^{t-1}\gamma_{t}$, the
ergodic iterates $\overline{{\boldsymbol{x}}}_{t}$ and
$\overline{{\boldsymbol{u}}}_{t}$ are defined as

$\overline{{\boldsymbol{x}}}_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}\Big{(}\sum_{i=0}^{t-1}\gamma_{i}{\boldsymbol{x}}_{i+1}\Big{)}/{s_{t}}\,,\;\;\overline{{\boldsymbol{u}}}_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}\Big{(}{\sum_{i=0}^{t-1}\gamma_{i}{\boldsymbol{u}}_{i+1}}\Big{)}/{s_{t}}\,.$
(15)

While results in this subsection will be stated in terms of this ergodic
sequence, in practice the last iterate gives most often a better empirical
convergence, see e.g., (Chambolle & Pock, 2015, §7.2.1) for a discussion of
this phenomenon. For a more theoretically-sound algorithm, one can compare the
objective at the ergodic and last iterate, and return the one with smallest
objective.

######  Theorem 2 (sublinear convergence rate).

For every ${t\geq 0}$ and any $({\boldsymbol{x}},{\boldsymbol{u}})$ in the
domain of $\mathcal{L}$ we have the following convergence rate for Algorithm 1
(both variants):

$\mathcal{L}(\overline{\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},\overline{\boldsymbol{u}}_{t+1})\leq\frac{\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+{\gamma_{0}^{2}\|{\boldsymbol{u}}_{0}-{\boldsymbol{u}}\|^{2}}}{2s_{t}}\,.$

#### Convergence in terms of function value suboptimality.

The previous result gives an $\mathcal{O}(1/t)$ convergence rate for arbitrary
convex functions in terms of the saddle point suboptimality. As we have
discussed previously, it is not possible to obtain similar rates in terms of
the function suboptimality without further assumptions. We will now show that
it is sufficient to assume Lipschitz continuity on $h$ to derive from the
previous theorem a convergence rate in terms of the primal function
suboptimality.

The following Corollary can be obtained by optimizing with respect to
${\boldsymbol{u}}$ the bound in the previous theorem and using the Lipschitz
continuity to bound $\|{\boldsymbol{u}}_{0}-{\boldsymbol{u}}\|^{2}$. This
gives an $\mathcal{O}(1/t)$ convergence rate for the primal function
suboptimality, roughly matching that of Davis & Yin (2015, Corollary D.5.2)
for the non adaptive variant:

######  Corollary 1.

Let $h$ be $\beta_{h}$-Lipschitz. Then, we have the following rate for the
weighted ergodic iterate on the objective
$P({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}f({\boldsymbol{x}})+g({\boldsymbol{x}})+h({\boldsymbol{x}})$:

$P(\overline{\boldsymbol{x}}_{t+1})-P({\boldsymbol{x}}^{\star})\leq\frac{\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}\\!+2\gamma_{0}^{2}(\|{\boldsymbol{u}}_{0}\|^{2}\\!+\beta_{h}^{2})}{2s_{t}}\,.$

###  3.4 Linear convergence

In this subsection we assume that $f$ is $\mu_{f}$-strongly convex and $h$ is
$L_{h}$-smooth (with $0<\mu_{f},0<L_{h}<+\infty$). We denote by
${\boldsymbol{x}}^{\star}$ the minimizer of the primal loss (unique by strong
convexity of $P$) and by ${\boldsymbol{u}}^{\star}$ the minimizer of the dual
loss (also unique by strong convexity of $D$, consequence of the
$L_{h}$-smoothness of $h$).

The convergence rates will be given in terms of the following quantities

$\displaystyle\rho\stackrel{{\scriptstyle\text{def}}}{{=}}\mu_{f}\min\\{\gamma_{0},\tau/L_{f}\\}~{},~{}\sigma\stackrel{{\scriptstyle\text{def}}}{{=}}1/(1+\gamma_{0}L_{h})$
(16)
$\displaystyle\qquad\qquad\quad\xi\stackrel{{\scriptstyle\text{def}}}{{=}}{\mu_{f}}/{(\mu_{f}+L_{h})}~{}.$

All these belong to the interval $(0,1)$. Assuming $\gamma_{0}\geq\tau/L_{f}$,
then $\rho$ is the inverse of $f$’s condition number, while $\sigma$ depends
the smoothness of $h$. $\xi$ is only used by variant 2 and is a less tight
bound that $\sigma$ that depends on both the strong convexity of $f$ and
smoothness of $h$. Note that by strong convexity, $\gamma_{0}<1/\mu_{f}$ as
otherwise the sufficient decrease condition would not succeed and so
$\sigma\geq\xi$.

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

Let ${\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t+1}$ be the iterates produced
by Algorithm 1 after $t$ iterations. Then we have the following linear
convergence for Variant 1 (V1) and Variant 2 (V2):

$\displaystyle\text{V1}\mathrel{\mathop{\ordinarycolon}}\|{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}\leq\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}D_{0}$
(17)
$\displaystyle\text{V2}\mathrel{\mathop{\ordinarycolon}}\|{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}E_{0}~{},$
(18)

with
$D_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}6\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\frac{6}{1-\sigma}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}}^{\star})\|^{2}$
and
${E_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}6\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\frac{6}{1-\xi}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}}^{\star})\|^{2}}$.

#### Discussion.

For $\gamma_{t}={1}/{L_{f}}$, the sufficient decrease condition is always
verified and the algorithm can be run with $\tau=1$. In this case, Variant 1
of Algorithm 1 defaults to TOS, and we can compare the obtained rates with
those in (Davis & Yin, 2015).

While the sublinear convergence rate obtained in Corollary 1 roughly matches
the rate obtained in (Davis & Yin (2015, Corollary D.5.2, see our Appendix
C.5)), the linear convergence rates are instead significantly different. The
linear convergence rate obtained in (Davis & Yin, 2015, Theorem D.6.6), after
optimizing for all parameters, yields a rate of $\rho\sigma^{2}$, which is
_strictly worse_ than the $\min\\{\rho,\sigma\\}$ rate that we obtained. This
difference can be quite large, e.g., for $\rho=\sigma$ this becomes $\rho$
versus $\rho^{3}$.

Finally, we note that the number of evaluations of the sufficient decrease
condition can be bounded as in (Nesterov et al., 2013, Lemma 4).

##  4 Applications

###  4.1 Learning with Multiple Penalties

In this subsection we discuss how some penalties with costly to compute
proximal operator can be decomposed as a sum of proximal terms and so fall
within the current framework. The exact expression of the proximal operators
is given in Appendix D.

#### Group lasso with overlap.

Jacob et al. (2009) generalized group $\ell_{1}$ norm by allowing each
variable to belong to more than one group, thereby introducing overlaps among
groups and allowing for more complex prior knowledge on the structure. For a
set of subindices $\mathcal{G}$ which we will call groups, this penalty is
defined as
$\|{\boldsymbol{x}}\|_{\mathcal{G}}=\sum_{G\in\mathcal{G}}\|[{\boldsymbol{x}}]_{G}\|_{2}$.
If each coefficient is at most in $s$ groups, then $\mathcal{G}$ can be
decomposed as $\mathcal{G}=\mathcal{G}_{1}\cup\ldots\cup\mathcal{G}_{s}$,
where the $\mathcal{G}_{i}$ are disjoint. This allows to express the group
lasso with overlap as a sum of $s$ non-overlapping group lasso penalties, for
which the proximal operator has a closed form expression.

#### Multidimensional total variation.

For the task of image restoration and denoising it is common to consider a
regularization term in the form of a total variation regularizer. For an image
${\boldsymbol{x}}$ of size $p\times q$, the 2-dimensional total variation norm
$\|{\boldsymbol{X}}\|_{\text{TV}}$ is defined as

$\underbrace{{\sum_{i=1}^{p}}{\sum_{j=1}^{q-1}}|{\boldsymbol{X}}_{i,j+1}-{\boldsymbol{X}}_{i,j}|}_{=g({\boldsymbol{X}})}+\underbrace{\sum_{j=1}^{q}{\sum_{ji=1}^{p-1}}|{\boldsymbol{X}}_{i+1,j}-{\boldsymbol{X}}_{i,j}|}_{=h({\boldsymbol{X}})}\,.$

From here we recognize that $g$ and $h$ are fused lasso (also known as
1D-total variation) penalties acting on the columns and rows of
${\boldsymbol{X}}$ respectively. Efficient methods to evaluate the proximal
operator of the fused lasso penalty have been developed by Condat (2013a);
Johnson (2013).

#### Isotonic and nearly isotonic penalties.

In some applications there exists a natural ordering between variables:
${\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}\leq\cdots\leq{\boldsymbol{x}}_{p}$.
This can be enforced through constraints, and the projection onto these is
known as isotonic regression (Best & Chakravarti, 1990). The indicator
function over the set of constraints can also be split into a sum of two
proximal terms (see Appendix D.2) as

$\displaystyle\imath\\{{\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}\leq{\boldsymbol{x}}_{3}\leq{\boldsymbol{x}}_{4}\leq\cdots\\}$
(19)
$\displaystyle=\underbrace{\imath\\{{\boldsymbol{x}}_{1}\\!\leq\\!{\boldsymbol{x}}_{2};{\boldsymbol{x}}_{3}\leq{\boldsymbol{x}}_{4};\cdots\\}}_{=g({\boldsymbol{x}})}+\underbrace{\imath\\{{\boldsymbol{x}}_{2}\leq{\boldsymbol{x}}_{3};{\boldsymbol{x}}_{4}\leq{\boldsymbol{x}}_{5};\cdots\\}}_{=h({\boldsymbol{x}})}\,.$

In cases in which the variables are only “mostly” non-decreasing, the
constraint can be relaxed via a nearly-isotonic penalty (Tibshirani et al.,
2011) of the form
$\sum_{i=1}^{p-1}\max\\{{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{i+1},0\\}$, in
which only the non-increasing coefficients are penalized. This penalty can be
split the same way as the isotonic constraints above.

####  $\ell_{1}$ trend filtering.

This penalty is defined by the absolute value of the second order differences
and promotes piecewise-linear coefficients (Kim et al., 2009). It is defined
as
$\|{\boldsymbol{x}}\|_{\text{TF}}\stackrel{{\scriptstyle\text{def}}}{{=}}\textstyle\sum_{i=1}^{p-2}|{\boldsymbol{x}}_{i}-2{\boldsymbol{x}}_{i+1}+{\boldsymbol{x}}_{i+2}|$.
We can split this sum into 3 proximal terms such that the resulting terms: the
$j$-th term contains the factors for which $i$ is congruent to 3 modulo $j$.

#### Constraints over doubly stochastic matrices.

Optimization problems with constraints on the set of doubly stochastic
matrices appear in many convex relaxations of combinatorial problems such as
seriation (Fogel et al., 2013), quadratic assignment (Lawler, 1963) and graph
matching (Conte et al., 2004; Aflalo et al., 2015). The set of double
stochastic matrices is composed of square matrices with nonnegative entries,
each of whose rows and columns sum to $1$, i.e.,
$\\{{\boldsymbol{X}}^{T}\boldsymbol{1}=\boldsymbol{1},{\boldsymbol{X}}\boldsymbol{1}=\boldsymbol{1},{\boldsymbol{X}}\geq\boldsymbol{0}\\}$.
The indicator function over this set can be split as

$\underbrace{\imath\\{{\boldsymbol{X}}^{T}\boldsymbol{1}=\boldsymbol{1},{\boldsymbol{X}}\boldsymbol{1}=\boldsymbol{1}\\}}_{=g({\boldsymbol{X}})}+\underbrace{\imath\\{{\boldsymbol{X}}\geq\boldsymbol{0}\\}}_{=h({\boldsymbol{X}})}\quad,$
(20)

and the projection onto both sets is available in closed form (Lu et al.,
2016, §4.3).

Figure 1: Comparison of different proximal splitting methods. The _top row_
gives result for two real datasets with an overlapping group lasso penalty.
The _second and third row_ show results on synthetic datasets for 4 different
penalties: overlapping group lasso (E, F), 2-dimensional total variation (G,
H), trace norm + $\ell_{1}$ (I, J) and nearly isotonic (K, L). The Adaptive
TOS (Variant 2, i.e., with growing step-sizes) is the best performing method
on 10 out of 12 experiments, and roughly equivalent to the best performing
method in the other 2 cases (G, J).

#### Dispersive sparsity.

In some applications it is desirable to encourage dispersion of the sparse
coefficients. This happens for example in the modeling of neural spiking, as
the spikes are assumed to be spaced across time (Hegde et al., 2009). El
Halabi & Cevher (2015) showed that this behavior can be promoted by
considering a penalty of the form
$\|{\boldsymbol{x}}\|_{1}+\imath\\{\boldsymbol{B}|{\boldsymbol{x}}|\leq c\\}$
for a matrix $\boldsymbol{B}$ and some predefined constant $c$, where
$|{\boldsymbol{x}}|$ denotes the component-wise absolute value. This penalty
can be split into three proximal terms by the introduction of a dummy variable
${\boldsymbol{z}}$, resulting in
$\|{\boldsymbol{x}}\|_{1}+\imath\\{\boldsymbol{B}{\boldsymbol{z}}\leq\boldsymbol{c}\\}+\imath\\{{\boldsymbol{z}}=|{\boldsymbol{x}}|\\}$.

Combination by addition. A popular method to promote the joint behavior of
different penalties is by adding them. This has been used to successfully
learn models with sparse and nonnegative coefficients (Yuan & Lin, 2007),
sparse and low rank matrices (Richard et al., 2012), sparse and piecewise
constant (Gramfort et al., 2013), to name a few.

###  4.2 Benchmarks

In this subsection we provide an empirical evaluation of the proposed method.
Due to space constraints we only give here a high level overview, deferring
details as well as an extended set of experiments to Appendix E. We consider
the following methods:

  * •

The proposed Adaptive TOS method (Algorithm 1), in its both variants.

  * •

The TOS method of Davis & Yin (2015), with step-sizes $1/L_{f}$ and
$1.99/L_{f}$ (the method is convergent for step-sizes $<2/L_{f}$).

  * •

The PDHG or Condat-Vũ algorithm (Condat, 2013b; Vũ, 2013), with step-sizes
$\tau$ and $\beta/\tau$, where $\beta$ was chosen as the one giving the best
overall performance over the grid $\beta=0.9,0.5,0.1$ (giving it a slight
advantage).

  * •

The adaptive PDHG of Malitsky & Pock (2018), with step-size hyperparameter
$\beta$ chosen by the same technique as for PDHG.

  * •

The averaged operator line search method of Giselsson et al. (2016) combined
with TOS, named TOS-AOLS.

We compared these methods on 4 different problems and show the results in
Figure 1. In the first row we show the benchmarks on a logistic regression
problem with overlapping group lasso penalty that we apply to two text
datasets (RCV1 and real-sim). Subfigures A and C were run with the
regularization parameter chosen to give $50\%$ of sparsity, while B, E are run
with higher levels of sparsity, chosen to give $5\%$ of sparsity.

In the second and third row we considered a battery of inverse problems with
different penalties on synthetic datasets. These consists of a least squares
(G, H, I, J) or logistic regression (rest) smooth term and 4 different
penalties specified in the title of each plot (overlapping group lasso, total
variation, trace norm $\ell_{1}$ and nearly isotonic, see Appendix E for a
precise formulation). For each problem, we show 2 different benchmarks,
corresponding to the low and high regularization regimes (denoted low reg and
high reg). We comment on a few trends from Fig. 1:

  * •

Best performing method. On 10 out of 12 experiments, the adaptive TOS
algorithm (Variant 2) is the best performing method, and in the other cases
(E, H) its performance is roughly the same as that of the best performing
method. In contrast, on 3 instances (A, I, K) it is an order of magnitude
faster than the next method.

  * •

Low vs high regularization regime. The advantage of the adaptive method is
highly correlated with the amount of regularization: in the low regularization
regime, on 3 out of 6 the adaptive TOS is an order of magnitude faster then
the fixed step-size method, while in the high regularization regime the
difference shrinks and in the same problems is never more than a factor 2.6.

  * •

Uniform curvature. The problems (G, H, I, J) in Fig. 1 use as smooth term a
quadratic loss (i.e., constant Hessian), while the other methods use a
logistic loss (non-constant Hessian). This suggests that the use of the
adaptive step-size strategy (and in particular Variant 2 with its growing
step-size) is more beneficial for smooth terms with non-uniform curvature.

##  5 Conclusion and Future Work

We have presented and analyzed an adaptive step-size method to solve
optimization problems consisting in a sum of a smooth term accessed through
its gradient and two or more potentially non-smooth terms accessed through
their proximal operator. The method does not rely on any step-size
hyperparameter (except for an initial estimate) and extensive empirical
evaluation has showed computational gains on a variety of problems. We mention
two possible extensions of this work.

First, existing convergence results fail to fully explain their surprisingly
good empirical convergence. To the best of our knowledge, no work so far has
derived linear convergence rates in absence of strong convexity and smoothness
of one of the proximal terms for these methods (as is however empirically
observed, see e.g. Figure 1).

Second, it is an open question whether this or other adaptive step-size
methods can be accelerated, as is the case of proximal gradient descent, which
admits the adaptive FISTA variant (Beck & Teboulle, 2009).

## Acknowledgements

FP would like to thank Wotao Yin, Walaa Moursi and Elvis Dohmatob for fruitful
discussions, Michael Eickenberg for feedback on the manuscript and Patrick R.
Johnstone for reporting a typo in an earlier version.

FP is funded through the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodorowska-Curie grant agreement 748900. Computing
time on was donated by Amazon through the program “AWS Cloud Credits for
Research”.

GG is funded by the Canada Excellence Research Chair in “Data Science for
Realtime Decision-making” and by the NSERC Discovery Grant RGPIN-2017-06936.

## References

  * Aflalo et al. (2015) Aflalo, Y., Bronstein, A., and Kimmel, R.  On convex relaxation of graph isomorphism.  _Proceedings of the National Academy of Sciences_ , 2015.
  * Agarwal et al. (2010) Agarwal, A., Negahban, S., and Wainwright, M. J.  Fast global convergence rates of gradient methods for high-dimensional statistical recovery.  In _Advances in Neural Information Processing Systems_ , 2010.
  * Barbero & Sra (2014) Barbero, Á. and Sra, S.  Modular proximal optimization for multidimensional total-variation regularization.  _preprint arXiv:1411.0589_ , 2014.
  * Bauschke & Combettes (2017) Bauschke, H. H. and Combettes, P. L.  _Convex analysis and monotone operator theory in Hilbert spaces_.  Springer Science & Business Media, 2017.
  * Bauschke et al. (2012) Bauschke, H. H., Boţ, R. I., Hare, W. L., and Moursi, W. M.  Attouch–Théra duality revisited: paramonotonicity and operator splitting.  _Journal of Approximation Theory_ , 2012.
  * Beck & Teboulle (2009) Beck, A. and Teboulle, M.  Gradient-based algorithms with applications to signal recovery.  _Convex Optimization in Signal Processing and Communications_ , 2009\.
  * Bertsekas (2015) Bertsekas, D. P.  _Convex optimization algorithms_.  Athena Scientific Belmont, 2015.
  * Best & Chakravarti (1990) Best, M. J. and Chakravarti, N.  Active set algorithms for isotonic regression; A unifying framework.  _Mathematical Programming_ , 1990.
  * Boyd et al. (2011) Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  _Foundations and Trends in Machine Learning_ , 2011.
  * Cai et al. (2010) Cai, J.-F., Candès, E. J., and Shen, Z.  A singular value thresholding algorithm for matrix completion.  _SIAM Journal on Optimization_ , 2010.
  * Candès et al. (2006) Candès, E. J., Romberg, J., and Tao, T.  Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information.  _IEEE Transactions on information theory_ , 2006.
  * Chambolle & Pock (2015) Chambolle, A. and Pock, T.  On the ergodic convergence rates of a first-order primal–dual algorithm.  _Mathematical Programming_ , 2015.
  * Chambolle & Pock (2016) Chambolle, A. and Pock, T.  An introduction to continuous optimization for imaging.  _Acta Numerica_ , 2016.
  * Combettes & Pesquet (2011) Combettes, P. L. and Pesquet, J.-C.  Proximal splitting methods in signal processing.  In _Fixed-point algorithms for inverse problems in science and engineering_. Springer, 2011.
  * Condat (2013a) Condat, L.  A direct algorithm for 1D total variation denoising.  _IEEE Signal Processing Letters_ , 2013a.
  * Condat (2013b) Condat, L.  A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms.  _Journal of Optimization Theory and Applications_ , 2013b.
  * Conte et al. (2004) Conte, D., Foggia, P., Sansone, C., and Vento, M.  Thirty years of graph matching in pattern recognition.  _International journal of pattern recognition and artificial intelligence_ , 2004.
  * Davis & Yin (2015) Davis, D. and Yin, W.  A three-operator splitting scheme and its optimization applications.  _preprint arXiv:1504.01032v1_ , 2015.
  * Davis & Yin (2017) Davis, D. and Yin, W.  A three-operator splitting scheme and its optimization applications.  _Set-Valued and Variational Analysis_ , 2017.
  * El Halabi & Cevher (2015) El Halabi, M. and Cevher, V.  A totally unimodular view of structured sparsity.  In _Proceedings of the 18th International Conference on Artificial Intelligence and Statistics (AISTATS)_ , 2015.
  * Fogel et al. (2013) Fogel, F., Jenatton, R., Bach, F., and d’Aspremont, A.  Convex relaxations for permutation problems.  In _Advances in Neural Information Processing Systems_ , 2013.
  * Gabay & Mercier (1976) Gabay, D. and Mercier, B.  A dual algorithm for the solution of nonlinear variational problems via finite element approximation.  _Computers & Mathematics with Applications_, 1976.
  * Gidel et al. (2017) Gidel, G., Jebara, T., and Lacoste-Julien, S.  Frank-Wolfe Algorithms for Saddle Point Problems.  In _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics_ , 2017.
  * Giselsson et al. (2016) Giselsson, P., Fält, M., and Boyd, S.  Line search for averaged operator iteration.  In _IEEE 55th Conference on Decision and Control (CDC)_. IEEE, 2016\.
  * Glowinski & Marroco (1975) Glowinski, R. and Marroco, A.  Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de Dirichlet non linéaires.  _Revue française d’automatique, informatique, recherche opérationnelle_ , 1975.
  * Gramfort et al. (2013) Gramfort, A., Thirion, B., and Varoquaux, G.  Identifying predictive regions from fMRI with TV-L1 prior.  In _International Workshop on Pattern Recognition in Neuroimaging_. IEEE, 2013.
  * Hegde et al. (2009) Hegde, C., Duarte, M. F., and Cevher, V.  Compressive sensing recovery of spike trains using a structured sparsity model.  In _SPARS’09-Signal Processing with Adaptive Sparse Structured Representations_ , 2009.
  * Hiriart-Urruty & Lemaréchal (1993) Hiriart-Urruty, J.-B. and Lemaréchal, C.  _Convex analysis and minimization algorithms I: Fundamentals_.  Springer science & business media, 1993.
  * Iusem (1998) Iusem, A. N.  On Some Properties of Generalized Proximal Point Methods for Variational Inequalities.  _Journal of Optimization Theory and Applications_ , 1998.
  * Jacob et al. (2009) Jacob, L., Obozinski, G., and Vert, J.-P.  Group lasso with overlap and graph lasso.  In _Proceedings of the 26th annual international conference on machine learning_. ACM, 2009.
  * Johnson (2013) Johnson, N.  A dynamic programming algorithm for the fused lasso and $L_{0}$-segmentation.  _Journal of Computational and Graphical Statistics_ , 2013.
  * Kim et al. (2009) Kim, S.-J., Koh, K., Boyd, S., and Gorinevsky, D.  $\ell_{1}$ trend filtering.  _SIAM review_ , 2009.
  * Lawler (1963) Lawler, E. L.  The quadratic assignment problem.  _Management science_ , 1963.
  * Lewis et al. (2004) Lewis, D. D., Yang, Y., Rose, T. G., and Li, F.  RCV1: A new benchmark collection for text categorization research.  _Journal of machine learning research_ , 2004.
  * Lions & Mercier (1979) Lions, P.-L. and Mercier, B.  Splitting algorithms for the sum of two nonlinear operators.  _SIAM Journal on Numerical Analysis_ , 1979.
  * Lu et al. (2016) Lu, Y., Huang, K., and Liu, C.-L.  A fast projected fixed-point algorithm for large graph matching.  _Pattern Recognition_ , 2016.
  * Malitsky & Pock (2018) Malitsky, Y. and Pock, T.  A first-order primal-dual algorithm with linesearch.  _SIAM Journal on Optimization_ , 2018.
  * Nesterov et al. (2013) Nesterov, Y. et al.  Gradient methods for minimizing composite objective function.  _Mathematical Programming_ , 2013.
  * Nocedal & Wright (2006) Nocedal, J. and Wright, S.  _Numerical optimization_.  Springer Science & Business Media, 2006.
  * Parikh & Boyd (2013) Parikh, N. and Boyd, S.  Proximal algorithms.  _Foundations and Trends in optimization_ , 2013.
  * Pedregosa (2018) Pedregosa, F.  C-OPT: composite optimization in Python.  2018\.  doi: 10.5281/zenodo.1283339.  URL http://openopt.github.io/copt/.
  * Raguet et al. (2013) Raguet, H., Fadili, J., and Peyré, G.  A generalized forward-backward splitting.  _SIAM Journal on Imaging Sciences_ , 2013.
  * Richard et al. (2012) Richard, E., Savalle, P.-a., and Vayatis, N.  Estimation of Simultaneously Sparse and Low Rank Matrices.  In _Proceedings of the 29th International Conference on Machine Learning_ , 2012.
  * Richtárik & Takáč (2014) Richtárik, P. and Takáč, M.  Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function.  _Mathematical Programming_ , 2014.
  * Rockafellar (1997) Rockafellar, R. T.  Convex analysis, 1997.
  * Rockafellar & Wets (1998) Rockafellar, R. T. and Wets, R. J.-B.  _Variational analysis_.  Springer Science & Business Media, 1998.
  * Rudin et al. (1992) Rudin, L. I., Osher, S., and Fatemi, E.  Nonlinear total variation based noise removal algorithms.  _Physica D: Nonlinear Phenomena_ , 1992.
  * Tibshirani et al. (2011) Tibshirani, R. J., Hoefling, H., and Tibshirani, R.  Nearly-isotonic regression.  _Technometrics_ , 2011.
  * Vũ (2013) Vũ, B. C.  A splitting algorithm for dual monotone inclusions involving cocoercive operators.  _Advances in Computational Mathematics_ , 2013.
  * Yan (2018) Yan, M.  A New Primal–Dual Algorithm for Minimizing the Sum of Three Functions with a Linear Operator.  _Journal of Scientific Computing_ , 2018.
  * Yuan & Lin (2006) Yuan, M. and Lin, Y.  Model selection and estimation in regression with grouped variables.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 2006.
  * Yuan & Lin (2007) Yuan, M. and Lin, Y.  On the non-negative garrotte estimator.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 2007.

Adaptive Three Operator Splitting

Supplementary material

##  Appendix Appendix A Basic definitions and properties

######  Definition 1 (proper function).

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

######  Definition 2 (Fenchel conjugate).

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

$f^{*}({\boldsymbol{x}}^{\star})=\sup_{{\boldsymbol{x}}\in\mathcal{X}}\,\langle{\boldsymbol{x}}^{\star},{\boldsymbol{x}}\rangle-f({\boldsymbol{x}})~{}.$
(21)

######  Definition 3 (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 4 (strong convexity).

A function $f$ is said $\mu$-strongly convex if it verifies the following
inequality for all ${\boldsymbol{x}},{\boldsymbol{y}}$ in the domain and any
${\boldsymbol{u}}\in\partial f({\boldsymbol{x}})$

$f({\boldsymbol{x}})\leq
f({\boldsymbol{y}})+\langle{\boldsymbol{u}},{\boldsymbol{x}}-{\boldsymbol{y}}\rangle-\frac{\mu}{2}\|{\boldsymbol{x}}-{\boldsymbol{y}}\|^{2}$
(22)

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

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

$\operatorname{relint}(C)=\\{{\boldsymbol{x}}\in
C\mathrel{\mathop{\ordinarycolon}}\forall{{\boldsymbol{y}}\in
C}\;\exists{\lambda>1}\mathrel{\mathop{\ordinarycolon}}\lambda{\boldsymbol{x}}+(1-\lambda){\boldsymbol{y}}\in
C\\}$ (23)

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

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

${\boldsymbol{z}}=\operatorname{\mathbf{prox}}_{\gamma
g}({\boldsymbol{x}})\iff\frac{1}{\gamma}({\boldsymbol{x}}-{\boldsymbol{z}})\in\partial
g({\boldsymbol{z}})$ (24)

###### Proof.

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

$\displaystyle{\boldsymbol{z}}\in\operatorname*{\mathbf{arg\,min}}_{{\boldsymbol{z}}^{\prime}}g({\boldsymbol{z}}^{\prime})+\frac{1}{2\gamma}\|{\boldsymbol{z}}^{\prime}-{\boldsymbol{x}}\|^{2}$
(25) $\displaystyle\iff 0\in\partial
g({\boldsymbol{z}})+\frac{1}{\gamma}({\boldsymbol{z}}-{\boldsymbol{x}})$ (26)
$\displaystyle\iff\frac{1}{\gamma}({\boldsymbol{x}}-{\boldsymbol{z}})\in\partial
g({\boldsymbol{z}})$ (27)

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}})~{}.$ (28)

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

###### Proof.

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

######  Lemma 3 (Improved Young’s inequality).

For all $\boldsymbol{a},\boldsymbol{b},\boldsymbol{c}$ we have,

$\|\boldsymbol{a}+\boldsymbol{b}+\boldsymbol{c}\|^{2}\leq
3(\|\boldsymbol{a}\|^{2}+\|\boldsymbol{b}\|^{2}+\|\boldsymbol{c}\|^{2})$ (29)

###### Proof.

We develop the squared norm to get,

$\displaystyle\|\boldsymbol{a}+\boldsymbol{b}+\boldsymbol{c}\|^{2}$
$\displaystyle=\|\boldsymbol{a}\|^{2}+2\langle\boldsymbol{a},\boldsymbol{b}\rangle+\|\boldsymbol{b}\|^{2}+2\langle\boldsymbol{b},\boldsymbol{c}\rangle+\|\boldsymbol{c}\|^{2}+2\langle\boldsymbol{a},\boldsymbol{c}\rangle$
(30)
$\displaystyle\leq\|\boldsymbol{a}\|^{2}+\|\boldsymbol{a}\|^{2}+\|\boldsymbol{b}\|^{2}+\|\boldsymbol{b}\|^{2}+\|\boldsymbol{b}\|^{2}+\|\boldsymbol{c}\|^{2}+\|\boldsymbol{c}\|^{2}+\|\boldsymbol{a}\|^{2}+\|\boldsymbol{c}\|^{2}$
(31) $\displaystyle\qquad(\text{Young's inequality on
}2\langle\boldsymbol{a},\boldsymbol{b}\rangle,2\langle\boldsymbol{b},\boldsymbol{c}\rangle\text{
and }2\langle\boldsymbol{a},\boldsymbol{c}\rangle)$
$\displaystyle=3(\|\boldsymbol{a}\|^{2}+\|\boldsymbol{b}\|^{2}+\|\boldsymbol{c}\|^{2})$
(32)

∎

##  Appendix Appendix B Adaptive (k+1) operator splitting

For completeness we write below the full adaptive $(k+1)$ operator splitting,
which is merely Algorithm 1 applied to the reformulation of (OPT-$k$) in §2.1.
It hence solves an optimization problem of the form

$\operatorname*{minimize}_{{\boldsymbol{x}}\in{\mathbb{R}}^{p}}\,\vphantom{\sum_{i}^{n}}\varphi({\boldsymbol{x}})+\textstyle\sum_{j=1}^{k}h_{j}({\boldsymbol{x}})~{},$
(33)

where $\varphi$ is $L_{f}$-smooth and each $h_{j}$ is proximal.

In this case we name the quadratic used in the sufficient decrease condition
$\widetilde{Q}_{t}$, and is defined as

$\displaystyle\widetilde{Q}_{t}({\boldsymbol{x}})\stackrel{{\scriptstyle\text{def}}}{{=}}\varphi(\overline{\boldsymbol{Z}}_{t})+\langle\nabla\varphi(\overline{\boldsymbol{Z}}_{t}),{\boldsymbol{x}}-\overline{\boldsymbol{Z}}_{t}\rangle+\frac{1}{2\gamma_{t}}\|{\boldsymbol{x}}\mathbf{1}_{k}^{T}-{\boldsymbol{Z}}_{t}\|_{F}^{2}~{},$
(34)

where $\|\cdot\|_{F}$ denotes the Frobenius norm and $\mathbf{1}_{k}$ is the
$k$-dimensional vector of ones.

Input: ${\boldsymbol{z}}_{0}\in{\mathbb{R}}^{p}$,
${\boldsymbol{u}}_{0}\in{\mathbb{R}}^{p}$, $\gamma_{0}>0,\tau\in(0,1)$

1 ${\boldsymbol{U}}_{0}={\boldsymbol{u}}_{0}\boldsymbol{1}_{k}^{T}$

2${\boldsymbol{Z}}_{0}={\boldsymbol{z}}_{0}\boldsymbol{1}_{k}^{T}$

3for  _$t=1,\ldots,T$ _ do

4

5 repeat $\triangleright$ step-size search loop

6

7 ${\boldsymbol{r}}_{t}=\frac{1}{k}\nabla\varphi({\boldsymbol{z}}_{t})$

8
${\boldsymbol{x}}_{t+1}=\overline{{\boldsymbol{Z}}}_{t}-\gamma_{t}\overline{\boldsymbol{U}}_{t}-\gamma_{t}{\boldsymbol{r}}_{t}$

9
$\delta_{t}=\widetilde{Q}_{t}({\boldsymbol{x}}_{t+1})-\varphi({\boldsymbol{x}}_{t+1})$

10 if  _$\delta_{t}\geq 0\,$_ then

11 break $\triangleright$ sufficient decrease verified

12 else

13 $\gamma_{t}=\tau\gamma_{t}$ $\triangleright$ decrease step-size

14

15

16 for  _$j=1,\ldots,k\,$_ do

17

18 ${\boldsymbol{Z}}_{j,t+1}=\operatorname{\mathbf{prox}}_{\gamma
h_{j}}({{\boldsymbol{x}}}_{t+1}+\gamma_{t}{{\boldsymbol{U}}}_{j,t})$

19
${\boldsymbol{U}}_{j,t+1}={\boldsymbol{U}}_{j,t}+({\boldsymbol{x}}_{t+1}-{\boldsymbol{z}}_{t+1})/\gamma_{t}$

20

21 Variant 1

22

23 $\gamma_{t+1}=\gamma_{t}$

24

25 Variant 2 $\triangleright$ only if each $h_{j}$ is $\beta_{h}$-Lipschitz

26

27 Choose any
$\gamma_{t+1}\in[\gamma_{t},\sqrt{\gamma_{t}^{2}+\gamma_{t}\delta_{t}(2\beta_{h})^{-2}}]$

28

return _${{\boldsymbol{x}}}_{t+1}$ , ${{\boldsymbol{U}}}_{t+1}$_

Algorithm 2 Algo 1 Algorithm 3 Adaptive $(k+1)$-Operator Splitting

#### Discussion.

The above algorithm requires the storage of a vector of size $p$ and two
matrices of size $k\times p$. Hence, for three operators ($k=2$) this requires
to store $2p$ more elements than the three previous formulation, which only
requires to store three arrays of size $p$. Hence, for $k=2$, the formulation
of Algorithm 1 should be preferred. Of course, for $k>2$, it is necessary to
use this formulation.

##  Appendix Appendix C Analysis

In this section we provide the proofs for the theorems stated in the Analysis
section without proof. The appendix is organized as follows:

  * •

In Appendix C.1 we prove the fixed point characterization of our three
operator splitting variant. Its proof is mostly independent of the other
results.

  * •

In Appendix C.2 we prove an inequality that relates the saddle point
suboptimality to the current and previous iterates. This inequality forms the
core of both linear and sublinear convergence proofs. Because of its
importance we name it “Key recursive inequality”.

  * •

Finally, in Appendix C.3 and Appendix C.4 we make use of the previous “key
recursive inequality” to prove the sublinear (Theorem 2 and Corollary 1) and
linear (Theorem 3) convergence results.

Throughout this section we assume that Assumptions 1 and 2 of §3 are verified
without explicit mention.

### Appendix C.1 Fixed point characterization

In this subsection we provide a proof for Theorem 1. Davis & Yin (2017, Lemma
2.2) proved a weaker result that characterized only the first coordinate: in
our notation they proved that if
$({\boldsymbol{x}},{\boldsymbol{u}})\in\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})$,
then ${\boldsymbol{x}}$ is a primal solution. Our theorem extends this results
into a full characterization of the operator.

######  Theorem 1.

Let $\mathcal{P}^{\star}$ denote the set of minimizers of the primal objective
(10) and $\mathcal{D}^{\star}$ the set of minimizers of the dual objective
(11) . Then the set of fixed points of ${\boldsymbol{T}}_{\gamma}$ is given by

$\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})=\mathcal{P}^{\star}\times\mathcal{D}^{\star}~{}.$
(35)

###### Proof.

We will find it useful to first characterize the fixed points of
$\boldsymbol{T}_{\gamma}$ by a subdifferential inclusion. Consider the
following set of equivalences

$\displaystyle({\boldsymbol{z}},{\boldsymbol{u}})\in\operatorname*{{Fix}}({\boldsymbol{T}}_{\gamma})$
$\displaystyle\iff({\boldsymbol{z}},{\boldsymbol{u}})={\boldsymbol{T}}_{\gamma}({\boldsymbol{z}},{\boldsymbol{u}})$
(36)
$\displaystyle\iff\begin{cases}{\boldsymbol{z}}=\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{x}}({\boldsymbol{z}},{\boldsymbol{u}})+\gamma{\boldsymbol{u}})\\\
{\boldsymbol{u}}={\boldsymbol{u}}+({\boldsymbol{x}}({\boldsymbol{z}},{\boldsymbol{u}})-{\boldsymbol{z}})/\gamma\\\
\text{ with
}{\boldsymbol{x}}({\boldsymbol{z}},{\boldsymbol{u}})=\operatorname{\mathbf{prox}}_{\gamma
g}({\boldsymbol{z}}-\gamma({\boldsymbol{u}}+\nabla
f({\boldsymbol{z}})))\end{cases}$ (by definition of $\boldsymbol{T}_{\gamma}$)
(37)
$\displaystyle\iff\begin{cases}{\boldsymbol{z}}=\operatorname{\mathbf{prox}}_{\gamma
g}({\boldsymbol{z}}-\gamma({\boldsymbol{u}}+\nabla f({\boldsymbol{z}})))\\\
{\boldsymbol{z}}=\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{z}}+\gamma{\boldsymbol{u}})\\\ \end{cases}$
(${\boldsymbol{z}}={\boldsymbol{x}}({\boldsymbol{z}},{\boldsymbol{u}})$ by
second equation) (38) $\displaystyle\iff\begin{cases}-{\boldsymbol{u}}-\nabla
f({\boldsymbol{z}})\in\partial g({\boldsymbol{z}})\\\
{\boldsymbol{u}}\in\partial h({\boldsymbol{z}})\end{cases}$ (by Lemma 1) (39)
$\displaystyle\iff\begin{cases}-{\boldsymbol{u}}\in\partial(f+g)({\boldsymbol{z}})\\\
{\boldsymbol{u}}\in\partial h({\boldsymbol{z}})\end{cases}$ (40)
$\displaystyle\iff\begin{cases}{\boldsymbol{z}}\in\partial(f+g)^{*}(-{\boldsymbol{u}})\\\
{\boldsymbol{z}}\in\partial h^{*}({\boldsymbol{u}})\end{cases}$ (by Lemma 2)
(41)

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

_Part 1_. Our goal is to prove
$\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})\subseteq{\mathcal{P}^{\star}\times\mathcal{D}^{\star}}\,$.
Let
$({\boldsymbol{z}},{\boldsymbol{u}})\in\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})$.
Adding together the equations in (40) one obtains

$0\in\partial h({\boldsymbol{z}})+\partial g({\boldsymbol{z}})+\nabla
f({\boldsymbol{z}})~{},$ (42)

and so ${\boldsymbol{z}}$ is a minimizer of the primal objective. Subtracting
the equations in (41) we also have

$0\in\partial
h^{*}({\boldsymbol{u}})-\partial(f+g)^{*}(-{\boldsymbol{u}})~{},$ (43)

and so ${\boldsymbol{u}}$ is a minimizer of the dual objective. We have proved
$\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})\subseteq\mathcal{P}^{\star}\times\mathcal{D}^{\star}\,$.

_Part 2_. Our goal now is to prove the inverse inclusion,
${\mathcal{P}^{\star}\times\mathcal{D}^{\star}}\subseteq\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})$.
Let
$({\boldsymbol{z}},{\boldsymbol{u}})\in\mathcal{P}^{\star}\times\mathcal{D}^{\star}$,
we will prove that $({\boldsymbol{z}},{\boldsymbol{u}})$ is a fixed point of
${\boldsymbol{T}}_{\gamma}$.

We start by recalling the notion of _paramonotonicity_ , which will play a key
role in this 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*}\widetilde{\boldsymbol{a}}\in{\boldsymbol{C}}\boldsymbol{a}\\\
\widetilde{\boldsymbol{b}}\in{\boldsymbol{C}}\boldsymbol{b}\\\
\langle\widetilde{\boldsymbol{a}}-\widetilde{\boldsymbol{b}},\boldsymbol{a}-\boldsymbol{b}\rangle=0\end{rcases*}\implies\widetilde{\boldsymbol{a}}\in{\boldsymbol{C}}\boldsymbol{b}\text{
and }\widetilde{\boldsymbol{b}}\in{\boldsymbol{C}}\boldsymbol{a}~{}.$ (44)

It is known 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}}))$ (45)
$\displaystyle{\boldsymbol{z}}_{\boldsymbol{u}}\in\partial
h^{*}({\boldsymbol{u}})\cap(\partial(f+g)^{*}(-{\boldsymbol{u}}))~{},$ (46)

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}}))~{}.$
(47)

Using Eq. (45) and (47) 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$ (48)

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}})$ (49)

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}})$
(50)

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 h({\boldsymbol{z}})~{}.\end{cases}$ (51)

These are the same subdifferential inclusions of (40), which are equivalent to
$({\boldsymbol{z}},{\boldsymbol{u}})\in\operatorname*{{Fix}}(\boldsymbol{T}_{\gamma})$
(Eq. (36)). This concludes the proof. ∎

### Appendix C.2 Key recursive inequality

In this subsection we present a lemma that will be key for further proofs, as
it relates the saddle point suboptimality to the current and previous
iterates.

######  Lemma 4 (Key recursive inequality).

Let $f$ be $\mu_{f}$-strongly convex and $h^{*}$ be $\mu_{h}$-strongly convex
(where we allow $\mu_{f}=0$ and/or $\mu_{h}=0$). Then after $t$ iterations of
Algorithm 1 we have the following inequality for all
$({\boldsymbol{x}},{\boldsymbol{u}})$ in the domain of $\mathcal{L}$, with
$\gamma_{-1}=\gamma_{0}\,,~{}\delta_{-1}=0$.

For Variant 1:

$2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1}))+\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}+(1-\sigma)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}\\!-\\!{\boldsymbol{u}})\|^{2}~{}.$
(52)

For Variant 2:

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1}))$
$\displaystyle+\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}+(1-\xi)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}$
(53)
$\displaystyle\qquad\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}\\!-\\!{\boldsymbol{u}})\|^{2}+\gamma_{t-1}\delta_{t-1}-2\gamma_{t}\delta_{t}~{}.$

###### Proof.

The proof is structured in two parts. In the first part, we will bound the
“primal” suboptimality
${\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t+1})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1})}$
and in the second part we will bound the “dual” suboptimality
${\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t+1})}$.
Finally, adding both will yield the desired result.

_Part 1_. By $\mu_{f}$-strong convexity $f$ verifies the following inequality
for an arbitrary ${\boldsymbol{x}}$

$f({\boldsymbol{z}}_{t})-f({\boldsymbol{x}})\leq\langle\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\rangle-\frac{\mu_{f}}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}~{}.$
(54)

By the line search condition and the definition of $\delta_{t}$ (we recall
$\delta_{t}=Q_{t}({\boldsymbol{x}}_{t+1})-f({\boldsymbol{x}}_{t+1})$) we have
at each iteration
$f({\boldsymbol{x}}_{t+1})=Q_{t}({\boldsymbol{x}}_{t+1})-\delta_{t}$. This
gives by definition of $Q_{t}$:

$f({\boldsymbol{x}}_{t+1})=f({\boldsymbol{z}}_{t})+\langle\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t+1}-{\boldsymbol{z}}_{t}\rangle+\frac{1}{2\gamma_{t}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t+1}\|^{2}-\delta_{t}~{}\\\
$ (55)

Adding the previous two equations gives the following inequality, which we
will make use of later on

$f({\boldsymbol{x}}_{t+1})-f({\boldsymbol{x}})\,\leq\langle\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}\rangle\,+\frac{1}{2\gamma_{t}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t+1}\|^{2}-\frac{\mu_{f}}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\delta_{t}~{}.$
(56)

From the subdifferential characterization of the proximal operator (Lemma 1),
the update
${\boldsymbol{x}}_{t+1}={\operatorname{\mathbf{prox}}_{\gamma_{t}g}({\boldsymbol{z}}_{t}-\gamma_{t}({\boldsymbol{u}}_{t}+\nabla
f({\boldsymbol{z}}_{t})))}$ of Line 1 implies the following subdifferential
inclusion

$\frac{1}{\gamma_{t}}({\boldsymbol{z}}_{t}-\gamma_{t}({\boldsymbol{u}}_{t}+\nabla
f({\boldsymbol{z}}_{t}))-{\boldsymbol{x}}_{t+1})\in\partial
g({\boldsymbol{x}}_{t+1})~{}.$ (57)

By of $g$ we then have the following inequality

$\displaystyle g({\boldsymbol{x}}_{t+1})-g({\boldsymbol{x}})$
$\displaystyle\leq\frac{1}{\gamma_{t}}\langle{\boldsymbol{z}}_{t}-\gamma_{t}({\boldsymbol{u}}_{t}+\nabla
f({\boldsymbol{z}}_{t}))-{\boldsymbol{x}}_{t+1},{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}\rangle$
(58)
$\displaystyle=\frac{1}{\gamma_{t}}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t+1},{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}\rangle-\langle{\boldsymbol{u}}_{t}+\nabla
f({\boldsymbol{z}}_{t}),{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}\rangle~{}.$
(59)

Adding together (56) and (59) we obtain

$\displaystyle
f({\boldsymbol{x}}_{t+1})+g({\boldsymbol{x}}_{t+1})-f({\boldsymbol{x}})-g({\boldsymbol{x}})$
$\displaystyle\leq\frac{1}{\gamma_{t}}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t+1},{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}\rangle+\frac{1}{2\gamma_{t}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t+1}\|^{2}$
(60)
$\displaystyle\qquad\qquad-\langle{\boldsymbol{u}}_{t},{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}\rangle-\frac{\mu_{f}}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\delta_{t}~{}.$

We will now use this last inequality to bound
$\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t})$.
For this, we will make extensive of the cosine identity
$2\langle\boldsymbol{a},\boldsymbol{b}\rangle={\|\boldsymbol{a}+\boldsymbol{b}\|^{2}}-\|\boldsymbol{a}\|^{2}-\|\boldsymbol{b}\|^{2}$
and the primal-dual relationship
$\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{t})={\boldsymbol{x}}_{t+1}-{\boldsymbol{z}}_{t}$,
which is an immediate consequence of the definition of
${\boldsymbol{u}}_{t+1}$ in Line 1.

$\displaystyle\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t+1})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1})=f({\boldsymbol{x}}_{t+1})+g({\boldsymbol{x}}_{t+1})-f({\boldsymbol{x}})-g({\boldsymbol{x}})+\langle{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}},{\boldsymbol{u}}_{t+1}\rangle$
(61)
$\displaystyle\qquad\stackrel{{\scriptstyle\eqref{eq:primal_part_ineq}}}{{\leq}}\frac{1}{\gamma_{t}}\langle{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t+1},{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}\rangle+\frac{1}{2\gamma_{t}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}_{t+1}\|^{2}+\langle{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}},{\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{t}\rangle-\frac{\mu_{f}}{2}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\delta_{t}$
(62)
$\displaystyle\qquad=\left(\frac{1-\gamma_{t}\mu_{f}}{2\gamma_{t}}\right)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\frac{1}{2\gamma_{t}}\|{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}\|^{2}-\frac{1}{\gamma_{t}}\langle{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}},{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}_{t+1}\rangle-\delta_{t}$
(63) (cosine identity on first term and primal-dual relationship on second-
last term)
$\displaystyle\qquad=\left(\frac{1-\gamma_{t}\mu_{f}}{2\gamma_{t}}\right)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\frac{1}{2\gamma_{t}}\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}+\frac{1}{2\gamma_{t}}\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}_{t+1}\|^{2}-\delta_{t}$
(64) (cosine identity on the second-last term)

_Part 2_. From the subdifferential characterization of the proximal operator
(Lemma 1), the update
${\boldsymbol{z}}_{t+1}={\operatorname{\mathbf{prox}}_{\gamma_{t}h}({\boldsymbol{x}}_{t+1}+\gamma_{t}{\boldsymbol{u}}_{t})}$
of Line 1 we have the following subdifferential inclusion

${\boldsymbol{u}}_{t+1}=\frac{1}{\gamma_{t}}({\boldsymbol{x}}_{t+1}-{\boldsymbol{z}}_{t+1})+{\boldsymbol{u}}_{t}\in\partial
h({\boldsymbol{z}}_{t+1})\implies{\boldsymbol{z}}_{t+1}\in\partial
h^{*}({\boldsymbol{u}}_{t+1})~{},$ (65)

where the implication is a consequence of the conjugate-inverse identity
(Lemma 2). As we did before, this inclusion can be used to obtain an
inequality in terms of the function values. By $\mu_{h}$-strong convexity of
$h^{*}$ we have

$h^{*}({\boldsymbol{u}}_{t+1})-h^{*}({\boldsymbol{u}})\leq\langle{\boldsymbol{z}}_{t+1},{\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}\rangle-\frac{\mu_{h}}{2}\|{\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}\|^{2}~{}.$
(66)

We can now use this inequality to to bound
$\mathcal{L}({\boldsymbol{x}}_{t},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}}_{t},{\boldsymbol{u}}_{t+1})$
as follows:

$\displaystyle\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t+1})=h^{*}({\boldsymbol{u}}_{t+1})-h^{*}({\boldsymbol{u}})+\langle{\boldsymbol{x}}_{t+1},{\boldsymbol{u}}-{\boldsymbol{u}}_{t+1}\rangle$
(67)
$\displaystyle\qquad\stackrel{{\scriptstyle\eqref{eq:dual_part_ineq}}}{{\leq}}\langle{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}\rangle-\frac{\mu_{h}}{2}\|{\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}\|^{2}$
(68)
$\displaystyle\qquad=\frac{1}{\gamma_{t}}\langle{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}_{t+1},\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\rangle-\frac{\mu_{h}}{2}\|{\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}\|^{2}$
(69)
$\displaystyle\qquad=\frac{1}{\gamma_{t}}\langle{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}_{t+1},\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\rangle-\frac{\gamma_{t}^{-1}\mu_{h}}{2\gamma_{t}}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}$
(70)
$\displaystyle\qquad=\frac{1}{2\gamma_{t}}\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}-\frac{1}{2\gamma_{t}}\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}_{t+1}\|^{2}-\Big{(}\frac{1+\gamma_{t}^{-1}\mu_{h}}{2\gamma_{t}}\Big{)}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}~{},$
(71)

where in the second line we have used again the primal-dual relationship
${\boldsymbol{x}}_{t+1}-{\boldsymbol{z}}_{t+1}=\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{t})$.

_Third part: putting it all together_. Adding the inequalities from Eq. (64)
and Eq. (71) and multiplying everything by by $2\gamma_{t}$ we obtain

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1}))$
$\displaystyle\leq(1-\gamma_{t}\mu_{f})\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}$
$\displaystyle\qquad+\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}-(1+\gamma_{t}^{-1}\mu_{h})\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}-2\gamma_{t}\delta_{t}~{}.$
(72)

We will now prove that $\gamma_{t}$ is lower bounded by $\min\\{\tau
L_{f}^{-1},\gamma_{0}\\}$ through a distinction of cases. If
$\gamma_{0}\geq\tau/L_{f}$, then by the properties of $L_{f}$-smooth
functions, the sufficient decrease condition is verified for all step-size
smaller than $1/L_{f}$. The resulting step-size can still be smaller than this
quantity if the sufficient decrease condition fails for some step-size
$\gamma>1/L_{f}$ and $\tau$ is small enough such that $\tau\gamma\leq
1/L_{f}$. Even in this worst case scenario we have $\gamma_{t}\geq\tau/L_{f}$,
which proves the bound for the case $\gamma_{0}\geq\tau/L_{f}$. If
$\gamma_{0}\leq\tau/L_{f}$, then the sufficient decrease condition is verified
at this first iterate, and the step-size can only increase.

We can hence bound $\gamma_{t}\mu_{f}$ by $\rho$, defined in Eq. (16), to
obtain

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1}))$
$\displaystyle\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}$
$\displaystyle\qquad+\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}-(1+\gamma_{t}^{-1}\mu_{h})\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}-2\gamma_{t}\delta_{t}~{}.$
(73)

For Variant 1, we can just drop the non-positive term $-2\gamma_{t}\delta_{t}$
and use the non-increasing step to bound the following terms as:

$\displaystyle\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}$
$\displaystyle\leq\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}\quad\text{
(with $\gamma_{-1}=\gamma_{0}$ by definition)}$ (74)
$\displaystyle-(1+\gamma_{t}^{-1}\mu_{h})$
$\displaystyle\leq-(1+\gamma_{0}^{-1}\mu_{h})$ (75)
$\displaystyle=-\frac{\gamma_{0}+\mu_{h}}{\gamma_{0}}=-(\frac{\gamma_{0}+\mu_{h}-\mu_{h}}{\gamma_{0}+\mu_{h}})^{-1}$
(76)
$\displaystyle=-(1-\frac{\mu_{h}}{\gamma_{0}+\mu_{h}})^{-1}=-(1-\sigma)^{-1}$
(77)

Replacing in (Appendix C.2) gives

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1}))$
$\displaystyle\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}$
$\displaystyle\qquad+\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}-(1-\sigma)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}~{},$
(78)

from where we obtain the desired bound (52) by reordering the terms.

We will now derive a similar bound for Variant 2. In this case, the
$\beta_{h}$-Lipschitz assumption on $h$ implies that the norm of every element
in $\operatorname*{dom}h^{*}$ is bounded by $\beta_{h}$ (see e.g.,
(Rockafellar, 1997, Corollary 13.3.3)). This way we bound
$\|{\boldsymbol{u}}_{t}-{\boldsymbol{u}}\|^{2}\leq
2\|{\boldsymbol{u}}_{t}\|^{2}+2\|{\boldsymbol{u}}\|^{2}\leq 4\beta_{h}$.

Assuming first $t>0$, we have the following sequence of inequalities for any
$\gamma_{t}\geq\gamma_{t-1}$:

$\displaystyle\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}-\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}-\gamma_{t-1}\delta_{t-1}$
$\displaystyle=(\gamma_{t}^{2}-\gamma_{t-1}^{2})\|{\boldsymbol{u}}_{t}-{\boldsymbol{u}}\|^{2}-\gamma_{t-1}\delta_{t-1}$
(79)
$\displaystyle\leq(\gamma_{t}^{2}-\gamma_{t-1}^{2})4\beta_{h}^{2}-\gamma_{t-1}\delta_{t-1}$
(80) (using $\gamma_{t}\geq\gamma_{t-1}$ and that $h$ is
$\beta_{h}$-Lipschitz)
$\displaystyle\leq(\gamma_{t-1}^{2}+\gamma_{t-1}\delta_{t-1}\beta_{h}^{-2}/4-\gamma_{t-1}^{2})4\beta_{h}^{2}-\gamma_{t-1}\delta_{t-1}$
(81) (by the choice of $\gamma_{t+1}$ in Line 1)
$\displaystyle=\gamma_{t-1}^{2}\beta_{h}^{2}+\gamma_{t-1}\delta_{t-1}-\gamma_{t-1}^{2}\beta_{h}^{2}-\gamma_{t-1}\delta_{t-1}$
(82) $\displaystyle=0~{},$ (83)

which reordering gives

$\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}\leq\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}+\gamma_{t-1}\delta_{t-1}~{}.$
(84)

This inequality is also trivially true when $\gamma_{t}\leq\gamma_{t-1}$
because of the non-negativity of $\delta_{t}$. If $t=0$, then we have by
definition $\gamma_{0}=\gamma_{-1},\delta_{-1}=0$ and so the above inequality
is also trivially verified.

Finally, plugging this lasts bound in (Appendix C.2) gives

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1}))$
$\displaystyle\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}$
$\displaystyle\quad+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}-(1+\gamma_{t}^{-1}\mu_{h})\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}-2\gamma_{t}\delta_{t}+\gamma_{t-1}\delta_{t-1}~{}.$

Now we will use the following bound on $\gamma_{t}^{-1}$:
$\mu_{f}\leq\gamma_{t}^{-1}$. This bound can be deduced from the the strong
convexity inequality (Definition 4), as otherwise the sufficient decrease
condition $\delta_{t}\geq 0$ would not hold. In all, we have the desired bound

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t+1}))$
$\displaystyle\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}\|^{2}-\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}})\|^{2}$
$\displaystyle\qquad-(1+\mu_{f}\mu_{h})\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}+\gamma_{t-1}\delta_{t-1}-2\gamma_{t}\delta_{t}~{}.$

Finally, by trivial algebraic manipulations we have

$(1-\xi)^{-1}=(1-\frac{\mu_{f}\mu_{h}}{1+\mu_{f}\mu_{h}})^{-1}=(\frac{1}{1+\mu_{f}\mu_{h}})^{-1}={1+\mu_{f}\mu_{h}}$
(85)

which is the term that multiplies
$\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}$ in the previous
equation. Replacing with $(1-\xi)^{-1}$ in that equation we have the desired
bound.

∎

### Appendix C.3 Sublinear convergence

######  Theorem 2 (sublinear convergence rate).

For every iteration ${t\geq 0}$ and any $({\boldsymbol{x}},{\boldsymbol{u}})$
in the domain of $\mathcal{L}$ we have the following convergence rate for
Algorithm 1 (both variants):

$\mathcal{L}(\overline{\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},\overline{\boldsymbol{u}}_{t+1})\leq\frac{{\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}}+{\gamma_{0}^{2}\|{\boldsymbol{u}}_{0}-{\boldsymbol{u}}\|^{2}}}{2s_{t}}~{}.$
(86)

###### Proof.

Adding the Equation of Lemma 4 with $\mu_{f}=\mu_{h}=0$ (which implies
$\rho=\sigma=\xi=0$) from $0$ to $t$ and dropping positive terms in the left
hand side we get for both Variant 1:

$\displaystyle\sum_{i=0}^{t}\gamma_{i}\left(\mathcal{L}({\boldsymbol{x}}_{i+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{i+1})\right)$
$\displaystyle\leq\frac{1}{2}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+\frac{1}{2}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}})\|^{2}$
$\displaystyle\qquad-\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}-\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}$
(87)
$\displaystyle\leq\frac{1}{2}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+\frac{1}{2}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}})\|^{2}~{}.$
(88)

Similarly, for Variant 2 we have

$\displaystyle\sum_{i=0}^{t}\gamma_{i}\left(\mathcal{L}({\boldsymbol{x}}_{i+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{i+1})\right)$
$\displaystyle\leq\frac{1}{2}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+\frac{1}{2}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}})\|^{2}-\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}})\|^{2}$
$\displaystyle\qquad\qquad-\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}\|^{2}+\gamma_{-1}\delta_{-1}-\sum_{i=0}^{t}\gamma_{t}\delta_{t}$
(89)
$\displaystyle\leq\frac{1}{2}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+\frac{1}{2}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}})\|^{2}+\gamma_{-1}\delta_{-1}-\sum_{i=0}^{t}\gamma_{t}\delta_{t}$
(90)
$\displaystyle\leq\frac{1}{2}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+\frac{1}{2}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}})\|^{2}~{},$
(91)

where the last line follows from the definition $\delta_{-1}=0$ and the non-
negativity of $\delta_{t}$.

In all, we have for both variants

$\sum_{i=0}^{t}\gamma_{i}\left(\mathcal{L}({\boldsymbol{x}}_{i+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{i+1})\right)\leq\frac{1}{2}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+\frac{1}{2}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}})\|^{2}~{}.$
(92)

For a fixed $({\boldsymbol{x}},{\boldsymbol{u}})$, the saddle point
suboptimality
$\mathcal{L}({\boldsymbol{x}}_{t},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{t})$
is convex in ${\boldsymbol{x}}_{t}$ and ${\boldsymbol{u}}_{t}$. By Jensens
inequality we can then bound the left hand side of Eq. (92) as

$\frac{1}{s_{t}}\sum_{i=0}^{t}\gamma_{i}(\mathcal{L}({\boldsymbol{x}}_{i+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},{\boldsymbol{u}}_{i+1}))\geq\mathcal{L}(\overline{\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},\overline{\boldsymbol{u}}_{t+1})~{}.$
(93)

Combining this last inequality with Eq. (92) we obtain

$\mathcal{L}(\overline{\boldsymbol{x}}_{t+1},{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}},\overline{\boldsymbol{u}}_{t+1})\leq\frac{1}{2s_{t}}\left(\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}})\|^{2}\right)~{},$
(94)

which is the desired result. ∎

######  Corollary 1.

Let $h$ be $\beta_{h}$-Lipschitz. Then, we have the following rate for the
weighted ergodic iterate

$P(\overline{\boldsymbol{x}}_{t+1})\\!-\\!P({\boldsymbol{x}}^{\star})\leq\frac{\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+2\gamma_{0}^{2}(\|{\boldsymbol{u}}_{0}\|^{2}+\beta_{h}^{2})}{2s_{t}}$
(95)

###### Proof.

Let
$\widehat{\boldsymbol{u}}\stackrel{{\scriptstyle\text{def}}}{{=}}\operatorname*{\mathbf{arg\,min}}_{{\boldsymbol{u}}}\mathcal{L}(\overline{\boldsymbol{x}}_{t+1},{\boldsymbol{u}})$
and $({\boldsymbol{x}}^{\star},{\boldsymbol{u}}^{\star})$ be a saddle point of
$\mathcal{L}$. Then
$\mathcal{L}(\overline{\boldsymbol{x}}_{t+1},\widehat{\boldsymbol{u}})=P(\overline{\boldsymbol{x}}_{t+1})$
and
$\mathcal{L}({\boldsymbol{x}}^{\star},{\boldsymbol{u}}^{\star})=P({\boldsymbol{x}}^{\star})$
by definition of Fenchel dual.

Using this and the previous theorem we can write the following set of
inequalities

$\displaystyle P(\overline{\boldsymbol{x}}_{t+1})-P({\boldsymbol{x}}^{\star})$
$\displaystyle=\mathcal{L}(\overline{\boldsymbol{x}}_{t+1},\widehat{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}}^{\star},{\boldsymbol{u}}^{\star})$
(96)
$\displaystyle\leq\mathcal{L}(\overline{\boldsymbol{x}}_{t+1},\widehat{\boldsymbol{u}})-\mathcal{L}({\boldsymbol{x}}^{\star},\widehat{\boldsymbol{u}})$
(97) (definition of saddle point, Eq. (9) with
${\boldsymbol{x}}={\boldsymbol{x}}^{\star}$)
$\displaystyle\leq\frac{1}{2s_{t}}\left(\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}\|^{2}+\|\gamma_{0}({\boldsymbol{u}}_{0}-\widehat{\boldsymbol{u}})\|^{2}\right)$
(98) (Theorem 2 with
${\boldsymbol{x}}={\boldsymbol{x}}^{\star},{\boldsymbol{u}}=\widehat{\boldsymbol{u}}$)

The $\beta_{h}$-Lipschitz assumption on $h$ implies that the norm of every
element in $\operatorname*{dom}h^{*}$ is bounded by $\beta_{h}$ (see e.g.,
(Rockafellar, 1997, Corollary 13.3.3)). This way we bound
$\|\gamma_{0}({\boldsymbol{u}}_{0}-\widehat{\boldsymbol{u}})\|^{2}\leq
2\gamma_{0}^{2}\|{\boldsymbol{u}}_{0}\|^{2}+2\gamma_{0}^{2}\|\widehat{\boldsymbol{u}}\|^{2}\leq
2\gamma_{0}^{2}(\|{\boldsymbol{u}}_{0}\|^{2}+\beta_{h}^{2})$. Plugging this
bound into the last inequality we have the desired bound

$P(\overline{\boldsymbol{x}}_{t+1})-P({\boldsymbol{x}}^{\star})\leq\frac{\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\gamma_{0}^{2}(\|{\boldsymbol{u}}_{0}\|^{2}+\beta_{h}^{2})}{2s_{t}}~{}.$
(99)

∎

### Appendix C.4 Linear convergence

In this subsection we assume that $f$ is $\mu_{f}$-strongly convex and $h$ is
$L_{h}$-smooth (with $\mu_{f}>0,0<L_{h}<+\infty$). We denote by
${\boldsymbol{x}}^{\star}$ the minimizer of the primal loss (unique by strong
convexity of $P$) and by ${\boldsymbol{u}}^{\star}$ the minimizer of the dual
loss (also unique by strong convexity of $D$, which is a consequence of the
duality between $L$-smoothness and strong convexity).

######  Theorem 3.

Let ${\boldsymbol{x}}_{t+1},{\boldsymbol{u}}_{t+1}$ be the iterates produced
by Algorithm 1 after $t$ iterations. Then we have the following linear
convergence for Variant 1 (V1) and Variant 2 (V2):

$\displaystyle\text{V1}\mathrel{\mathop{\ordinarycolon}}\|{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}\leq\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}D_{0}$
(100)
$\displaystyle\text{V2}\mathrel{\mathop{\ordinarycolon}}\|{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}E_{0}~{},$
(101)

with
$D_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}6\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\frac{6}{1-\sigma}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}}^{\star})\|^{2}$
and
$E_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}6\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\frac{6}{1-\xi}\|\gamma_{0}({\boldsymbol{u}}_{0}-{\boldsymbol{u}}^{\star})\|^{2}$.

###### Proof.

Variant 1: By the duality between Lipschitz gradient and strong convexity, $h$
being $L_{h}$-smooth implies that $h^{*}$ is $L_{h}^{-1}$-strongly convex.
Applying Lemma 4 with ${\boldsymbol{x}}={\boldsymbol{x}}^{\star}$,
${\boldsymbol{u}}={\boldsymbol{u}}^{\star}$ we obtain the following inequality
for Variant 1:

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}^{\star})-$
$\displaystyle\mathcal{L}({\boldsymbol{x}}^{\star},{\boldsymbol{u}}_{t+1}))+\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+(1-\sigma)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}$
(102)
$\displaystyle\qquad\qquad\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}~{},$
(103)

where we note that the assumption $0<L_{h}$ implies $\sigma<1$ and so
$(1-\sigma)^{-1}$ is well defined. For convenience we introduce the notation
$\alpha_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+(1-\sigma)^{-1}\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}$.
With this, the previous inequality can be simplified to

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}^{\star})-$
$\displaystyle\mathcal{L}({\boldsymbol{x}}^{\star},{\boldsymbol{u}}_{t+1}))+\alpha_{t+1}\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}~{}.$
(104)

By the definition of saddle point we have
$\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}^{\star})-\mathcal{L}({\boldsymbol{x}}^{\star},{\boldsymbol{u}}_{t+1})\geq
0$. Dropping this non-negative term gives

$\displaystyle\alpha_{t+1}\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}~{}.$
(105)

We now make a distinction of cases based on the relative magnitude of $\rho$
and $\sigma$.

  * •

If $\rho\leq\sigma$ then $(1-\rho)\geq(1-\sigma)\implies
1\leq(1-\rho)(1-\sigma)^{-1}$ and so from the previous equation we have

$\displaystyle\alpha_{t+1}$
$\displaystyle\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}$
(106)
$\displaystyle\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+(1-\rho)(1-\sigma)^{-1}\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}$
(107) $\displaystyle=(1-\rho)\alpha_{t}$ (108)

  * •

Otherwise, if $\sigma<\rho$, then $(1-\rho)<(1-\sigma)$ and we have

$\displaystyle\alpha_{t+1}$
$\displaystyle\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}$
(109)
$\displaystyle\leq(1-\sigma)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}$
(110) $\displaystyle=(1-\sigma)\alpha_{t}$ (111)

Combining the two previous equations we have

$\alpha_{t+1}\leq\Big{(}1-\min\\{\rho,\sigma\\}\Big{)}\alpha_{t}~{},$ (112)

which leads by recurrence to

$\alpha_{t+1}\leq\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}\alpha_{0}~{}.$
(113)

By definition of $\alpha_{t+1}$ the previous inequality gives

$\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+(1-\sigma)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}\leq\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}\alpha_{0}.$
(114)

From where we have the following geometric bounds for the primal and dual
variables

$\displaystyle\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}$
$\displaystyle\leq\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}\alpha_{0}$
(115)
$\displaystyle\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}$
$\displaystyle\leq(1-\sigma)\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}\alpha_{0}$
(116)
$\displaystyle\leq\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+2}\alpha_{0}$
(117)

This would be sufficient to derive a convergence rate in terms of
$\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}$. However, since our
sublinear convergence result was in terms of ${\boldsymbol{x}}_{t+1}$, we
would like to state this result in terms of ${\boldsymbol{x}}_{t+1}$ too. This
is possible with a small amount of work and loosing a constant factor $6$.

Using the primal-dual relationship
${\boldsymbol{x}}_{t+1}={\boldsymbol{z}}_{t+1}+\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{t})$
and an improved version of Young’s inequality (Lemma 3)
$\|\boldsymbol{a}+\boldsymbol{b}+\boldsymbol{c}\|^{2}\leq
3(\|\boldsymbol{a}\|^{2}+\|\boldsymbol{b}\|^{2}+\|\boldsymbol{c}\|^{2})$ we
have

$\displaystyle\|{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}$
$\displaystyle=\|{\boldsymbol{z}}_{t+1}+\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{*})+\gamma_{t}({\boldsymbol{u}}_{*}-{\boldsymbol{u}}_{t})-{\boldsymbol{x}}^{\star}\|^{2}\quad\text{
(adding and substracting ${\boldsymbol{u}}^{\star}$)}$ (118)
$\displaystyle\leq
3\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+3\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{*})\|^{2}+3\|\gamma_{t}({\boldsymbol{u}}_{*}-{\boldsymbol{u}}_{t})\|^{2}\quad\text{
(Young's inequality)}$ (119) $\displaystyle\leq
3\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+3\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+3\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}$
(120) ( $\gamma_{t}\leq\gamma_{t-1}$) $\displaystyle\leq
3\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+3(1-\sigma)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+3\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}$
(121) ($1\leq(1-\sigma)^{-1}$) $\displaystyle\leq
3\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}\alpha_{0}\,+\,3\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}$
(122) (by Eq.(113)) $\displaystyle\leq
3\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}\alpha_{0}\,+\,3\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}\alpha_{0}\quad\text{
(by Eq. \eqref{eq:recursive_bound_u})}$ (123)
$\displaystyle\leq\Big{(}1-\min\big{\\{}\rho,\sigma\big{\\}}\Big{)}^{t+1}6\alpha_{0}~{}.$
(124)

The claimed rate then follows by definition of $\alpha_{0}$.

Variant 2: As in Variant 1, by the duality between Lipschitz gradient and
strong convexity, $h$ being $L_{h}$-smooth implies that $h^{*}$ is
$L_{h}^{-1}$-strongly convex. Using Lemma 4 with
${\boldsymbol{x}}={\boldsymbol{x}}^{\star}$,
${\boldsymbol{u}}={\boldsymbol{u}}^{\star}$ we obtain the following inequality
for Variant 2:

$\displaystyle
2\gamma_{t}(\mathcal{L}({\boldsymbol{x}}_{t+1},{\boldsymbol{u}}^{\star})-$
$\displaystyle\mathcal{L}({\boldsymbol{x}}^{\star},{\boldsymbol{u}}_{t+1}))+\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+(1-\xi)^{-1}\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}+2\gamma_{t}\delta_{t}$
(125)
$\displaystyle\qquad\qquad\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}+\gamma_{t-1}\delta_{t-1}~{}.,$
(126)

where we note that the assumption $0<L_{h}$ implies $\xi<1$ and so
$(1-\xi)^{-1}$ well defined. Once again by the definition of saddle point we
can drop the first non-negative term to get

$\widetilde{\alpha}_{t+1}\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}+\gamma_{t-1}\delta_{t-1}~{},$
(127)

where we noted
$\widetilde{\alpha}_{t}\stackrel{{\scriptstyle\text{def}}}{{=}}\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+(1-\xi)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t-1}-{\boldsymbol{u}}^{\star})\|^{2}+2\gamma_{t-1}\delta_{t-1}$.
Similarly as for the first variant, we can make a distinction of cases to
prove the contraction:

  * •

If $\rho<\xi$, then $(1-\rho)\implies(1-\rho)(1-\xi)\geq 1$ and we have

$\displaystyle\widetilde{\alpha}_{t+1}$
$\displaystyle\leq(1-\rho)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+(1-\rho)(1-\xi)^{-1}\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}+\frac{1}{2}(2\gamma_{t-1}\delta_{t-1})$
(128)
$\displaystyle=\left(1-\min\big{\\{}\rho,\mfrac{1}{2}\big{\\}}\right)\widetilde{\alpha}_{t}$
(129)

  * •

If $\xi<\rho$, then $(1-\rho)<(1-\xi)$ and we have

$\displaystyle\widetilde{\alpha}_{t+1}$
$\displaystyle\leq(1-\xi)\|{\boldsymbol{z}}_{t}-{\boldsymbol{x}}^{\star}\|^{2}+\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}+\frac{1}{2}(2\gamma_{t-1}\delta_{t-1})$
(130)
$\displaystyle=\left(1-\min\big{\\{}\xi,\mfrac{1}{2}\big{\\}}\right)\widetilde{\alpha}_{t}$
(131)

Combining the previous two bounds, we have the following geometric decrease on
the sequence $(\widetilde{\alpha}_{t})$,

$\widetilde{\alpha}_{t+1}\leq\Big{(}1-\min\big{\\{}\rho,\mfrac{1}{2},\xi\big{\\}}\Big{)}\widetilde{\alpha}_{t}~{}.$
(132)

Then by definition of $\widetilde{\alpha}_{t+1}$ we have

$\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+(1-\xi)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+2\gamma_{t}\delta_{t}\leq\Big{(}1-\min\big{\\{}\rho,\sigma,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}\widetilde{\alpha}_{0}.$
(133)

From where we can derive the following geometric bounds for the primal and
dual variables:

$\displaystyle\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}$
$\displaystyle\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}\tilde{\alpha}_{0}$
(134)
$\displaystyle(1-\xi)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+2\gamma_{t}\delta_{t}$
$\displaystyle\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}\widetilde{\alpha}_{0}$
(135)

We will now make a distinction of cases on $\xi$ to derive a more convenient
bound for this last inequality. If
$\xi\leq\frac{1}{2}\implies(1-\xi)^{-1}\leq\frac{1}{2}$ and so we have

$\displaystyle(1-\xi)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+(1-\xi)^{-1}\gamma_{t}\delta_{t}$
$\displaystyle\leq(1-\xi)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+4(1-\xi)^{-1}\gamma_{t}\delta_{t}$
(136) (by non-negativity of $\delta_{t}$)
$\displaystyle\leq(1-\xi)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+2\gamma_{t}\delta_{t}$
(137)
$\displaystyle\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}\widetilde{\alpha}_{0}$
(138)

Multiplying both sides by $(1-\xi)$ and using
$(1-\xi)\leq(1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}})$ we have

$\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+\gamma_{t}\delta_{t}\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+2}\widetilde{\alpha}_{0}~{}.$
(139)

If on the other hand,
$\xi\geq\frac{1}{2}\implies(1-\frac{1}{2})(1-\xi)^{-1}\leq 1$ and so we have

$\displaystyle\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+\gamma_{t}\delta_{t}$
$\displaystyle\leq(1-\mfrac{1}{2})\left((1-\xi)^{-1}\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+2\gamma_{t}\delta_{t}\right)$
(140)
$\displaystyle\leq(1-\mfrac{1}{2})\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}\widetilde{\alpha}_{0}$
(141)
$\displaystyle\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+2}\widetilde{\alpha}_{0}$
(142)

By (138) and (142), we see that in both cases we have the following bound:

$\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+\gamma_{t}\delta_{t}\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+2}\widetilde{\alpha}_{0}~{}.$
(143)

To obtain a convergence rate in terms of
$\|{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}$, we will use the
primal-dual relationship
${{\boldsymbol{x}}_{t+1}={\boldsymbol{z}}_{t+1}+\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{t})}$
and an improved version of Young’s inequality (Lemma 3)
$\|\boldsymbol{a}+\boldsymbol{b}+\boldsymbol{c}\|^{2}\leq
3(\|\boldsymbol{a}\|^{2}+\|\boldsymbol{b}\|^{2}+\|\boldsymbol{c}\|^{2})$. We
have the following sequence of inequalities:

$\displaystyle\|{\boldsymbol{x}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}$
$\displaystyle=\|{\boldsymbol{z}}_{t+1}+\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{*})+\gamma_{t}({\boldsymbol{u}}_{*}-{\boldsymbol{u}}_{t}))-{\boldsymbol{x}}^{\star}\|^{2}$
(144) (adding and substracting ${\boldsymbol{u}}^{\star}$) $\displaystyle\leq
3\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+3\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}_{*})\|^{2}+3\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}_{*})\|^{2}$
(145) (Young’s inequality) $\displaystyle\leq
3\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+3\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+3\|\gamma_{t-1}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}+3(\gamma_{t-1}\delta_{t-1}(2\beta_{h})^{-2})\|{\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star}\|^{2}$
(by the maximum step-size increase
$\gamma_{t}^{2}\leq\gamma_{t-1}^{2}+\gamma_{t-1}\delta_{t-1}(2\beta_{h})^{-2}$)
(146) $\displaystyle\leq
3\|{\boldsymbol{z}}_{t+1}-{\boldsymbol{x}}^{\star}\|^{2}+3\|\gamma_{t}({\boldsymbol{u}}_{t+1}-{\boldsymbol{u}}^{\star})\|^{2}+3\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}+3\gamma_{t-1}\delta_{t-1}$
(147) (by the Lipschitz assumption on $h$ have
$\|{\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star}\|^{2}\leq 4\beta_{h}^{2}$)
(148) $\displaystyle\leq
3\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}\widetilde{\alpha}_{0}\,+\,3\|\gamma_{t}({\boldsymbol{u}}_{t}-{\boldsymbol{u}}^{\star})\|^{2}+3\gamma_{t-1}\delta_{t-1}$
(149) (by Eq.(133) and using the bound $(1-\xi)^{-1}>1$) $\displaystyle\leq
3\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}\widetilde{\alpha}_{0}\,+\,3\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}\widetilde{\alpha}_{0}$
(150) (by Eq. (143))
$\displaystyle\leq\Big{(}1-\min\big{\\{}\rho,\xi,\mfrac{1}{2}\big{\\}}\Big{)}^{t+1}6\widetilde{\alpha}_{0}.$
(151)

The desired rate is then a consequence of the definition of
$\widetilde{\alpha}_{0}$, using that by definition $\delta_{-1}=0$. ∎

### Appendix C.5 Convergence rates comparison

We compare the obtained convergence rate against those in Davis & Yin (2015)
for the fixed step-size strategy.

#### Sublinear convergence.

Our own Corollary 1 and (Davis & Yin, 2015, Corollary D.5.2) provide a
convergence rate of the form

$P(\overline{\boldsymbol{x}}_{t+1})-P({\boldsymbol{x}}^{\star})\leq\frac{L_{f}}{2(t+1)}Q_{0}~{},$
(152)

with a different definition of $Q_{0}$ in both cases:

$\displaystyle\text{This work: }Q_{0}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\frac{2}{L_{f}^{2}}\left(\|{\boldsymbol{u}}_{0}\|^{2}+\beta_{h}^{2}\right)$
(153)
$\displaystyle\text{\cite[citep]{(\@@bibref{AuthorsPhrase1Year}{davis2015three}{\@@citephrase{,
}}{})}: }Q_{0}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\|{\boldsymbol{y}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\frac{1}{2}\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}^{*}\|^{2}+\frac{4}{L_{f}}\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}^{*}\|\|\nabla
f({\boldsymbol{x}}^{\star})\|+\frac{4\beta_{h}}{L_{f}}\|{\boldsymbol{y}}_{0}-{\boldsymbol{y}}^{*}\|~{},$
(154)

where
${\boldsymbol{y}}_{0}\stackrel{{\scriptstyle\text{def}}}{{=}}{\boldsymbol{z}}_{0}+\frac{1}{L_{f}}{\boldsymbol{u}}_{0}$,
${\boldsymbol{y}}^{*}\stackrel{{\scriptstyle\text{def}}}{{=}}{\boldsymbol{x}}^{\star}+\frac{1}{L}{\boldsymbol{u}}^{\star}$
the value of $Q_{0}$ for (Davis & Yin, 2015) was obtained by plugging the
step-size $\gamma=1/L_{f}$ in their Corollary D.5.2, and optimizing with
respect to their $\varepsilon$ parameter.

Both quantities are difficult to compare, but it is instructive to compare
them on specific cases. If $h=0$ and
${\boldsymbol{u}}_{0}={\boldsymbol{u}}^{\star}=0$ we have

$\displaystyle\text{This work: }Q_{0}$
$\displaystyle=\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}$ (155)
$\displaystyle\text{\cite[citep]{(\@@bibref{AuthorsPhrase1Year}{davis2015three}{\@@citephrase{,
}}{})}: }Q_{0}$
$\displaystyle=\frac{3}{2}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|^{2}+\frac{4}{L_{f}}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|\|\nabla
f({\boldsymbol{x}}^{\star})\|+\frac{4\beta_{h}}{L_{f}}\|{\boldsymbol{z}}_{0}-{\boldsymbol{x}}^{\star}\|~{},$
(156)

And in this case it is clear that our bound is better. Consider now the case
$f=g=0$ and ${\boldsymbol{z}}_{0}={\boldsymbol{x}}^{\star}=0$. Then we have
the following value of $Q_{0}$:

$\displaystyle\text{This work: }Q_{0}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{2}{L_{f}^{2}}\left(\|{\boldsymbol{u}}_{0}\|^{2}+\beta_{h}^{2}\right)$
(157)
$\displaystyle\text{\cite[citep]{(\@@bibref{AuthorsPhrase1Year}{davis2015three}{\@@citephrase{,
}}{})}: }Q_{0}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\frac{1}{L_{f}^{2}}\|{\boldsymbol{u}}_{0}\|^{2}+\frac{1}{2L_{f}^{2}}\|{\boldsymbol{u}}_{0}-{\boldsymbol{u}}^{\star}\|^{2}+\frac{4\beta_{h}}{L_{f}^{3/2}}\|{\boldsymbol{u}}_{0}-{\boldsymbol{u}}^{\star}\|~{},$
(158)

in this case however the bound of (Davis & Yin, 2015) can be better as it
avoids the quadratic depency on $\beta_{h}$, although the quantity
$\|{\boldsymbol{u}}_{0}-{\boldsymbol{u}}^{\star}\|$ is also in the worse case
in the order of $2\beta_{h}$.

#### Linear convergence.

For a step-size $\gamma=1/L_{f}$, the convergence rate of (Davis & Yin, 2017,
Theorem D.6, point 6) gives a convergence of the form (using our notation and
the step-size $\gamma=1/L_{f}$):

$\frac{2\mu_{f}(1-\eta)}{L_{f}(1+\gamma L_{h})^{2}}$ (159)

where $\eta$ is a quantity that verifies the inequalities

$\eta>\frac{1}{\varepsilon},\text{ with }\varepsilon>\frac{1}{2}$ (160)

by the choice of step-size (see their Theorem D.6) for a definition of these
quantities. Optimizing with respect to $\eta$ gives $\eta<1/2$ and so the rate
can be arbitrarily close to

$\frac{\mu_{f}}{L_{f}(1+\gamma L_{h})^{2}}=\rho\sigma^{2}~{}.$ (161)

This bound is clearly _strictly worse_ than our $\min\\{\rho,\sigma\\}$ of
Theorem 3 for Variant 1 with $\tau=1$. The constant factor $D_{0}$ is worse
than the one obtained by Davis & Yin (2015), but only by a factor
$6(1-\sigma)^{-1}$, and note also that the rate of Davis & Yin (2015) is not
given in terms of the primal variables but in terms of the less interpretable
quantity
$\|{\boldsymbol{x}}_{t}+\gamma{\boldsymbol{u}}_{t+1}-{\boldsymbol{x}}^{\star}-\gamma{\boldsymbol{u}}^{\star}\|^{2}$.

The difference between our convergence rate and that of Davis & Yin (2015) can
be quite large, as illustrated in the plot below for values in the interval
$[0,\frac{1}{2}]^{2}$:

##  Appendix Appendix D Proximal operators

In this section we provide a more detailed description of the different
splitting methods. We start with a reminder of some penalties that have a
known closed form proximal operator.

  * •

$\ell_{1}$ norm. The proximal operator of the $\ell_{1}$ norm
$\|{\boldsymbol{x}}\|_{1}\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{j=1}^{p}|{\boldsymbol{x}}_{j}|$
norm is given by the soft-thresholding operator. More precisely, for
${\boldsymbol{x}}\in{\mathbb{R}}^{p}$ it is given componentwise as

$[\operatorname{\mathbf{prox}}_{\gamma\|\cdot\|}({\boldsymbol{x}})]_{j}=\begin{cases}\left(1-\frac{\gamma}{|{\boldsymbol{x}}_{j}|}\right){\boldsymbol{x}}_{j}~{},&\text{
if $|{\boldsymbol{x}}_{j}|\geq\gamma$}\\\ 0&\text{otherwise}\end{cases}$ (162)

for $j\in\\{1,\ldots,p\\}$.

  * •

Group lasso norm. Let $\mathcal{G}$ be a partition of $\\{1,\ldots,p\\}$. Then
the group lasso or $\ell_{1}/\ell_{2}$ norm for the $\mathcal{G}$ partition is
defined as
$\|{\boldsymbol{x}}\|_{\mathcal{G}}\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{g\in\mathcal{G}}\|[{\boldsymbol{x}}]_{g}\|_{2}$.
Its proximal operator is given as

$[\operatorname{\mathbf{prox}}_{\gamma\|\cdot\|_{\mathcal{G}}}({\boldsymbol{x}})]_{g}=\begin{cases}\left(1-\frac{\gamma}{\|[{\boldsymbol{x}}]_{g}\|}\right)_{+}[{\boldsymbol{x}}]_{g}~{},&\text{
if $\|[{\boldsymbol{x}}]_{g}\|\geq\gamma$}\\\ 0&\text{otherwise}\end{cases}$
(163)

A proof of this can be found in (Yuan & Lin, 2006, Proposition 1).

### Appendix D.1 $\ell_{1}$ trend filtering.

This penalty can be split as the sum of three proximal terms $h_{1},h_{2}$ and
$h_{3}$ as follows:

$\displaystyle\|{\boldsymbol{x}}\|_{\text{TF}}$
$\displaystyle\stackrel{{\scriptstyle\text{def}}}{{=}}\textstyle\sum_{i=1}^{p-2}|{\boldsymbol{x}}_{i}-2{\boldsymbol{x}}_{i+1}+{\boldsymbol{x}}_{i+2}|$
(164)
$\displaystyle=\underbrace{\sum_{i=1}^{(p-2)/3}|{\boldsymbol{x}}_{i}-2{\boldsymbol{x}}_{i+1}+{\boldsymbol{x}}_{i+2}|}_{=h_{1}({\boldsymbol{x}})}+\underbrace{\sum_{i=1}^{(p-2)/3}|{\boldsymbol{x}}_{i+1}-2{\boldsymbol{x}}_{i+2}+{\boldsymbol{x}}_{i+3}|}_{=h_{2}({\boldsymbol{x}})}+\underbrace{\sum_{i=1}^{(p-2)/3}|{\boldsymbol{x}}_{i+2}-2{\boldsymbol{x}}_{i+3}+{\boldsymbol{x}}_{i+4}|}_{=h_{3}({\boldsymbol{x}})}$
(165)

The term $h_{1}$ can be written as $\|\boldsymbol{L}{\boldsymbol{x}}\|_{1}$
for an $\boldsymbol{L}$ matrix of the form

$\boldsymbol{L}=\begin{bmatrix}1&-2&1&0&0&0&\dots\\\ 0&0&0&1&-2&1&\ddots&\\\
\vdots&\ddots&\\\ 0&0&0&0&0&0&1&-2&1\\\ \end{bmatrix}$ (166)

Usually, penalties of the form $\|\boldsymbol{L}{\boldsymbol{x}}\|_{1}$ do not
have a closed form proximal operator. However, this matrix $\boldsymbol{L}$
has a very special property that will allow to compute its proximal operator
in closed form. In particular, this matrix is a _semi-orthogonal_ matrix,
i.e., it verifies $\boldsymbol{L}\boldsymbol{L}^{T}=\nu\boldsymbol{I}$ for
some $\nu>0$ ($\nu=6$ in our case). For these matrices, the proximal operator
of $\varphi(\boldsymbol{L}{\boldsymbol{x}})$ can be computed as (see e.g.
(Combettes & Pesquet, 2011, Table 1)):

$\operatorname{\mathbf{prox}}_{\varphi(\boldsymbol{L}\cdot)}={\boldsymbol{x}}+\nu^{-1}\boldsymbol{L}^{T}(\operatorname{\mathbf{prox}}_{\nu\varphi}(\boldsymbol{L}{\boldsymbol{x}})-\boldsymbol{L}{\boldsymbol{x}}),$
(167)

where in this case $\varphi$ is the $\ell_{1}$ norm and so its proximal
operator is the soft thresholding operator of Eq. (162).

### Appendix D.2 Isotonic and nearly-isotonic penalties

The isotonic constraint
$\\{{\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}\leq\ldots\leq{\boldsymbol{x}}_{p}\\}$
can be enforced through the use of the indicator function
${\imath\\{{\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}\leq{\boldsymbol{x}}_{3}\leq{\boldsymbol{x}}_{4}\leq\cdots\
\leq{\boldsymbol{x}}_{p}\\}}$. Suppose that $p=2$, i.e., we only have the
constraint $\\{{\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}\\}$. In this case,
the projection of $({\boldsymbol{z}}_{1},{\boldsymbol{z}}_{2})$ can be
computed as follows: if ${\boldsymbol{z}}_{1}\leq{\boldsymbol{z}}_{2}$ the
projection does obviously nothing, otherwise it can be computed as the
projection onto the line generated by the vector $(1,1)$. This gives:

$\operatorname{\mathbf{prox}}_{\imath\\{{\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}\\}}({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})=\begin{cases}({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})&\text{
if ${\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}$}\\\
(({\boldsymbol{x}}_{1}+{\boldsymbol{x}}_{2})/{2},({{\boldsymbol{x}}_{1}+{\boldsymbol{x}}_{2}})/{2})&\text{
otherwise.})\end{cases}$ (168)

Now lets consider the general case. The full indicator function can be
decomposed into a sum of the two following terms

$\imath\\{{\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}\leq{\boldsymbol{x}}_{3}\leq{\boldsymbol{x}}_{4}\leq\cdots\\}=\underbrace{\imath\\{{\boldsymbol{x}}_{1}\\!\leq\\!{\boldsymbol{x}}_{2};{\boldsymbol{x}}_{3}\leq{\boldsymbol{x}}_{4};\cdots\\}}_{=g({\boldsymbol{x}})}+\underbrace{\imath\\{{\boldsymbol{x}}_{2}\leq{\boldsymbol{x}}_{3};{\boldsymbol{x}}_{4}\leq{\boldsymbol{x}}_{5};\cdots\\}}_{=h({\boldsymbol{x}})}~{}.$
(169)

where the both terms are block separable with blocks of size 2. We can hence
use the proximal operator of Eq. (168) to derive the proximal operator of $g$
and $h$, which by block-separability is merely a concatenation of the previous
one. These proximal operators are given block-wise as

$\displaystyle[\operatorname{\mathbf{prox}}_{\gamma
g}({\boldsymbol{z}})]_{\\{2i,2i+1\\}}$
$\displaystyle=\begin{cases}({\boldsymbol{x}}_{2i},{\boldsymbol{x}}_{2i+1})&\text{
if ${\boldsymbol{x}}_{2i}\leq{\boldsymbol{x}}_{2i+1}$}\\\
(({\boldsymbol{x}}_{2i}+{\boldsymbol{x}}_{2i+1})/{2},({{\boldsymbol{x}}_{2i}+{\boldsymbol{x}}_{2i+1}})/{2})&\text{
otherwise.})\end{cases}$ (170)
$\displaystyle[\operatorname{\mathbf{prox}}_{\gamma
h}({\boldsymbol{z}})]_{\\{2i+1,2i+2\\}}$
$\displaystyle=\begin{cases}({\boldsymbol{x}}_{2i+1},{\boldsymbol{x}}_{2i+2})&\text{
if ${\boldsymbol{x}}_{2i+1}\leq{\boldsymbol{x}}_{2i+2}$}\\\
(({\boldsymbol{x}}_{2i+1}+{\boldsymbol{x}}_{2i+2})/{2},({{\boldsymbol{x}}_{2i+1}+{\boldsymbol{x}}_{2i+2}})/{2})&\text{
otherwise.})\end{cases}$ (171)

#### The nearly isotonic penalty.

The nearly isotonic penalty was proposed by Tibshirani et al. (2011) as a
relaxation of the aforementioned isotonic constraints and is defined as

$\|{\boldsymbol{x}}\|_{\text{iso}}\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{i=1}^{p-1}\max\\{{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{i+1},0\\}~{}.$
(172)

The penalty is zero if ${\boldsymbol{x}}_{i}\leq{\boldsymbol{x}}_{i+1}$ and so
it encourages the coefficients to be non-decreasing. It can be split and
similarly as we did for the isotonic constraints:

$\sum_{i=1}^{p-1}\max\\{{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{i+1},0\\}=\underbrace{\sum_{i=1}^{\lfloor
p/2\rfloor}\max\\{{\boldsymbol{x}}_{2i}-{\boldsymbol{x}}_{2i+1},0\\}}_{=g({\boldsymbol{x}})}+\underbrace{\sum_{i=1}^{\lfloor(p-1)/2\rfloor}\max\\{{\boldsymbol{x}}_{2i+1}-{\boldsymbol{x}}_{2i+2},0\\}}_{=h({\boldsymbol{x}})}~{}.$
(173)

Let
$\psi({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})\stackrel{{\scriptstyle\text{def}}}{{=}}\max\\{{\boldsymbol{x}}_{1}-{\boldsymbol{x}}_{2},0\\}$.
If we denote by $\partial_{i}$ the subgradient with respect to the $i$-th
coordinate we have

$\displaystyle\partial_{1}\psi({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})=\begin{cases}0&\text{
if }{\boldsymbol{x}}_{1}<{\boldsymbol{x}}_{2}\\\ [0,1]&\text{ if
}{\boldsymbol{x}}_{1}={\boldsymbol{x}}_{2}\\\ 1&\text{ if
}{\boldsymbol{x}}_{1}>{\boldsymbol{x}}_{2}\end{cases}\quad\text{ and
}\quad\partial_{2}\psi({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})=\begin{cases}0&\text{
if }{\boldsymbol{x}}_{1}<{\boldsymbol{x}}_{2}\\\ [-1,0]&\text{ if
}{\boldsymbol{x}}_{1}={\boldsymbol{x}}_{2}\\\ -1&\text{ if
}{\boldsymbol{x}}_{1}>{\boldsymbol{x}}_{2}\end{cases}$ (174)

then it is easy to verify using the subdifferential inclusion of Lemma 1 that
its proximal operator is given by

$\operatorname{\mathbf{prox}}_{\gamma\psi}({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})=\begin{cases}({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2})&\text{
if ${\boldsymbol{x}}_{1}\leq{\boldsymbol{x}}_{2}$}\\\
({\boldsymbol{x}}_{1}-\gamma,{\boldsymbol{x}}_{2}+\gamma)&\text{ if
${\boldsymbol{x}}_{1}-\gamma\geq{\boldsymbol{x}}_{2}+\gamma$}\\\
(({\boldsymbol{x}}_{1}+{\boldsymbol{x}}_{2})/2,({\boldsymbol{x}}_{1}+{\boldsymbol{x}}_{2})/2)&\text{
otherwise }.\end{cases}$ (175)

The proximal operator of the nearly isotonic penalty can then be computed by
applying block-wise the above proximal operator.

### Appendix D.3 Doubly stochastic constraints

The set of double stochastic matrices is composed of square matrices with
nonnegative entries, each of whose rows and columns sum to $1$, i.e.,
$\\{{\boldsymbol{X}}^{T}\boldsymbol{1}=\boldsymbol{1},\boldsymbol{1}^{T}{\boldsymbol{X}}=\boldsymbol{1},{\boldsymbol{X}}\geq\boldsymbol{0}\\}$.
The indicator function over this set can be split as

$\underbrace{\imath\\{{\boldsymbol{X}}^{T}\boldsymbol{1}=\boldsymbol{1},\boldsymbol{1}^{T}{\boldsymbol{X}}=\boldsymbol{1}\\}}_{=g({\boldsymbol{X}})}+\underbrace{\imath\\{{\boldsymbol{X}}\geq\boldsymbol{0}\\}}_{=h({\boldsymbol{X}})}\quad.$
(176)

For any ${\boldsymbol{X}}\in{\mathbb{R}}^{p\times p}$, the proximal operator
of $g$ and $h$ is then given in closed form as

$\displaystyle\operatorname{\mathbf{prox}}_{\gamma g}({\boldsymbol{X}})$
$\displaystyle={\boldsymbol{X}}+\left(\frac{1}{n}\boldsymbol{I}+\frac{\boldsymbol{1}^{T}{\boldsymbol{X}}\boldsymbol{1}}{n^{2}}\boldsymbol{I}-\frac{1}{n}{\boldsymbol{X}}\right)\boldsymbol{1}\boldsymbol{1}^{T}-\frac{1}{n}\boldsymbol{1}\boldsymbol{1}^{T}{\boldsymbol{X}}$
(177) $\displaystyle\operatorname{\mathbf{prox}}_{\gamma h}({\boldsymbol{X}})$
$\displaystyle=({\boldsymbol{X}}+|{\boldsymbol{X}}|)/2~{},$ (178)

where $|\cdot|$ takes the absolute value componentwise and $\boldsymbol{1}$
denotes the $p$-dimensional vector of ones. The proof of this result can be
found in (Lu et al., 2016).

##  Appendix Appendix E Benchmarks

In this appendix we provide a more detailed description of the experiments, as
well as present an extended set of benchmarks.

  * •

Appendix E.1 discusses the implementation of the different algorithms.

  * •

Appendix E.2 details the synthetic data generation process that we used in the
synthetic benchmarks.

  * •

Each of the experiments we run has a subsection in this appendix with an
extended discussion and set of results. These are organized in the order they
appear in Figure 1: Overlapping group lasso (Appendix E.3), total variation
(Appendix E.4), sparse and low rank matrix recovery with trace norm $\ell_{1}$
penalty (Appendix E.5), and nearly isotonic penalty (Appendix E.6).

### Appendix E.1 Implementation details

Particular care has been taken to ensure a fair and extensive comparison. We
have implemented all algorithms in Python, and used the just-in-time compiler
Numba for some non-vectorizable computationally demanding parts such as
computing the fused lasso proximal operator in the total variation-regularized
problems. We profiled and optimized all algorithms equally, and in some cases
made modifications that depart from the canonical description when this
improved significantly performance (see TOS-AOLS below).

All methods have been implemented in Python and except TOS-AOLS are part of
the C-OPT optimization package (Pedregosa, 2018).

#### PDHG.

The PDHG or Condat-Vũ algorithm (Condat, 2013b; Vũ, 2013) depends on two step-
size parameters $\tau$ and $\sigma$. We parametrize the second step-size as
$\sigma=\beta/\tau$ (this gave a better results than the parametrization
$\sigma=\beta\tau$ used by Malitsky & Pock (2018) for the range $\beta$
considered) then by (Condat, 2013b, Theorem 3.1) $\tau$ needs to be of the
form

$\tau\leq\frac{2(1-\beta)}{L_{f}}~{},~{}\text{ with $\beta<1$ }.$ (179)

With the above value of $\tau$, we tested three different values for $\beta$,
$\beta=0.1,0.5,0.9$, and selected the one that had overall better performance
(we found this to be $\beta=0.5$) hence giving it a slight advantage with
respect to other methods that do not require step-size tuning.

#### Adaptive PDHG.

Although one of the step-sizes is computed by a line search method, due to the
two step-sizes it still requires to select the parametrization constant
$\beta$. This parameter was computed the same way than for PDHG, except in
this case we found the best performing value was $\beta=0.1$. In this method
the step-size is allowed to grow, and the initial value for the step-size in
next iterate needs to be selected from an interval of admissible values.
Following the experimental section of (Malitsky & Pock, 2018), we always take
the maximum admissible step-size.

As it is visible from some of the benchmarks (see e.g., Figures 3, 6, 7), the
method seems to suffer more than other methods from numerical instabilities in
the high precision regime (typically when the suboptimality has reached
$10^{-10}$). We investigated this issue and developed an implementation that
uses Moreau’s decomposition to replaces $\operatorname{\mathbf{prox}}_{\tau
h^{*}}$ by $\operatorname{\mathbf{prox}}_{h/\tau}$, thinking that perhaps the
different scales between $\tau$ and $\beta/\tau$ would be the cause of the
instability, but this did not solve the issue. Interestingly, these
instabilities did not appear on the non adaptive variant, despite both
algorithms share much of the same code.

#### TOS-AOLS.

We implemented the averaged operator line search method of Giselsson et al.
(2016), using TOS as the averaged operator and name it TOS-AOLS. Our
implementation has a crucial improvement with respect to the version detailed
in the original reference that we found was crucial to obtain competitive
results. In our implementation, the line search multiplicative factor
$\alpha_{t}$ (in their notation) is not bound to start from an upper bound
$\alpha_{\text{max}}$ as is described in the reference but is instead
incremented as $\alpha_{t+1}=1.05\alpha_{t}$ in case of success of the line
search condition. This gave much better empirical results.

#### Adaptive TOS.

As described in §2.

### Appendix E.2 Synthetic data generation

For the experiments with synthetic data (except for the image deblurring
task), we followed the same data generation process as in (Agarwal et al.,
2010), which generates a design matrix $\\{\boldsymbol{a}_{i}\\}_{i=1}^{n}$
according to the recursion

$\displaystyle\boldsymbol{a}_{1}$ $\displaystyle={\boldsymbol{z}}_{1}$ (180)
$\displaystyle\boldsymbol{a}_{i}$
$\displaystyle={\boldsymbol{z}}_{i}+p\boldsymbol{a}_{j-1}\text{ for $i>1$}$
(181)

with ${\boldsymbol{z}}_{i}$ sampled from a standardized Gaussian distribution.
$p\in[0,1)$ is a correlation parameter that makes the problem more ill-
conditioned as $p\to 1$. This is necessary as otherwise the resulting dataset
has an unrealistic perfect spectrum. In the experiments we chose $p=0.95$,
which for a matrix of size $65\times 65$ gives a condition number of around
60. Although this value of $p$ might seem high, the resulting condition number
is in fact smaller than the one found in the real datasets. For example, the
real-sim dataset has a condition over 200.

### Appendix E.3 Overlapping group lasso benchmarks