∎

11institutetext: 1 School of Science, Lanzhou University of Technology,
Lanzhou, Gansu, P. R. China. Email<EMAIL_ADDRESS>
2 Department of Mathematics, Nanjing University, Nanjing, P. R. China. This
author was supported by the NSFC grants 11771208 and 11922111. Email:
<EMAIL_ADDRESS>
# A golden ratio primal-dual algorithm for structured convex optimization

Xiaokai Chang1 Junfeng Yang2

(Received: date / Accepted: date)

###### Abstract

We design, analyze and test a golden ratio primal-dual algorithm (GRPDA) for
solving structured convex optimization problem, where the objective function
is the sum of two closed proper convex functions, one of which involves a
composition with a linear transform. GRPDA preserves all the favorable
features of the classical primal-dual algorithm (PDA), i.e., the primal and
the dual variables are updated in a Gauss-Seidel manner, and the per iteration
cost is dominated by the evaluation of the proximal point mappings of the two
component functions and two matrix-vector multiplications. Compared with the
classical PDA, which takes an extrapolation step, the novelty of GRPDA is that
it is constructed based on a convex combination of essentially the whole
iteration trajectory. We show that GRPDA converges within a broader range of
parameters than the classical PDA, provided that the reciprocal of the convex
combination parameter is bounded above by the golden ratio, which explains the
name of the algorithm. An ${\mathcal{O}}(1/N)$ ergodic convergence rate result
is also established based on the primal-dual gap function, where $N$ denotes
the number of iterations. When either the primal or the dual problem is
strongly convex, an accelerated GRPDA is constructed to improve the ergodic
convergence rate from ${\mathcal{O}}(1/N)$ to ${\mathcal{O}}(1/N^{2})$.
Moreover, we show for regularized least-squares and linear equality
constrained problems that the reciprocal of the convex combination parameter
can be extended from the golden ratio to $2$ and meanwhile a relaxation step
can be taken. Our preliminary numerical results on LASSO, nonnegative least-
squares and minimax matrix game problems, with comparisons to some state-of-
the-art relative algorithms, demonstrate the efficiency of the proposed
algorithms.

###### Keywords:

Structured convex optimization saddle point problem primal-dual algorithm
golden ratio acceleration convergence rate fixed point iteration

###### MSC:

49M29 65K10 65Y20 90C25

##  1 Introduction

Let $\mathbb{R}^{p}$ and $\mathbb{R}^{q}$ be finite-dimensional Euclidean
spaces, each endowed with an inner product and the induced norm denoted by
$\langle\cdot,\cdot\rangle$ and $\|\cdot\|=\sqrt{\langle\cdot,\cdot\rangle}$,
respectively. In this paper, we consider structured convex optimization
problem of the form

$\displaystyle\min_{x\in\mathbb{R}^{q}}f(Kx)+g(x),$ (1)

where $f:\mathbb{R}^{p}\rightarrow(-\infty,+\infty]$ and
$g:\mathbb{R}^{q}\rightarrow(-\infty,+\infty]$ are extended real-valued closed
proper convex functions Rockafellar1970Convex , and $K\in\mathbb{R}^{p\times
q}$ is a linear operator from $\mathbb{R}^{q}$ to $\mathbb{R}^{p}$. Problem
(1) arises from numerous applications, including signal and image processing,
machine learning, statistics, mechanics and economics, to name a few, see,
e.g., Chambolle2011A ; Bouwmans2016Handbook ; Yang2011Alternating ;
Hayden2013A ; Bertsekas1982Projection  and the references therein. The saddle
point or primal-dual form of (1) reads

$\displaystyle\min_{x\in\mathbb{R}^{q}}\max_{y\in\mathbb{R}^{p}}~{}~{}\mathcal{L}(x,y):=g(x)+\langle
Kx,y\rangle-f^{*}(y),$ (2)

where $f^{*}(y)=\sup_{u\in\mathbb{R}^{p}}\langle y,u\rangle-f(u)$,
$y\in\mathbb{R}^{p}$, is the Legendre-Fenchel conjugate of $f$. Saddle point
problems are ubiquitous in optimization as they provide a very convenient way
to represent many nonsmooth problems. In particular, (2) is intrinsically
related to the Fenchel dual problem of (1), which is given by

$\displaystyle\max_{y\in\mathbb{R}^{p}}-f^{*}(y)-g^{*}(-K^{\top}y),$ (3)

where $K^{\top}$ denotes the matrix transpose or adjoint operator of $K$.

We note that in this paper we restrict the decision variable $x$ in
$\mathbb{R}^{q}$, yet all the analysis can be extended to any finite
dimensional real Euclidean spaces, e.g., the matrix space $\mathbb{R}^{m\times
n}$ with trace inner product and the induced Frobenius norm. In the rest of
this section, we define some notation, make our assumptions, review some
algorithms for solving (1)-(3) that are closely related to this work, and
summarize our contributions and the organization of this paper.

###  1.1 Notation and assumptions

The matrix and vector transpose operation is denoted by superscript “$\top$”.
The operator norm of $K$ is denoted by $L$, i.e.,
$L:=\|K\|=\|K^{\top}\|=\sup\\{\langle
Kx,y\rangle:\|x\|=\|y\|=1,\,x\in\mathbb{R}^{q},\,y\in\mathbb{R}^{p}\\}$.
Throughout this paper, we denote the golden ratio by $\phi$, i.e.,
$\phi={\sqrt{5}+1\over 2}$. Let $h$ be any extended real-valued closed proper
convex function defined on a finite dimensional Euclidean space
$\mathbb{R}^{m}$. The effective domain of $h$ is denoted by
$\text{dom}(h):=\\{x\in\mathbb{R}^{m}:h(x)<+\infty\\}$, and the
subdifferential of $h$ at $x\in\mathbb{R}^{m}$ is denoted by $\partial
h(x):=\\{\xi\in\mathbb{R}^{m}:\,h(y)\geq
h(x)+\langle\xi,y-x\rangle\text{~{}for all~{}}y\in\mathbb{R}^{m}\\}$.
Furthermore, for $\lambda>0$, the proximal point mapping of $\lambda h$ is
given by

$\displaystyle{\rm Prox}_{\lambda
h}(x):=\arg\min_{y\in\mathbb{R}^{m}}\Big{\\{}h(y)+{1\over
2\lambda}\|y-x\|^{2}\Big{\\}},\quad x\in\mathbb{R}^{m}.$

Since $h$ is closed proper and convex, for any $\lambda>0$, ${\rm
Prox}_{\lambda h}$ is uniquely well defined everywhere. The indicator function
of a set $C$ is denoted by $\iota_{C}(x)$, i.e., $\iota_{C}(x)=0$ if $x\in C$
and $+\infty$ if otherwise. The relative interior of $C$ is denoted by
$\text{ri}(C)$. The identity operator or identity matrix is denoted by $I$,
whose domain or order is clear from the context. The zero vector or matrix is
simply denoted by $0$. The composition of two operators is denoted by
“$\circ$”. The stagnation of two column vectors $u$ and $v$ is also denoted by
$(u;v)$, i.e., $(u;v)=(u^{\top},v^{\top})^{\top}$. The sequence of positive
natural numbers is denoted by ${\mathbb{N}}=\\{1,2,3,\ldots\\}$. Other
notation will be specified later.

Throughout the paper, we make the following blanket assumptions.

###### Assumption 1.1

Assume that the set of solutions of (1) is nonempty and, in addition, there
exists $\tilde{x}\in\text{ri}(\text{dom}(g))$ such that
$K\tilde{x}\in\text{ri}(\text{dom}(f))$.

Under Assumption 1.1, it follows from (Rockafellar1970Convex, , Corollaries
28.2.2 and 28.3.1) that ${\bar{x}}\in\mathbb{R}^{q}$ is a solution of (1) if
and only if there exists ${\bar{y}}\in\mathbb{R}^{p}$ such that
$({\bar{x}},{\bar{y}})$ is a saddle point of ${\mathcal{L}}(x,y)$, i.e.,
${\mathcal{L}}({\bar{x}},y)\leq{\mathcal{L}}({\bar{x}},{\bar{y}})\leq{\mathcal{L}}(x,{\bar{y}})$
for all $(x,y)\in\mathbb{R}^{q}\times\mathbb{R}^{p}$, and furthermore, such
$\bar{y}$ is an optimal solution of the dual problem (3). Throughout this
paper, we denote the set of solutions of (2) by ${\mathcal{S}}$, which is
nonempty under Assumption 1.1 and characterized by

${\mathcal{S}}:=\\{({\bar{x}},{\bar{y}})\in\mathbb{R}^{q}\times\mathbb{R}^{p}:\;0\in\partial
g({\bar{x}})+K^{\top}{\bar{y}}\text{~{}~{}and~{}~{}}0\in\partial
f^{*}({\bar{y}})-K{\bar{x}}\\}.$ (4)

In many applications including signal and image processing and machine
learning, the two component functions in (1) enforce, respectively, data
fitting and regularization. In such cases, $f$ and $g$ usually preserve simple
structures so that their proximal point mappings can be evaluated efficiently.
Examples of such functions are abundant, see, e.g., (Beck2017book, , Chapter
6). We therefore make the following assumption.

###### Assumption 1.2

Assume that the proximal point mappings of the component functions $f$ and $g$
either have closed form formulas or can be evaluated efficiently.

###  1.2 Related algorithms

To solve (1)-(3) simultaneously, one may resort to the well known alternating
direction method of multipliers (ADMM) GM75 ; GabM76 ; Jonathan1992On ;
Lions1979Splitting , the primal-dual algorithm (PDA) Chambolle2011A ;
He2012Convergence ; Pock2011Diagonal ; Esser2010General  and their accelerated
and generalized variants Liu2018Acceleration ; Malitsky2018A . Since the
literature on numerical algorithms for solving (1)-(3) has become so vast, a
thorough overview is far beyond the focus of this paper. Instead, we next
review only some primal-dual type algorithms that are most closely related to
our work.

A main feature of primal-dual type algorithms is that both the primal and the
dual variables are updated at each iteration, and thus the primal and the dual
problems are solved simultaneously. Among others, ADMM GM75 ; GabM76  has been
well studied in the literature Lions1979Splitting ; Jonathan1992On ;
Fazel2013Hankel  and widely used in practice Boyd2010Distributed . The main
difficulty encountered by ADMM when applied to (a reformulation of) (1) is
that a subproblem of the form $\min_{x\in\mathbb{R}^{q}}\frac{1}{2}\|Kx-
b_{n}\|^{2}+g(x)$ needs to be solved at each iteration for some
$b_{n}\in\mathbb{R}^{p}$ varying with the iteration counter $n$. We note that
even though ${\rm Prox}_{g}$ is easy to evaluate, this problem, which involves
a linear operator $K$, has to be solved iteratively in general. On the other
hand, when $f(Kx)$ has a least-squares structure, i.e., $f(Kx)={1\over
2}\|Kx-b\|^{2}$ for some $b\in\mathbb{R}^{p}$, ADMM can be applied to the
equivalent problem $\min_{x,y}\\{f(Kx)+g(y):\;x=y\\}$. In this case, the
$x$-subproblem appears as a least-squares problem and is equivalent to solving
a linear system of equations with coefficient matrix of the form $I+\rho
K^{\top}K$ for some $\rho>0$, which could be expensive or even prohibitive for
large scale problems.

The simplest primal-dual type algorithm for solving (1)-(3), which does not
require to solve any subproblem iteratively or any linear system of equations,
is probably the classical Arrow-Hurwicz method Uzawa58 , which, started at
$(x_{0},y_{0})\in\mathbb{R}^{q}\times\mathbb{R}^{p}$, iterates as

$\displaystyle\left\\{\begin{array}[]{l}x_{n}={\rm Prox}_{\tau g}(x_{n-1}-\tau
K^{\top}y_{n-1}),\vskip 3.0pt plus 1.0pt minus 1.0pt\\\ y_{n}={\rm
Prox}_{\sigma f^{*}}(y_{n-1}+\sigma Kx_{n}),\end{array}\right.$ (7)

for $n\geq 1$. Here $\tau>0$ and $\sigma>0$ are step size parameters. As a
conjugate function, $f^{*}$ is always closed and convex. Furthermore, $f^{*}$
is also proper since $f$ is proper and convex. Thus, ${\rm Prox}_{\sigma
f^{*}}$ is uniquely well defined everywhere. The iterative scheme (7) is also
known as primal dual hybrid gradient method in image processing community, see
ZhC08cam ; Esser2010General ; Chambolle2011A . The main computational cost per
iteration is the evaluation of two proximal point mappings and two matrix-
vector multiplications. The convergence of the Arrow-Hurwicz method was
studied in Esser2010General  with very small step sizes, and
${\mathcal{O}}(1/\sqrt{N})$ rate of convergence, measured by primal-dual gap,
was obtained in Chambolle2011A ; Nedic2009Subgradient  when the domain of
$f^{*}$ is assumed to be bounded. However, the Arrow-Hurwicz method does not
converge in general, see He2012Convergence  for a divergent example.

Chambolle and Pock Chambolle2011A ; Pock2011Diagonal  adopted an extrapolation
step after obtaining $x_{n}$ in (7) to obtain
${\bar{x}}_{n}=x_{n}+\delta(x_{n}-x_{n-1})$ for some $\delta\in[0,1]$, which
is then used to replace $x_{n}$ in (7) in the update of $y_{n}$. The resulting
scheme is nowadays widely accepted as PDA and appears as

$\displaystyle\left\\{\begin{array}[]{l}x_{n}={\rm Prox}_{\tau g}(x_{n-1}-\tau
K^{\top}y_{n-1}),\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
{\bar{x}}_{n}=x_{n}+\delta(x_{n}-x_{n-1}),\vskip 3.0pt plus 1.0pt minus
1.0pt\\\ y_{n}={\rm Prox}_{\sigma f^{*}}(y_{n-1}+\sigma
K{\bar{x}}_{n}).\end{array}\right.$ (11)

In the case $\delta=1$, the convergence of (11) was established in
Chambolle2011A  under the condition $\tau\sigma L^{2}<1$, where $L=\|K\|$.
Note that, for $\delta=1$ the scheme (11) is referred to as a split inexact
Uzawa method by Esser, Zhang and Chan Esser2010General , where the connection
of PDA with preconditioned or linearized ADMM has been revealed, see also
Chambolle2011A ; Shefi2014Rate . Later, it was shown in He2012Convergence
that PDA (11) can be viewed as a weighted proximal point method applied to the
variational inequality (VI) representation of the optimality conditions of
(2). The overrelaxed, inertial and accelerated versions of (11) were
investigated in Chambolle2016ergodic . See also Chambolle2018STOCHASTIC  for a
stochastic variant of PDA.

The effectiveness of the PDA scheme (11) is essentially guaranteed by the
extrapolation or inertial step ${\bar{x}}_{n}=x_{n}+\delta(x_{n}-x_{n-1})$
since the Arrow-Hurwicz method (7), which corresponds to $\delta=0$, fails to
converge in general. For a fixed $\delta\in[0,1]$, the scheme (11) can be
written abstractly as $(x_{n},y_{n})=T(x_{n-1},y_{n-1})$ for some suitably
defined mapping $T$. Therefore, PDA can be viewed as a one-step fixed point
iterative method. Recently, Malitsky Malitsky2019Golden  proposed a fully
adaptive forward-backward type splitting algorithm, called golden ratio
algorithm, for solving mixed VI problem, where a novel convex combination
technique is introduced. It is shown that the golden ratio algorithm preserves
${\mathcal{O}}(1/N)$ ergodic convergence and R-linear convergence under an
error bound condition. The mixed VI problem is to find
$z^{*}\in\mathbb{R}^{m}$ such that

$\theta(z)-\theta(z^{*})+\langle z-z^{*},F(z^{*})\rangle\geq 0,\quad\forall
z\in\mathbb{R}^{m},$ (12)

where $\theta:\mathbb{R}^{m}\rightarrow(-\infty,+\infty]$ is a closed proper
convex function, and $F:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ is a monotone
mapping. Given an initial point $z_{0}\in\mathbb{R}^{m}$ and let
${\bar{z}}_{0}=z_{0}$, a basic version of the golden ratio algorithm
Malitsky2019Golden  applied to (12) iterates for $n\geq 1$ as

$\left\\{\begin{array}[]{rcl}\bar{z}_{n}&=&{\phi-1\over\phi}z_{n-1}+{1\over\phi}{\bar{z}}_{n-1},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ z_{n}&=&{\rm
Prox}_{\tau\theta}(\bar{z}_{n}-\tau F(z_{n-1})),\end{array}\right.$ (13)

where $\tau>0$ is a step size parameter and $\phi$ is the golden ratio. By
induction, $\bar{z}_{n}$ is a convex combination of
$\\{z_{i}:i=0,1,\ldots,n-1\\}$. As a result, the update from $z_{n-1}$ to
$z_{n}$ in (13) is essentially dependent on the whole iteration trajectory,
which is vastly different from the one-step iterative scheme (11). Note that
the primal-dual problem (2) is equivalent to (12) with
$\mathbb{R}^{m}=\mathbb{R}^{q}\times\mathbb{R}^{p}$, endowed with the natural
inner product $\langle(x,y),(u,v)\rangle:=\langle x,u\rangle+\langle
y,v\rangle$ for $(x,y),(u,v)\in\mathbb{R}^{q}\times\mathbb{R}^{p}$,

$\theta(z):=g(x)+f^{*}(y)\text{~{}~{}and~{}~{}}F(z):=(K^{\top}y,-Kx).$ (14)

Given the special structure of $\theta$ and $F$ in (14), the golden ratio
algorithm (13) breaks into

$\left\\{\begin{array}[]{rclrcl}\bar{x}_{n}&=&{\phi-1\over\phi}x_{n-1}+{1\over\phi}{\bar{x}}_{n-1},\;\;x_{n}&=&{\rm
Prox}_{\tau g}(\bar{x}_{n}-\tau K^{\top}y_{n-1}),\vskip 3.0pt plus 1.0pt minus
1.0pt\\\
\bar{y}_{n}&=&{\phi-1\over\phi}y_{n-1}+{1\over\phi}{\bar{y}}_{n-1},\;\;y_{n}&=&{\rm
Prox}_{\tau f^{*}}(\bar{y}_{n}+\tau Kx_{n-1}).\end{array}\right.$ (15)

Apparently, (15) is a Jacobian type algorithm, which fails to fully utilize
the latest available information. Furthermore, the primal and dual step sizes
are identically $\tau$, which offers less flexibility. Given the nice
convergence properties and the promising numerical results of golden ratio
algorithm Malitsky2019Golden , it is desirable to construct a Gauss-Seidel
type golden ratio PDA for the saddle point problem (2), which is able to fully
take advantage of the problem structure. This motivated the current work.

###  1.3 Contributions

We adapt the convex combination technique Malitsky2019Golden  into the Arrow-
Hurwicz scheme (7) to construct a Gauss-Seidel type golden ratio PDA
(abbreviated as GRPDA), which not only preserves nice convergence properties
but also performs favorably in practice. Our main contributions are summarized
below.

  * •

Using the convex combination technique Malitsky2019Golden , we construct a
GRPDA with fixed step size parameters $\tau>0$ and $\sigma>0$. Under the
condition $\tau\sigma L^{2}<\psi\in(1,\phi]$, global convergence and
${\mathcal{O}}(1/N)$ ergodic rate of convergence measured by primal-dual gap
function are established. Since $\tau\sigma L^{2}>1$ is permitted, GRPDA
converges in a broader range of parameters compared to PDA, which typically
requires $\tau\sigma L^{2}<1$. It is further explained that GRPDA is an ADMM-
like method with a proximal term $\frac{1}{2}\|x-x_{n-1}\|_{M}^{2}$, where
$M=\frac{1}{\tau}I-\sigma K^{\top}K$, plus an extra linear term. As a
consequence of the extra linear term, the weighting matrix $M$ is allowed to
be indefinite.

  * •

When either $g$ or $f^{*}$ is strongly convex, an accelerated GRPDA is
constructed to improve the ergodic convergence rate from ${\mathcal{O}}(1/N)$
to ${\mathcal{O}}(1/N^{2})$.

  * •

Let $b\in\mathbb{R}^{p}$. For regularized least-squares problem, i.e.,
$f(\cdot)={1\over 2}\|\cdot-b\|^{2}$, or linear equality constrained problem,
i.e., $f(\cdot)=\iota_{\\{b\\}}(\cdot)$, we show via spectral analysis that
the permitted range of $\psi$ can be extended from $(1,\phi]$ to $(1,2]$ and
meanwhile a relaxation step can be taken as well, where the relaxation
parameter lies in $(0,3/2)$. The analysis for these two special cases is based
on the fixed point theory of averaged operators.

  * •

We carry out numerical experiments on LASSO, nonnegative least-squares and
minimax matrix game problems, with comparisons to some state-of-the-art
algorithms, to demonstrate the favorable performance of the proposed
algorithms.

###  1.4 Organization

The rest of this paper is organized as follows. In Section 2, we summarize
some useful facts and identities and define the primal-dual gap function used
in subsequent analysis. The GRPDA with fixed step sizes is presented in
Section 3, where its convergence and ${\mathcal{O}}(1/N)$ ergodic rate of
convergence are established as well. In Section 4, we present an accelerated
GRPDA that enjoys a faster ${\mathcal{O}}(1/N^{2})$ ergodic rate of
convergence under the assumption that either $g$ or $f^{*}$ is strongly
convex. Section 5 is devoted to the analysis of two special cases, i.e.,
regularized least-squares problem and linear equality constrained problem.
Numerical results in comparison with some state-of-the-art relative algorithms
on LASSO, nonnegative least-squares and minimax matrix game problems are given
in Section 6 to demonstrate the efficiency of the proposed algorithms.
Finally, some concluding remarks are drawn in Section 7.

##  2 Preliminaries

In this section, we summarize some useful facts and identities and define the
primal-dual gap function, which will be useful in our analysis.

###### Fact 2.1

For any extended real-valued closed proper convex function $h$ defined on an
Euclidean space $\mathbb{R}^{m}$, $\lambda>0$ and $x\in\mathbb{R}^{m}$, it
holds that $p={\rm Prox}_{\lambda h}(x)$ if and only if

$\displaystyle\langle
p-x,y-p\rangle\geq\lambda\big{(}h(p)-h(y)\big{)},~{}~{}\forall
y\in\mathbb{R}^{m}.$

###### Fact 2.2

Let $\\{a_{n}\\}_{n\in{\mathbb{N}}}$ and $\\{b_{n}\\}_{n\in{\mathbb{N}}}$ be
two nonnegative real sequences. If there exists an integer $N>0$ such that
$a_{n+1}\leq a_{n}-b_{n}$ for all $n>N$, then $\lim_{n\rightarrow\infty}a_{n}$
exists and $\lim_{n\rightarrow\infty}b_{n}=0$.

The proofs of Facts 2.1 and 2.2 are easily derived and thus are omitted. The
following elementary identities will be used in our analsis. For any
$x,y,z\in\mathbb{R}^{m}$ and $\alpha\in{\mathbb{R}}$, there hold

$\displaystyle\langle x-y,x-z\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\|x-y\|^{2}+\frac{1}{2}\|x-z\|^{2}-\frac{1}{2}\|y-z\|^{2},$
(16) $\displaystyle\|\alpha x+(1-\alpha)y\|^{2}$ $\displaystyle=$
$\displaystyle\alpha\|x\|^{2}+(1-\alpha)\|y\|^{2}-\alpha(1-\alpha)\|x-y\|^{2}.$
(17)

We next define the primal-dual gap function that will be used in our analysis.
Let $(\bar{x},\bar{y})\in{\mathcal{S}}$ be any saddle point of the primal-dual
problem (2). Then, there hold $K\bar{x}\in\partial f^{*}(\bar{y})$ and
$-K^{\top}\bar{y}\in\partial g(\bar{x})$, which are equivalent to

$\displaystyle\left\\{\begin{array}[]{l}P(x):=g(x)-g(\bar{x})+\langle
K^{\top}\bar{y},x-\bar{x}\rangle\geq 0,\quad~{}\forall
x\in\mathbb{R}^{q},\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
D(y):=f^{*}(y)-f^{*}(\bar{y})-\langle K\bar{x},y-\bar{y}\rangle\geq
0,\quad~{}\forall y\in\mathbb{R}^{p}.\end{array}\right.$

The primal-dual gap function is defined by

$\displaystyle G(x,y):=P(x)+D(y)\geq
0,~{}\forall(x,y)\in\mathbb{R}^{q}\times\mathbb{R}^{p}.$ (19)

This primal-dual gap function is also used in, e.g., Malitsky2018A ;
Chambolle2016ergodic . Note that, for fixed
$(\bar{x},\bar{y})\in{\mathcal{S}}$, $P(x)$ and $D(y)$, and thus $G(x,y)$, are
convex. Although $P(\cdot)$, $D(\cdot)$ and $G(\cdot,\cdot)$ depend on
$(\bar{x},\bar{y})$, we do not indicate this dependence in our notation since
it is always clear from the context.

##  3 Golden ratio primal-dual algorithm

In this section, we present our GRPDA with fixed step sizes and establish its
convergence and ${\mathcal{O}}(1/N)$ ergodic convergence rate.

###  3.1 GRPDA with fixed step sizes

Recall that $\phi$ represents the golden ratio and $L=\|K\|$ is the operator
norm of $K$. Below, we introduce our GRPDA with fixed step sizes $\tau>0$ and
$\sigma>0$.

######  Algorithm 3.1 (GRPDA with fixed step sizes)

Step 0.

Let $\tau,\sigma>0$ and $\psi\in(1,\phi]$ be such that $\tau\sigma
L^{2}<\psi$. Choose $x_{0}\in\mathbb{R}^{q}$, $y_{0}\in\mathbb{R}^{p}$. Set
$z_{0}=x_{0}$ and $n=1$.

Step 1.

Compute

$\displaystyle\left\\{\begin{array}[]{rcl}z_{n}&=&\frac{\psi-1}{\psi}x_{n-1}+\frac{1}{\psi}z_{n-1},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ x_{n}&=&{\rm Prox}_{\tau g}(z_{n}-\tau
K^{\top}y_{n-1}),\vskip 3.0pt plus 1.0pt minus 1.0pt\\\ y_{n}&=&{\rm
Prox}_{\sigma f^{*}}(y_{n-1}+\sigma Kx_{n}).\end{array}\right.$ (23)

Step 2.

Set $n\leftarrow n+1$ and return to Step 1.

The GRPDA (23) and the PDA (11) are quite similar, both of which are
modifications of the Arrow-Hurwicz scheme (7). The difference is that PDA (11)
adopts an inertial or extrapolation technique, while GRPDA (23) uses a convex
combination $z_{n}=\frac{\psi-1}{\psi}x_{n-1}+\frac{1}{\psi}z_{n-1}$. To our
knowledge, this convex combination technique with $\psi$ upper bounded by the
golden ratio $\phi$ was initially introduced by Malitsky Malitsky2019Golden
to solve monotone mixed VI problem (12). A direct adaptation of
(Malitsky2019Golden, , Eq. (10)) to the primal-dual problem (2) would give the
algorithm (15), which is a Jacobian type algorithm and has less flexibility
since the primal and the dual step sizes are required to be identical. In
contrast, GRPDA (23) is a combination of Malitsky2019Golden  with the Arrow-
Hurwicz scheme (7), and it is of Gauss-Seidel type since it utilizes the
latest available information. Compared with the direct adaptation (15), the
primal and the dual step sizes in (23) are not necessarily identical, which
offers more flexibility in practice. Moreover, compared with PDA (11), which
requires $\tau\sigma L^{2}<1$, the condition $\tau\sigma
L^{2}<\psi\in(1,\phi]$ required by GRPDA (23) permits larger step sizes.
Apparently, GRPDA has the same per iteration cost as those of the Arrow-
Hurwicz method and the PDA.

###  3.2 Connection with ADMM and PDA

It is well known that the classical PDA can be interpreted as a preconditioned
or linearized ADMM with positive definition proximal term, see Chambolle2011A
; Shefi2014Rate ; Esser2010General . As the main difference between GRPDA and
PDA is that the extrapolation step in (11) is replaced by the convex
combination in (23), it is natural to connect GRPDA with the ADMM. To be
specific, we define the augmented Lagrangian function of
$\min_{x,w}\\{g(x)+f(w):Kx=w\\}$, a reformulation of (1), as

${{\mathcal{L}}}_{\sigma}(x,w,y):=g(x)+f(w)+\langle y,Kx-w\rangle+{\sigma\over
2}\|Kx-w\|^{2},$

where $y\in\mathbb{R}^{p}$ is the Lagrange multiplier and $\sigma>0$ is the
penalty parameter. Then, by using the Moreau decomposition $y={\rm
Prox}_{\lambda f}(y)+\lambda{\rm
Prox}_{{1\over\lambda}f^{*}}({y\over\lambda})$ for any $\lambda>0$ and
$y\in\mathbb{R}^{p}$ and following Chambolle2011A ; Shefi2014Rate ;
Esser2010General , it is easy to show that GRPDA (23) is equivalent to

$\displaystyle\left\\{\begin{array}[]{rcl}z_{n}&=&\frac{\psi-1}{\psi}x_{n-1}+{1\over\psi}z_{n-1},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
x_{n}&=&\arg\min\limits_{x}\left\\{{{\mathcal{L}}}_{\sigma}(x,w_{n-1},y_{n-1})+\ell(x,z_{n},x_{n-1},w_{n-1})+\frac{1}{2}\|x-x_{n-1}\|^{2}_{M}\right\\},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
w_{n}&=&\arg\min\limits_{w}{{\mathcal{L}}}_{\sigma}(x_{n},w,y_{n-1}),\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
y_{n}&=&y_{n-1}+\sigma(Kx_{n}-w_{n}),\end{array}\right.$ (28)

where $\ell(x,z_{n},x_{n-1},w_{n-1})={1\over\tau}\langle
x,z_{n}-x_{n-1}\rangle-\sigma\langle Kx,Kx_{n-1}-w_{n-1}\rangle$,
$M={1\over\tau}I-\sigma K^{\top}K$ and $w_{0}\in\mathbb{R}^{p}$ is arbitrary.
Here the proximal term $\frac{1}{2}\|x-x_{n-1}\|^{2}_{M}$ is used to cancel
out ${\sigma\over 2}\|Kx\|^{2}$ in the $x$-subproblem. Since, by following
Chambolle2011A ; Shefi2014Rate ; Esser2010General , the derivation of (28)
from (23) is standard, we omit the details. By discarding the linear term
$\ell(x,z_{n},x_{n-1},w_{n-1})$ in (28) and using the Moreau decomposition to
carry out a similar reduction, we obtain

$\displaystyle\left\\{\begin{array}[]{rcl}{\bar{y}}_{n-1}&=&2y_{n-1}-y_{n-2},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ x_{n}&=&{\rm Prox}_{\tau g}\big{(}x_{n-1}-\tau
K^{\top}{\bar{y}}_{n-1}\big{)},\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
y_{n}&=&{\rm Prox}_{\sigma f^{*}}(y_{n}+\sigma Kx_{n}),\end{array}\right.$

which is a form of PDA with $\delta=1$. Note that PDA requires $\tau\sigma
L^{2}<1$ to guarantee the positive definiteness of $M$. In comparison, by
introducing the linear term $\ell(x,z_{n},x_{n-1},w_{n-1})$, GRPDA (28) allows
larger step sizes since $\tau\sigma L^{2}<\psi\in(1,\phi]$ suffices for global
convergence. In this case, $M$ is permitted to be indefinite. This broader
convergence region of GRPDA is beneficial in practice.

###  3.3 Convergence results

We next establish the convergence and ${\mathcal{O}}(1/N)$ ergodic convergence
rate of GRPDA (23). For convenience, we let $\beta:=\sigma/\tau$ in the rest
of this section. First, we present a useful lemma.

###### Lemma 3.1

Let $\\{(x_{n},y_{n},z_{n})\\}_{n\in{\mathbb{N}}}$ be the sequence generated
by Algorithm 3.1 from any initial point
$(x_{0},y_{0})\in\mathbb{R}^{q}\times\mathbb{R}^{p}$ and $z_{0}=x_{0}$, and
let $G(\cdot,\cdot)$ be defined as in (19) with any fixed
$(\bar{x},\bar{y})\in{\mathcal{S}}$. Then, it holds for any
$(x,y)\in\mathbb{R}^{q}\times\mathbb{R}^{p}$ that

$\displaystyle\tau G(x_{n},y_{n})$ $\displaystyle\leq$ $\displaystyle\langle
x_{n+1}-z_{n+1},\bar{x}-x_{n+1}\rangle+{1\over\beta}\left\langle
y_{n}-y_{n-1},\bar{y}-y_{n}\right\rangle+\psi\left\langle
x_{n}-z_{n+1},x_{n+1}-x_{n}\right\rangle$ (30) $\displaystyle+\tau\langle
K^{\top}(y_{n}-y_{n-1}),x_{n}-x_{n+1}\rangle.$

Proof. It follows from (23) and Fact 2.1 that

$\displaystyle\langle x_{n+1}-z_{n+1}+\tau
K^{\top}y_{n},~{}\bar{x}-x_{n+1}\rangle\geq\tau\big{(}g(x_{n+1})-g(\bar{x})\big{)},$
(31) $\displaystyle\big{\langle}{1\over\beta}(y_{n}-y_{n-1})-\tau
Kx_{n},~{}\bar{y}-y_{n}\big{\rangle}\geq\tau\big{(}f^{*}(y_{n})-f^{*}(\bar{y})\big{)}.$
(32)

Similar to (31), it holds that

$\displaystyle\langle x_{n}-z_{n}+\tau
K^{\top}y_{n-1},~{}x_{n+1}-x_{n}\rangle\geq\tau\big{(}g(x_{n})-g(x_{n+1})\big{)}.$
(33)

Elementary calculations show that the addition of (31), (32) and (33) gives

$\displaystyle\langle
x_{n+1}-z_{n+1},~{}\bar{x}-x_{n+1}\rangle+{1\over\beta}\big{\langle}y_{n}-y_{n-1},~{}\bar{y}-y_{n}\big{\rangle}+\left\langle
x_{n}-z_{n},~{}x_{n+1}-x_{n}\right\rangle$ $\displaystyle+\tau\langle
K^{\top}(y_{n}-y_{n-1}),~{}x_{n}-x_{n+1}\rangle-\tau\langle
K^{\top}\bar{y},~{}x_{n}-\bar{x}\rangle+\tau\langle
K\bar{x},~{}y_{n}-\bar{y}\rangle$ $\displaystyle\geq$
$\displaystyle\tau\big{(}f^{*}(y_{n})-f^{*}(\bar{y})\big{)}+\tau\big{(}g(x_{n})-g(\bar{x})\big{)},$

which is the same as (30) by considering $x_{n}-z_{n}=\psi(x_{n}-z_{n+1})$
(implied by the first equality in (23)) and following the definition of
$G(\cdot,\cdot)$ in (19). This completes the proof. $\Box$

Now, we are ready to establish the convergence of GRPDA (23).

###### Theorem 3.1

Let $\\{(x_{n},y_{n},z_{n})\\}_{n\in{\mathbb{N}}}$ be the sequence generated
by Algorithm 3.1 from any initial point
$(x_{0},y_{0})\in\mathbb{R}^{q}\times\mathbb{R}^{p}$ and $z_{0}=x_{0}$, and
let $G(\cdot,\cdot)$ be defined as in (19) with any fixed
$(\bar{x},\bar{y})\in{\mathcal{S}}$. Then,
$\\{(x_{n},y_{n})\\}_{n\in{\mathbb{N}}}$ converges to a solution of (2), i.e.,
an element of ${\mathcal{S}}$.

Proof. From Lemma 3.1 and (16), we obtain

$\displaystyle\|x_{n+1}-\bar{x}\|^{2}+{1\over\beta}\|y_{n}-\bar{y}\|^{2}+2\tau
G(x_{n},y_{n})$ (34) $\displaystyle\leq$
$\displaystyle\|z_{n+1}-\bar{x}\|^{2}+{1\over\beta}\|y_{n-1}-\bar{y}\|^{2}+2\tau\langle
K^{\top}(y_{n}-y_{n-1}),~{}x_{n}-x_{n+1}\rangle$
$\displaystyle-\psi\|z_{n+1}-x_{n}\|^{2}+(\psi-1)\|x_{n+1}-z_{n+1}\|^{2}-\psi\|x_{n+1}-x_{n}\|^{2}-{1\over\beta}\|y_{n}-y_{n-1}\|^{2}.$

Since $x_{n+1}={\psi\over\psi-1}z_{n+2}-{1\over\psi-1}z_{n+1}$, it follows
from (17) that

$\displaystyle\|x_{n+1}-\bar{x}\|^{2}$ $\displaystyle=$
$\displaystyle{\psi\over\psi-1}\|z_{n+2}-\bar{x}\|^{2}-{1\over\psi-1}\|z_{n+1}-\bar{x}\|^{2}+{\psi\over(\psi-1)^{2}}\|z_{n+2}-z_{n+1}\|^{2}$
(35) $\displaystyle=$
$\displaystyle{\psi\over\psi-1}\|z_{n+2}-\bar{x}\|^{2}-{1\over\psi-1}\|z_{n+1}-\bar{x}\|^{2}+\frac{1}{\psi}\|x_{n+1}-z_{n+1}\|^{2},$

where the second equality is due to
$z_{n+2}-z_{n+1}={\psi-1\over\psi}(x_{n+1}-z_{n+1})$. By plugging (35) into
(34) and combining terms, we obtain

$\displaystyle{\psi\over\psi-1}\|z_{n+2}-\bar{x}\|^{2}+{1\over\beta}\|y_{n}-\bar{y}\|^{2}+2\tau
G(x_{n},y_{n})$ (36) $\displaystyle\leq$
$\displaystyle{\psi\over\psi-1}\|z_{n+1}-\bar{x}\|^{2}+{1\over\beta}\|y_{n-1}-\bar{y}\|^{2}+2\tau\langle
K^{\top}(y_{n}-y_{n-1}),~{}x_{n}-x_{n+1}\rangle$
$\displaystyle-\psi\|z_{n+1}-x_{n}\|^{2}+(\psi-1-\frac{1}{\psi})\|x_{n+1}-z_{n+1}\|^{2}-\psi\|x_{n+1}-x_{n}\|^{2}-{1\over\beta}\|y_{n}-y_{n-1}\|^{2}$
$\displaystyle\leq$
$\displaystyle{\psi\over\psi-1}\|z_{n+1}-\bar{x}\|^{2}+{1\over\beta}\|y_{n-1}-\bar{y}\|^{2}+2\tau\langle
K^{\top}(y_{n}-y_{n-1}),~{}x_{n}-x_{n+1}\rangle$
$\displaystyle-\psi\|z_{n+1}-x_{n}\|^{2}-\psi\|x_{n+1}-x_{n}\|^{2}-{1\over\beta}\|y_{n}-y_{n-1}\|^{2},$

where the second “$\leq$” follows from $\psi-1-\frac{1}{\psi}\leq 0$ as
$\psi\leq\phi$. Recall that $\beta=\sigma/\tau$. Then, $\tau\sigma L^{2}<\psi$
is equivalent to $\tau<\sqrt{{\psi/\beta}}/L$. Let
$\tau=\zeta\sqrt{{\psi/\beta}}/L$ with $0<\zeta<1$. It follows from Cauchy-
Schwarz inequality and $L=\|K^{\top}\|$ that

$\displaystyle 2\tau\langle K^{\top}(y_{n}-y_{n-1}),~{}x_{n}-x_{n+1}\rangle$
$\displaystyle\leq$ $\displaystyle
2\zeta\sqrt{{\psi/\beta}}\|y_{n}-y_{n-1}\|\|x_{n+1}-x_{n}\|$
$\displaystyle\leq$
$\displaystyle\zeta\Big{(}\psi\|x_{n+1}-x_{n}\|^{2}+{1\over\beta}\|y_{n}-y_{n-1}\|^{2}\Big{)},$

which together with (36) leads to

$\displaystyle 2\tau G(x_{n},y_{n})+a_{n+1}\leq a_{n}-b_{n},$ (37)

with

$\displaystyle\left\\{\begin{array}[]{rcl}a_{n}&=&{\psi\over\psi-1}\|z_{n+1}-\bar{x}\|^{2}+{1\over\beta}\|y_{n-1}-\bar{y}\|^{2},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
b_{n}&=&\psi\|z_{n+1}-x_{n}\|^{2}+(1-\zeta)\big{(}\psi\|x_{n+1}-x_{n}\|^{2}+{1\over\beta}\|y_{n}-y_{n-1}\|^{2}\big{)}.\end{array}\right.$
(40)

Since $G(x_{n},y_{n})\geq 0$, $a_{n}\geq 0$ and $b_{n}\geq 0$, it follows from
Fact 2.2 that $\lim_{n\rightarrow\infty}a_{n}$ exists and
$\lim_{n\rightarrow\infty}b_{n}=0$. By the definition of $b_{n}$, this implies

$\lim\limits_{n\rightarrow\infty}\|z_{n+1}-x_{n}\|=\lim\limits_{n\rightarrow\infty}\|x_{n+1}-x_{n}\|=\lim\limits_{n\rightarrow\infty}\|y_{n+1}-y_{n}\|=0.$

Furthermore, we also have $\lim\limits_{n\rightarrow\infty}\|x_{n}-z_{n}\|=0$
since $x_{n}-z_{n}=\psi(x_{n}-z_{n+1})$.

It follows from (40) and the existence of $\lim_{n\rightarrow\infty}a_{n}$
that the sequences $\\{z_{n+1}\\}_{n\in{\mathbb{N}}}$ and
$\\{y_{n}\\}_{n\in{\mathbb{N}}}$, and thus
$\\{(z_{n+1},y_{n})\\}_{n\in{\mathbb{N}}}$, are bounded. Since
$\lim\limits_{n\rightarrow\infty}\|x_{n}-z_{n}\|=0$,
$\\{x_{n}\\}_{n\in{\mathbb{N}}}$ is also bounded. Let
$\\{(x_{n_{k}+1},y_{n_{k}})\\}_{k\in{\mathbb{N}}}$ be a subsequence of
$\\{(z_{n+1},y_{n})\\}_{n\in{\mathbb{N}}}$, which is convergent to
$(x^{*},y^{*})$. Then
$\lim_{n\rightarrow\infty}x_{n_{k}}=\lim_{n\rightarrow\infty}x_{n_{k}+1}=x^{*}$.
Similar to (31) and (32), for any
$(x,y)\in\mathbb{R}^{q}\times\mathbb{R}^{p}$, there hold

$\displaystyle\left\\{\begin{array}[]{l}\langle x_{n_{k}+1}-z_{n_{k}+1}+\tau
K^{\top}y_{n_{k}},~{}x-x_{n_{k}+1}\rangle\geq\tau\big{(}g(x_{n_{k}+1})-g(x)\big{)},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
\big{\langle}{1\over\beta}(y_{n_{k}}-y_{n_{k}-1})-\tau
Kx_{n_{k}},~{}y-y_{n_{k}}\big{\rangle}\geq\tau\big{(}f^{*}(y_{n_{k}})-f^{*}(y)\big{)}.\end{array}\right.$

By driving $k\rightarrow\infty$, taking into account that both $g$ and $f^{*}$
are closed (and thus lower semicontinuous) and cancelling out $\tau>0$, we
obtain

$\displaystyle\langle K^{\top}y^{*},~{}x-x^{*}\rangle\geq
g(x^{*})-g(x)\text{~{}~{}and~{}~{}}-\langle Kx^{*},~{}y-y^{*}\rangle\geq
f^{*}(y^{*})-f^{*}(y),$

which hold for any $(x,y)\in\mathbb{R}^{q}\times\mathbb{R}^{p}$. This implies
that $(x^{*},y^{*})\in{\mathcal{S}}$. Since all the discussions remain valid
for any $(\bar{x},\bar{y})\in{\mathcal{S}}$, $(\bar{x},\bar{y})$ can be
replaced by $(x^{*},y^{*})$ in the definition of
$\\{a_{n}\\}_{n\in{\mathbb{N}}}$ in the first place. As such, we have
$\lim_{k\rightarrow\infty}a_{n_{k}}=0$ since

$\displaystyle\lim_{k\rightarrow\infty}\|z_{n_{k}+1}-x^{*}\|=\lim_{k\rightarrow\infty}\|x_{n_{k}}-x^{*}\|=0\text{~{}~{}and~{}~{}}\lim_{k\rightarrow\infty}\|y_{n_{k}-1}-x^{*}\|=\lim_{k\rightarrow\infty}\|y_{n_{k}}-x^{*}\|=0.$

Since $\\{a_{n}\\}_{n\in{\mathbb{N}}}$ is monotonically nonincreasing, it
follows that $\lim_{n\rightarrow\infty}a_{n}=0$ and thus

$\displaystyle\lim_{n\rightarrow\infty}x_{n}=\lim_{n\rightarrow\infty}z_{n+1}=x^{*}\text{~{}~{}and~{}~{}}\lim_{n\rightarrow\infty}y_{n}=y^{*}.$

This completes the proof. $\Box$

###### Remark 3.1

Our analysis is motivated by Malitsky Malitsky2019Golden . In particular, to
get rid of the term $\|x_{n+1}-z_{n+1}\|^{2}$ in (36), it is required that
$\psi-1-{1\over\psi}\leq 0$, and apparently the golden ratio is an upper bound
of such $\psi$. As long as $\tau\sigma L^{2}<\psi$ is satisfied, larger $\psi$
allows larger step sizes $\tau$ and $\sigma$, which is helpful to improve
numerical performance. On the other hand, $\psi>1$ is necessary to guarantee
the nonnegativity of $\\{a_{n}\\}_{n\in{\mathbb{N}}}$.

###  3.4 Ergodic convergence rate

For convex-concave saddle point problems, many algorithms exhibit
${\mathcal{O}}(1/N)$ ergodic rate of convergence, see Chambolle2016ergodic ;
Monteiro2011Complexity ; Nemirovski2006Prox-Method . It is shown by Nemirovski
Nemirovski2006Prox-Method  that this rate is in fact tight. We next establish
the same ergodic rate of convergence for GRPDA, where the measure is the
primal-dual gap function defined in (19).

###### Theorem 3.2

Let $\\{(x_{n},y_{n},z_{n})\\}_{n\in{\mathbb{N}}}$ be the sequence generated
by Algorithm 3.1 from any initial point
$(x_{0},y_{0})\in\mathbb{R}^{q}\times\mathbb{R}^{p}$ and $z_{0}=x_{0}$, and
let $G(\cdot,\cdot)$ be defined as in (19) with any fixed
$(\bar{x},\bar{y})\in{\mathcal{S}}$. Let $N\geq 1$ be any integer and define

$\displaystyle
X_{N}=\frac{1}{N}\sum\nolimits_{n=1}^{N}x_{n}~{}~{}\mbox{and}~{}~{}Y_{N}=\frac{1}{N}\sum\nolimits_{n=1}^{N}y_{n}.$
(42)

Then, it holds that

$\displaystyle G(X_{N},Y_{N})\leq\frac{1}{2\tau
N}\Big{(}{\psi\over\psi-1}\|z_{2}-\bar{x}\|^{2}+{1\over\beta}\|y_{0}-\bar{y}\|^{2}\Big{)}.$

Proof. It follows from (37) and (40) that $2\tau G(x_{n},y_{n})\leq
a_{n}-a_{n+1}$ for all $n\geq 1$. Sum this up for $n=1,\ldots,N$, we obtain
$2\tau\sum_{n=1}^{N}G(x_{n},y_{n})\leq a_{1}-a_{N+1}\leq a_{1}$. Since
$G(x,y)=P(x)+D(y)$ and $P(x)$ and $D(y)$ are, respectively, convex in $x$ and
$y$, it follows from (42) that

$\displaystyle G(X_{N},Y_{N})=P(X_{N})+D(Y_{N})\leq{1\over
N}\sum_{n=1}^{N}P(x_{n})+{1\over N}\sum_{n=1}^{N}D(y_{n})={1\over
N}\sum_{n=1}^{N}G(x_{n},y_{n}).$

The conclusion follows since
$a_{1}={\psi\over\psi-1}\|z_{2}-\bar{x}\|^{2}+{1\over\beta}\|y_{0}-\bar{y}\|^{2}$.
$\Box$

##  4 Accelerated GRPDA

When either $g$ or $f^{*}$ is strongly convex, it was shown in Chambolle2011A
that one can adaptively choose the primal and the dual step sizes, as well as
the inertial parameter, so that PDA achieves a faster ${\mathcal{O}}(1/N^{2})$
convergence rate. Similar results have been achieved in Malitsky2018A  for PDA
with variable step sizes from linesearch. In this section, we show that the
same is true for GRPDA provided that the parameters are chosen adaptively.
Similar to Chambolle2011A ; Malitsky2018A , the ratio $\beta=\sigma/\tau$
plays an important role.

Under the regularity condition specified in Assumption 1.1, the “$\min$” and
the “$\max$” in the primal dual problem (2) can be switched
(Rockafellar1970Convex, , Corollary 31.2.1), and thus (2) is equivalent to

$\displaystyle\min_{y\in\mathbb{R}^{p}}\max_{x\in\mathbb{R}^{q}}~{}~{}f^{*}(y)-\langle
K^{\top}y,x\rangle-g(x),$ (43)

As such, by switching the roles of $(g,K,x,q)$ and $(f^{*},-K^{T},y,p)$, (43)
is reducible to (2). Therefore, we will only treat the case when $g$ is
strongly convex for succinctness, and it is completely analogues when $f^{*}$
is strongly convex. In this section, we thus make the following assumption.

###### Assumption 4.1

Assume that $g$ is $\gamma$-strongly convex, i.e., for some $\gamma>0$ it
holds that

$\displaystyle g(y)\geq g(x)+\langle
u,y-x\rangle+\frac{\gamma}{2}\|y-x\|^{2},~{}~{}\forall\,x,y\in\mathbb{R}^{q},\;\forall\,u\in\partial
g(x).$

Our accelerated GRPDA is presented below. Recall that $\phi$ denotes the
golden ratio.

######  Algorithm 4.1 (Accelerated GRPDA when $g$ is $\gamma$-strongly convex)

Step 0.

Let $\psi_{0}=1.3247...$ be the unique real root of $\psi^{3}-\psi-1=0$.
Choose $\psi\in(\psi_{0},\phi)$, $\beta_{0}>0$, $x_{0}\in\mathbb{R}^{q}$ and
$y_{0}\in\mathbb{R}^{p}$, and set $z_{0}=x_{0}$,
$\varphi=\frac{1+\psi}{\psi^{2}}$,
$\tau_{0}=\frac{1}{L}\sqrt{{\psi\over\beta_{0}}}$ and $n=1$.

Step 1.

Compute

$\displaystyle z_{n}$ $\displaystyle=$
$\displaystyle{\psi-1\over\psi}x_{n-1}+\frac{1}{\psi}z_{n-1},$ (44)
$\displaystyle x_{n}$ $\displaystyle=$ $\displaystyle{\rm
Prox}_{\tau_{n-1}g}(z_{n}-\tau_{n-1}K^{\top}y_{n-1}).$ (45)

Step 2.

Compute

$\displaystyle\omega_{n}$ $\displaystyle=$
$\displaystyle\frac{\psi-\varphi}{\psi+\varphi\gamma\tau_{n-1}},$ (46)
$\displaystyle\beta_{n}$ $\displaystyle=$
$\displaystyle\beta_{n-1}(1+\omega_{n}\gamma\tau_{n-1}),$ (47)
$\displaystyle\tau_{n}$ $\displaystyle=$
$\displaystyle\min\Big{\\{}\varphi\tau_{n-1},~{}~{}\frac{\psi}{\tau_{n-1}\beta_{n}L^{2}}\Big{\\}},$
(48) $\displaystyle y_{n}$ $\displaystyle=$ $\displaystyle{\rm
Prox}_{\beta_{n}\tau_{n}f^{*}}(y_{n-1}+\beta_{n}\tau_{n}Kx_{n}).$ (49)

Step 3.

Set $n\leftarrow n+1$ and return to Step 1.

###### Remark 4.1

We present the following remarks on Algorithm 4.1.

  1. 1.

In Algorithm 4.1, $\omega_{n}$ is used to update $\beta_{n}$, which plays a
key role in establishing the ${\mathcal{O}}(1/N^{2})$ convergence rate. The
condition $\psi>\psi_{0}$, where $\psi_{0}$ is the unique real root of
$\psi^{3}-\psi-1=0$, ensures that $\omega_{n}>0$.

  2. 2.

To achieve the ${\mathcal{O}}(1/N^{2})$ convergence rate, a balance between
the increase rate $\tau_{n}/\tau_{n-1}$ and the size of $\tau_{n}$ is
maintained in (48). In particular, the ratio
$\tau_{n}/\tau_{n-1}\leq\varphi\in(1,\psi)$ for $\psi\in(\psi_{0},\phi)$. It
follows from $\varphi=\frac{1+\psi}{\psi^{2}}$ that larger $\psi$ gives
smaller bound $\varphi$, and vice versa.

  3. 3.

Assume $\gamma=0$. Then $\beta_{n}=\beta_{0}$ for all $n\geq 1$. Furthremore,
it is easy to show by induction that, for all $n\geq 1$, the “$\min$” in (48)
is always attained by the second term on the right hand side and thus
$\tau_{n}=\tau_{0}$ for all $n\geq 1$. Therefore, Algorithm 4.1 reduces to
Algorithm 3.1, i.e., GRPDA with fixed step sizes $\tau=\tau_{0}$ and
$\sigma=\sigma_{0}:=\beta_{0}\tau_{0}$, in which case
$\tau_{0}\sigma_{0}L^{2}=\psi<\phi$ mets the requirement of Algorithm 3.1 as
we require $\psi<\phi$ in Algorithm 4.1.

We next establish some useful properties of $\tau_{n}$, $\beta_{n}$ and
$\omega_{n}$ generated in Algorithm 4.1. For convenience, we let
$\delta_{n}={\tau_{n}\over\tau_{n-1}}$ in the rest of this section.

###### Lemma 4.1

Let $\\{(\omega_{n},\beta_{n},\tau_{n})\\}_{n\in{\mathbb{N}}}$ be generated by
Algorithm 4.1 and define $\delta_{n}={\tau_{n}\over\tau_{n-1}}$ for $n\geq 1$.
Then, for all $n\geq 1$, there hold $0<\beta_{n-1}\leq\beta_{n}$,
$\frac{\psi-\varphi}{\psi+\varphi\gamma\sqrt{\varphi}\tau_{0}}\leq\omega_{n}<1$
and

$\displaystyle\frac{1}{\sqrt{\varphi(1+\gamma\sqrt{\varphi}\tau_{0})}}\frac{\sqrt{\psi}}{L\sqrt{\beta_{n}}}\leq\tau_{n}\leq\frac{\sqrt{\delta_{n}\psi}}{L\sqrt{\beta_{n}}}\leq\frac{\sqrt{\varphi\psi}}{L\sqrt{\beta_{n}}}\leq\frac{\sqrt{\varphi\psi}}{L\sqrt{\beta_{0}}}=\sqrt{\varphi}\tau_{0}.$
(50)

Moreover, there exists a constant $c>0$ such that $\beta_{n}\geq cn^{2}$ for
all $n\geq 1$.

Proof. Since $\beta_{0}>0$, $\omega_{n}>0$ for $n\geq 1$ and $\tau_{n}>0$ for
$n\geq 0$, it is clear from (47) that $0<\beta_{n-1}\leq\beta_{n}$ for all
$n\geq 1$. Recall that $\delta_{n}={\tau_{n}\over\tau_{n-1}}$. It follows from
(48) that
$\tau_{n}^{2}\beta_{n}\leq\frac{\delta_{n}\psi}{L^{2}}\leq\frac{\varphi\psi}{L^{2}}$.
Further considering $\beta_{n}\geq\beta_{0}>0$ for all $n\geq 1$ and
$\tau_{0}=\sqrt{\psi/\beta_{0}}/L$, we see that all the relations in (50),
except the first one, follow. Now, the relation
$\frac{\psi-\varphi}{\psi+\varphi\gamma\sqrt{\varphi}\tau_{0}}\leq\omega_{n}$
follows from (46) and $\tau_{n-1}\leq\sqrt{\varphi}\tau_{0}$ for all $n\geq
1$. On the other hand, $\omega_{n}<1$ is obvious. We next prove the first
inequality in (50) by induction.

It follows from $\tau_{0}=\frac{\sqrt{\psi}}{\sqrt{\beta_{0}}L}$ and
$\varphi>1$ that
$\tau_{0}\geq\frac{1}{\sqrt{\varphi(1+\gamma\sqrt{\varphi}\tau_{0})}}\frac{\sqrt{\psi}}{L\sqrt{\beta_{0}}}$.
Now, suppose that
$\tau_{k}\geq\frac{1}{\sqrt{\varphi(1+\gamma\sqrt{\varphi}\tau_{0})}}\frac{\sqrt{\psi}}{L\sqrt{\beta_{k}}}$
for $0\leq k\leq n$. Then, (48) implies that either

$\displaystyle\tau_{n+1}=\varphi\tau_{n}\geq\tau_{n}\geq\frac{1}{\sqrt{\varphi(1+\gamma\sqrt{\varphi}\tau_{0})}}\frac{\sqrt{\psi}}{L\sqrt{\beta_{n}}}\geq\frac{1}{\sqrt{\varphi(1+\gamma\sqrt{\varphi}\tau_{0})}}\frac{\sqrt{\psi}}{L\sqrt{\beta_{n+1}}}$

since $\beta_{n+1}\geq\beta_{n}$ and $\varphi>1$, or

$\displaystyle\tau_{n+1}=\frac{\psi}{\tau_{n}\beta_{n+1}L^{2}}$
$\displaystyle\geq$
$\displaystyle\frac{L\sqrt{\beta_{n}}}{\sqrt{\varphi\psi}}\frac{\psi}{\beta_{n+1}L^{2}}=\sqrt{\frac{\beta_{n}}{\varphi\beta_{n+1}}}\frac{\sqrt{\psi}}{\sqrt{\beta_{n+1}}L}$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{{\varphi(1+\omega_{n+1}\gamma\tau_{n})}}}\frac{\sqrt{\psi}}{\sqrt{\beta_{n+1}}L}$
$\displaystyle\geq$
$\displaystyle\frac{1}{\sqrt{\varphi(1+\gamma\sqrt{\varphi}\tau_{0})}}\frac{\sqrt{\psi}}{\sqrt{\beta_{n+1}}L},$

where the first inequality follows from
$\tau_{n}\leq\frac{\sqrt{\varphi\psi}}{L\sqrt{\beta_{n}}}$ and the second is
due to $\omega_{n+1}<1$ and $\tau_{n}\leq\sqrt{\varphi}\tau_{0}$ for all
$n\geq 0$. Therefore, the first inequality in (50) also holds for all $n\geq
1$.

To complete the proof, it only remains to show that there exists $c>0$ such
that $\beta_{n}\geq cn^{2}$ for all $n\geq 1$. For simplicity, we let
$\omega:=\frac{\psi-\varphi}{\psi+\varphi\gamma\sqrt{\varphi}\tau_{0}}$. Then,
$\omega_{n}\geq\omega$ for all $n\geq 1$ and thus

$\displaystyle\beta_{n+1}=\beta_{n}(1+\omega_{n+1}\gamma\tau_{n})$
$\displaystyle\geq$
$\displaystyle\beta_{n}\Big{(}1+\frac{\omega\gamma}{\sqrt{\varphi(1+\gamma\sqrt{\varphi}\tau_{0})}}\frac{\sqrt{\psi}}{L\sqrt{\beta_{n}}}\Big{)}=\beta_{n}+c_{0}\sqrt{\beta_{n}},$
(51)

where
$c_{0}:=\frac{\omega\gamma\sqrt{\psi}}{\sqrt{\varphi(1+\gamma\sqrt{\varphi}\tau_{0})}L}>0$.
By induction, it is easy to show that $\beta_{n}\geq cn^{2}$ for all $n\geq 1$
with $c:=\min(c_{0}^{2}/9,\beta_{1})>0$. This completes the proof. $\Box$

Now, we are ready to establish the promised ${\mathcal{O}}(1/N^{2})$ ergodic
convergence rate.

###### Theorem 4.1

Let
$\\{(z_{n},x_{n},y_{n},\omega_{n},\beta_{n},\tau_{n})\\}_{n\in{\mathbb{N}}}$
be the sequence generated by Algorithm 4.1 and
$(\bar{x},\bar{y})\in{\mathcal{S}}$ be any solution of (2). For any integer
$N\geq 1$, define

$\displaystyle
S_{N}=\sum_{n=1}^{N}\beta_{n}\tau_{n},~{}~{}X_{N}=\frac{1}{S_{N}}\sum_{n=1}^{N}\beta_{n}\tau_{n}x_{n}\text{~{}~{}and~{}~{}}Y_{N}=\frac{1}{S_{N}}\sum_{n=1}^{N}\beta_{n}\tau_{n}y_{n}.$

Then, there hold $\|z_{N+2}-\bar{x}\|={\mathcal{O}}(1/N)$ and
$G(X_{N},Y_{N})={\mathcal{O}}(1/N^{2})$.

Proof. Fix $n\geq 1$. Since $g$ is $\gamma$-strongly convex, an inequality
stronger than the one stated in Fact 2.1 can be used. It follows from (45)
that $z_{n}-\tau_{n-1}K^{\top}y_{n-1}-x_{n}\in\tau_{n-1}\partial g(x_{n})$,
and thus

$\displaystyle\langle
x_{n}-z_{n}+\tau_{n-1}K^{\top}y_{n-1},~{}x-x_{n}\rangle\geq\tau_{n-1}\Big{(}g(x_{n})-g(x)+\frac{\gamma}{2}\|x_{n}-x\|^{2}\Big{)},\;\forall\,x.$
(52)

By passing $n+1$ to $n$ and $\bar{x}$ to $x$ in (52), we obtain

$\displaystyle\langle
x_{n+1}-z_{n+1}+\tau_{n}K^{\top}y_{n},~{}~{}\bar{x}-x_{n+1}\rangle\geq\tau_{n}\Big{(}g(x_{n+1})-g(\bar{x})+\frac{\gamma}{2}\|x_{n+1}-\bar{x}\|^{2}\Big{)}.$
(53)

Similarly, by passing $x_{n+1}$ to $x$ in (52) and multiplying both sides by
$\delta_{n}={\tau_{n}\over\tau_{n-1}}$, we obtain

$\displaystyle\langle\delta_{n}(x_{n}-z_{n})+\tau_{n}K^{\top}y_{n-1},~{}x_{n+1}-x_{n}\rangle$
$\displaystyle\geq$
$\displaystyle\tau_{n}\Big{(}g(x_{n})-g(x_{n+1})+\frac{\gamma}{2}\|x_{n+1}-x_{n}\|^{2}\Big{)}.$
(54)

Similar to (32), it follows from (49) that

$\displaystyle\Big{\langle}{1\over\beta_{n}}(y_{n}-y_{n-1})-\tau_{n}Kx_{n},~{}\bar{y}-y_{n}\Big{\rangle}\geq\tau_{n}\big{(}f^{*}(y_{n})-f^{*}(\bar{y})\big{)}$
(55)

From (44), it is easy to derive $x_{n}-z_{n}=\psi(x_{n}-z_{n+1})$. Then, by
adding (53)-(55) and using similar arguments as in Lemma 3.1, we obtain

$\displaystyle\tau_{n}G(x_{n},y_{n})$ $\displaystyle\leq$
$\displaystyle\langle
x_{n+1}-z_{n+1},\bar{x}-x_{n+1}\rangle+\frac{1}{\beta_{n}}\big{\langle}y_{n}-y_{n-1},\bar{y}-y_{n}\big{\rangle}+\psi\delta_{n}\big{\langle}x_{n}-z_{n+1},x_{n+1}-x_{n}\big{\rangle}$
(56) $\displaystyle+\tau_{n}\langle
K^{\top}(y_{n}-y_{n-1}),x_{n}-x_{n+1}\rangle-\gamma\tau_{n}\|x_{n+1}-\bar{x}\|^{2}.$

By using (16) and Cauchy-Schwartz inequality and combining terms, (56) implies

$\displaystyle(1+2\gamma\tau_{n})\|x_{n+1}-\bar{x}\|^{2}+\frac{1}{\beta_{n}}\|y_{n}-\bar{y}\|^{2}+2\tau_{n}G(x_{n},y_{n})$
$\displaystyle\leq$
$\displaystyle\|z_{n+1}-\bar{x}\|^{2}+\frac{1}{\beta_{n}}\|y_{n-1}-\bar{y}\|^{2}-\psi\delta_{n}\|z_{n+1}-x_{n}\|^{2}+(\psi\delta_{n}-1)\|x_{n+1}-z_{n+1}\|^{2}$
$\displaystyle-\psi\delta_{n}\|x_{n+1}-x_{n}\|^{2}-\frac{1}{\beta_{n}}\|y_{n}-y_{n-1}\|^{2}+2\tau_{n}\|K^{\top}(y_{n}-y_{n-1})\|\|x_{n+1}-x_{n}\|.$

Plugging in (35) and recalling $L=\|K\|=\|K^{\top}\|$, we obtain

$\displaystyle(1+2\gamma\tau_{n})\frac{\psi}{\psi-1}\|z_{n+2}-\bar{x}\|^{2}+\frac{1}{\beta_{n}}\|y_{n}-\bar{y}\|^{2}+2\tau_{n}G(x_{n},y_{n})$
(57) $\displaystyle\leq$
$\displaystyle\frac{\psi+2\gamma\tau_{n}}{\psi-1}\|z_{n+1}-\bar{x}\|^{2}+\frac{1}{\beta_{n}}\|y_{n-1}-\bar{y}\|^{2}+\Big{(}\psi\delta_{n}-1-\frac{1+2\gamma\tau_{n}}{\psi}\Big{)}\|x_{n+1}-z_{n+1}\|^{2}$
$\displaystyle-\psi\delta_{n}\|z_{n+1}-x_{n}\|^{2}-\psi\delta_{n}\|x_{n+1}-x_{n}\|^{2}-\frac{1}{\beta_{n}}\|y_{n}-y_{n-1}\|^{2}$
$\displaystyle+2\tau_{n}L\|y_{n}-y_{n-1}\|\|x_{n+1}-x_{n}\|.$

It follows from $\tau_{n}\leq\frac{\sqrt{\psi\delta_{n}}}{L\sqrt{\beta_{n}}}$
(the second inequality in (50)) that

$\displaystyle 2\tau_{n}L\|y_{n}-y_{n-1}\|\|x_{n+1}-x_{n}\|$
$\displaystyle\leq$
$\displaystyle\psi\delta_{n}\|x_{n+1}-x_{n}\|^{2}+\frac{1}{\beta_{n}}\|y_{n}-y_{n-1}\|^{2}.$

Furthermore,
$\psi\delta_{n}-1-\frac{1+2\gamma\tau_{n}}{\psi}\leq\psi\varphi-1-\frac{1}{\psi}=0$
since $\delta_{n}\leq\varphi$. Therefore, (57) implies

$\displaystyle(1+{2}\gamma\tau_{n})\frac{\psi}{\psi-1}\|z_{n+2}-\bar{x}\|^{2}+\frac{1}{\beta_{n}}\|y_{n}-\bar{y}\|^{2}+2\tau_{n}G(x_{n},y_{n})$
(58) $\displaystyle\leq$
$\displaystyle\frac{\psi+{2}\gamma\tau_{n}}{\psi-1}\|z_{n+1}-\bar{x}\|^{2}+\frac{1}{\beta_{n}}\|y_{n-1}-\bar{y}\|^{2}.$

Since
$(1+{2}\gamma\tau_{n})\frac{\psi}{\psi-1}=\frac{\psi(1+{2}\gamma\tau_{n})}{\psi+{2}\gamma\tau_{n+1}}\frac{\psi+{2}\gamma\tau_{n+1}}{\psi-1}$
and

$\displaystyle\frac{\psi(1+{2}\gamma\tau_{n})}{\psi+{2}\gamma\tau_{n+1}}$
$\displaystyle\geq$
$\displaystyle\frac{\psi(1+{2}\gamma\tau_{n})}{\psi+{2}\gamma\varphi\tau_{n}}=1+\frac{\psi-\varphi}{{\psi\over
2}+\gamma\varphi\tau_{n}}\gamma\tau_{n}\geq
1+\frac{\psi-\varphi}{\psi+\gamma\varphi\tau_{n}}\gamma\tau_{n}=1+\omega_{n+1}\gamma\tau_{n},$

where the first inequality follows from $\tau_{n+1}\leq\varphi\tau_{n}$ and
the second is due to $\psi^{3}-\psi-1>0$, it follows that

$\displaystyle(1+{2}\gamma\tau_{n})\frac{\psi}{\psi-1}\geq(1+\omega_{n+1}\gamma\tau_{n})\frac{\psi+2\gamma\tau_{n+1}}{\psi-1}=\frac{\beta_{n+1}}{\beta_{n}}\frac{\psi+2\gamma\tau_{n+1}}{\psi-1}.$
(59)

Define
$A_{n}:=\frac{\psi+2\gamma\tau_{n}}{2(\psi-1)}\|z_{n+1}-\bar{x}\|^{2}+\frac{1}{2\beta_{n}}\|y_{n-1}-\bar{y}\|^{2}$.
Combining (58) and (59), we deduce

$\displaystyle\beta_{n+1}A_{n+1}+\beta_{n}\tau_{n}G(x_{n},y_{n})\leq\beta_{n}A_{n}.$
(60)

By summing (60) for $n=1,\ldots,N$, we obtain

$\displaystyle\beta_{N+1}A_{N+1}+\sum_{n=1}^{N}\beta_{n}\tau_{n}G(x_{n},y_{n})\leq\beta_{1}A_{1}.$
(61)

The convexity of $G(x,y)$, (61) and the definition of $A_{n}$ imply that

$\displaystyle G(X_{N},Y_{N})$ $\displaystyle\leq$ $\displaystyle{1\over
S_{N}}\sum_{n=1}^{N}\beta_{n}\tau_{n}G(x_{n},y_{n})\leq{\beta_{1}A_{1}\over
S_{N}},$ (62) $\displaystyle\|z_{N+2}-\bar{x}\|^{2}$ $\displaystyle\leq$
$\displaystyle\frac{2(\psi-1)}{\psi+2\gamma\tau_{N+1}}{\beta_{1}A_{1}\over\beta_{N+1}}\leq{2\beta_{1}A_{1}\over\beta_{N+1}}.$
(63)

From Lemma 4.1, there exists $c>0$ such that $\beta_{n}\geq cn^{2}$ for all
$n\geq 1$. Then, (63) gives

$\displaystyle\|z_{N+2}-\bar{x}\|$ $\displaystyle\leq$
$\displaystyle\frac{\sqrt{2\beta_{1}A_{1}/c}}{N+1}.$

It follows from (51) that $\beta_{n+1}-\beta_{n}\geq c_{0}\sqrt{c}n$, which
together with (47) and $\omega_{n+1}<1$ implies

$\beta_{n}\tau_{n}=\frac{\beta_{n+1}-\beta_{n}}{\omega_{n+1}\gamma}\geq\frac{c_{0}\sqrt{c}n}{\gamma}.$

As a result, $S_{N}=\sum_{n=1}^{N}\beta_{n}\tau_{n}={\mathcal{O}}(N^{2})$, and
thus $G(X_{N},Y_{N})={\mathcal{O}}(1/N^{2})$ follows from (62). $\Box$

##  5 GRPDA for two special cases

This section is devoted to GRPDA for two special cases of the structured
convex optimization problem (1), which are specified in the following
assumption.

###### Assumption 5.1

Let $b\in\mathbb{R}^{p}$ be given. Assume that either
$f(\cdot)=\iota_{\\{b\\}}(\cdot)$, the indicator function of the singleton
$\\{b\\}$, or $f(\cdot)={1\over 2}\|\cdot-b\|^{2}$.

The two cases specified in Assumption 5.1 correspond respectively to

$\displaystyle\min_{x\in\mathbb{R}^{q}}\\{g(x):Kx=b\\}\text{~{}~{}and~{}~{}}\min_{x\in\mathbb{R}^{q}}\frac{1}{2}\|Kx-b\|^{2}+g(x),$

which are abundant in practice. For example, linear inverse problems, which
include many signal and image processing applications, usually enforce data
fitting via $f(Kx)$, which takes the form $\frac{1}{2}\|Kx-b\|^{2}$ for noisy
data and $\iota_{\\{b\\}}(Kx)$, or equivalently $Kx=b$, for ideal noiseless
data. A well known application is compressive signal/image sensing
Dono2006Compressed ; Needell2013compressed , which recovers sparse or
compressible signals from a small number of linear measurements via solving
the basis pursuit and/or the LASSO problems Tibshirani1996Regression ;
Chen2001Atomic .

The aim of this section is to propose a relaxed GRPDA when $f$ satisfies
Assumption 5.1, which is a modification of Algorithm 3.1 in two aspects: (i)
extending the permitted range of $\psi$ from $(1,\phi]$ to $(1,2]$, and (ii)
introducing a relaxation step. In the rest of this section, we always assume
that $f$ satisfies Assumption 5.1. Then, for any $u\in\mathbb{R}^{p}$ and
$\sigma>0$, it is easy to show that ${\rm Prox}_{\sigma f^{*}}(u)=\eta
u+\varrho b$, where

$\displaystyle\begin{array}[]{l}(\eta,\varrho):=(\eta(\sigma),\varrho(\sigma))=\left\\{\begin{array}[]{ll}(1,-\sigma)&\hbox{\quad
if $f(\cdot)=\iota_{\\{b\\}}(\cdot)$,}\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
{(1,\sigma)\over 1+\sigma}&\hbox{\quad if $f(\cdot)={1\over
2}\|\cdot-b\|^{2}$.}\end{array}\right.\end{array}$ (67)

In both cases, it holds that $\eta=\eta(\sigma)\in(0,1]$ for any $\sigma>0$.

###  5.1 GRPDA as fixed point iteration

The iterative scheme GRPDA (23) applied to (1) appears as

$\displaystyle\left\\{\begin{array}[]{rcl}y_{n-1}&=&\eta(y_{n-2}+\sigma
Kx_{n-1})+\varrho b,\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
z_{n}&=&\frac{\psi-1}{\psi}x_{n-1}+{1\over\psi}z_{n-1},\vskip 3.0pt plus 1.0pt
minus 1.0pt\\\ x_{n}&=&{\rm Prox}_{\tau g}(z_{n}-\tau
K^{\top}y_{n-1}),\end{array}\right.$ (71)

where $\eta$ and $\varrho$ are determined by (67). Note that, since $f$
satisfies Assumption 5.1, the schemes (71) and (23) are equivalent in the
sense that they, if initialized properly, generalize exactly the same sequence
of iterates. Our choice of (71), which states explicitly the computing
formulas of $(y_{n-1},z_{n},x_{n})$, instead of $(z_{n},x_{n},y_{n})$ as in
(23), is mainly for convenience of analysis.

To present our relaxed GRPDA, we first represent (71) as a fixed point
iterative scheme, which is summarized below.

###### Lemma 5.1

Let $\xi_{n}=(z_{n};\,x_{n};\,y_{n-1})$. Then, the iterative scheme (71) can
be represented as the fixed point iteration
$\xi_{n}=\mathbb{S}\circ\mathbb{T}(\xi_{n-1})$, where

$\displaystyle\left\\{\begin{array}[]{ccl}\mathbb{T}(\cdot)=T(\cdot)+\vartheta_{1},&&\vartheta_{1}=(0;\,-\varrho\tau
K^{\top}b;\,0),\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\mathbb{S}(\cdot)=S(\cdot)+\vartheta_{2},&&\vartheta_{2}=(0;\,0^{\top};\,\varrho
b),\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
S=\left[\begin{array}[]{ccccc}I&&&0&0\\\ 0&&&{\rm Prox}_{\tau g}&0\\\
0&&&0&I\end{array}\right]&\text{~{}~{}and~{}~{}}&T=\left[\begin{array}[]{ccccc}\frac{1}{\psi}I&{}{}\hfil&(1-\frac{1}{\psi})I&{}{}\hfil&0\\\
\frac{1}{\psi}I&{}{}\hfil&(1-\frac{1}{\psi})I-\eta\tau\sigma
K^{\top}K&{}{}\hfil&-\eta\tau K^{\top}\\\ 0&{}{}\hfil&\eta\sigma
K&{}{}\hfil&\eta I\end{array}\right].\end{array}\right.$ (81)

###### Proof

Direct verification from (71).

Let ${\mathbb{C}}$ be the set of complex numbers and $\Lambda(M)$ be the set
of eigenvalues of a matrix $M$. The following theorem is our key to extend
$\psi$ from $(1,\phi]$ to $(1,2]$ and to take a relaxation step.

###### Theorem 5.1

Let $\tau,\sigma>0$ and $\psi\in(1,2]$ be such that $\tau\sigma L^{2}<\psi$
and $\eta\in(0,1]$. Then, the matrix $T$, and thus the affine mapping
$\mathbb{T}$, defined in (81) is firmly nonexpansive.

Proof. It follows from, e.g., (Bauschke2011book, , Proposition 4.2) that $T$
is firmly nonexpansive if and only if $2T-I$ is nonexpansive, which is clearly
equivalent to $|\lambda|\leq 1$ for any $\lambda\in\Lambda(2T-I)$. Since
$\eta\in(0,1]$, we have $2\eta-1\in(-1,1]$. Therefore, to guarantee
nonexpansiveness of $2T-I$, it is sufficient to show that $|\lambda|\leq 1$
for any $\lambda\in\Lambda(2T-I)\setminus\\{\pm 1,2\eta-1\\}$. The reason that
we exclude $\pm 1$ and $2\eta-1$ from the spectrum $\Lambda(2T-I)$ will be
clear below.

Let $\lambda\in\Lambda(2T-I)\setminus\\{\pm 1,2\eta-1\\}$ and
$\xi=(z;\,x;\,y)\in\mathbb{R}^{q}\times\mathbb{R}^{q}\times\mathbb{R}^{p}$ be
any eigen-pair of $2T-I$, i.e., $(2T-I)\xi=\lambda\xi$ and
$\lambda\notin\\{\pm 1,2\eta-1\\}$. By the definition of $T$,
$(2T-I)\xi=\lambda\xi$ appears as

$\displaystyle\Big{(}\frac{2}{\psi}-1\Big{)}z+\Big{(}2-\frac{2}{\psi}\Big{)}x=\lambda
z,$ (82)
$\displaystyle\frac{2}{\psi}z+\Big{(}1-\frac{2}{\psi}\Big{)}x-2\eta\tau\sigma
K^{\top}Kx-2\eta\tau K^{\top}y=\lambda x,$ (83) $\displaystyle 2\eta\sigma
Kx+(2\eta-1)y=\lambda y.$ (84)

Combining (82) and (83), we obtain

$\displaystyle(\lambda+1)(z-x)-2\eta\tau\sigma K^{\top}Kx-2\eta\tau
K^{\top}y=0.$ (85)

First, we show that $KK^{\top}y\neq 0$ by contradiction. Assume that
$KK^{\top}y=0$. Then, it holds that $KK^{\top}Kx=0$ from (84), combining which
with (85) gives $(\lambda+1)K(z-x)=0$. Since $\lambda\neq-1$, there must hold
$Kz=Kx$. Then, multiplying both sides of (83) by $K$ shows that
$(1-\lambda)Kx=0$, which further implies that $Kx=0$ since $\lambda\neq 1$. As
a result, it follows from (84) and $\lambda\neq 2\eta-1$ that there must hold
$y=0$. Thus, along with $Kx=0$ and $\lambda\neq-1$, (85) would lead to $x=z$.
Since $\lambda\neq 1$, it then follows from (82) that $x=z=0$, which
contradicts to the fact that $\xi$ is an eigenvector and has to be nonzero.
Therefore, there must hold $KK^{\top}y\neq 0$.

Since $KK^{\top}y\neq 0$, we thus have $y\neq 0$. For simplicity, we let
$\lambda_{1}:=\frac{\lambda+1}{2}$. Since $\lambda\notin\\{\pm 1,2\eta-1\\}$,
we have $\lambda_{1}\notin\\{0,1,\eta\\}$. Then, (82) and (84) can be
restated, respectively, as

$\displaystyle\Big{(}1-\frac{1}{\psi}\Big{)}x=\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}z\text{~{}~{}and~{}~{}}\sigma
Kx=\frac{\lambda_{1}-\eta}{\eta}y.$ (86)

Multiplying both sides of (84) by $\tau K^{\top}$, adding to (83) and using
(82), we obtain

$(\lambda+1)(z-x-\tau K^{\top}y)=0.$

Since $\lambda\neq-1$, we thus obtain $z-x=\tau K^{\top}y$, wihch together
with (86) gives

$\displaystyle-\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}\frac{\lambda_{1}-\eta}{\eta}y$
$\displaystyle=$ $\displaystyle-\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}\sigma
Kx$ $\displaystyle=$
$\displaystyle\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}\sigma K(\tau
K^{\top}y-z)$ $\displaystyle=$
$\displaystyle\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}\tau\sigma
KK^{\top}y-\sigma\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}Kz$ $\displaystyle=$
$\displaystyle\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}\tau\sigma
KK^{\top}y-\sigma\Big{(}1-\frac{1}{\psi}\Big{)}Kx$ $\displaystyle=$
$\displaystyle\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}\tau\sigma
KK^{\top}y-\Big{(}1-\frac{1}{\psi}\Big{)}\frac{\lambda_{1}-\eta}{\eta}y.$

By combining the two terms of $y$ on both sides, we obtain

$\displaystyle(\lambda_{1}-1)\frac{\eta-\lambda_{1}}{\eta}y$ $\displaystyle=$
$\displaystyle\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}\tau\sigma KK^{\top}y,$

which simplifies to

$\displaystyle\frac{\eta-\lambda_{1}}{\eta}y=\frac{(\overline{\lambda_{1}}-1)\Big{(}\lambda_{1}-\frac{1}{\psi}\Big{)}}{|\lambda_{1}-1|^{2}}\tau\sigma
KK^{\top}y.$ (87)

Here $\overline{\lambda_{1}}$ denotes the complex conjugate of $\lambda_{1}$.
Recall that $y\neq 0$. Then (87) implies that $y$ is an eigenvector of
$KK^{\top}$. Since $KK^{\top}$ is real symmetric, positive semidefinite and
$\|KK^{\top}\|=L^{2}$, its eigenvalues must be real and lie in $[0,L^{2}]$. It
is then implied by $KK^{T}y\neq 0$ and $0<\tau\sigma L^{2}<\psi$ that there
exists $\psi_{1}\in(0,\psi)$ such that $\tau\sigma KK^{\top}y=\psi_{1}y$.

Let $s,t\in\mathbb{R}$ and $\lambda=s+t\sqrt{-1}\in{\mathbb{C}}$. Then
$\lambda_{1}=\frac{1+s}{2}+\frac{t}{2}\sqrt{-1}$ and it is elementary to
deduce from (87) that

$\displaystyle\frac{4\eta-2-2s}{\eta}$ $\displaystyle=$
$\displaystyle\frac{(s+1-\frac{2}{\psi})(s-1)+t^{2}}{|\lambda_{1}-1|^{2}}\psi_{1},\vskip
3.0pt plus 1.0pt minus 1.0pt$ (88) $\displaystyle\frac{t}{\eta}$
$\displaystyle=$
$\displaystyle\frac{t\Big{(}1-\frac{1}{\psi}\Big{)}\psi_{1}}{|\lambda_{1}-1|^{2}}.$
(89)

We split the discussions in two cases, (i) $t=0$, and (ii) $t\neq 0$.

  * •

Case (i). If $t=0$, we have $s^{2}-2\eta(1-\psi_{1})s-1+2\eta\varsigma=0$ from
(88), where $\varsigma:=1+\psi_{1}-\frac{2\psi_{1}}{\psi}$. Let
$u:=\eta(1-\psi_{1})$ and $v:=\eta^{2}(1-\psi_{1})^{2}-2\eta\varsigma+1$.
Then, we have $\lambda=s=u\pm\sqrt{v}$. Since $s\in{\mathbb{R}}$, there must
hold $\eta^{2}(1-\psi_{1})^{2}-2\eta\varsigma+1\geq 0$, which, by further
considering $\eta\leq 1$, implies that

$\displaystyle
0<\eta\leq\left\\{\begin{array}[]{ll}\min\Big{\\{}1,~{}~{}\frac{\varsigma-\sqrt{\varsigma^{2}-(1-\psi_{1})^{2}}}{(1-\psi_{1})^{2}}\Big{\\}},&~{}~{}\psi_{1}\neq
1,\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\min\Big{\\{}1,~{}~{}\frac{1}{4(1-\frac{1}{\psi})}\Big{\\}},&~{}~{}\psi_{1}=1.\end{array}\right.$

By $\psi_{1}\in(0,\psi)$ and $\psi>1$, it is elementary to show that
$\sqrt{v}<1-u$ and $\sqrt{v}<1+u$. Consequently, we have $\lambda\in(-1,1)$.

  * •

Case (ii). If $t\neq 0$, we have
$|\lambda_{1}-1|^{2}=\eta\psi_{1}\Big{(}1-\frac{1}{\psi}\Big{)}$ from (89).
Then, (88) implies

$(4\eta-2-2s)\Big{(}1-{1\over\psi}\Big{)}=\Big{(}s+1-{2\over\psi}\Big{)}(s-1)+t^{2},$

from which we obtain

$\displaystyle\Big{(}s+1-\frac{2}{\psi}\Big{)}^{2}+t^{2}$ $\displaystyle=$
$\displaystyle(4\eta-2-2s)\Big{(}1-{1\over\psi}\Big{)}+\Big{(}s+1-{2\over\psi}\Big{)}\Big{(}2-{2\over\psi}\Big{)}$
(91) $\displaystyle=$ $\displaystyle
4\Big{(}1-{1\over\psi}\Big{)}\Big{(}\eta-{1\over\psi}\Big{)}$ $\displaystyle=$
$\displaystyle\Big{(}1-\frac{2}{\psi}\Big{)}^{2}+4\eta\Big{(}1-\frac{1}{\psi}\Big{)}-1$
$\displaystyle\leq$
$\displaystyle\Big{(}1-\frac{2}{\psi}\Big{)}^{2}+4\Big{(}1-\frac{1}{\psi}\Big{)}-1$
$\displaystyle=$ $\displaystyle 4\Big{(}1-\frac{1}{\psi}\Big{)}^{2},$

where the inequality follows from $\eta\in(0,1]$ and $\psi>1$. Since
$\psi\in(1,2]$, a simple geometric argument based on (91) shows that
$|\lambda|=\sqrt{s^{2}+t^{2}}\leq 1$.

Combining Cases (i) and (ii), we have shown that $|\lambda|\leq 1$ for any
$\tau,\sigma>0$ and $\psi\in(1,2]$ such that $\tau\sigma L^{2}<\psi$. This
completes the proof. $\Box$

###  5.2 Relaxed GRPDA

Now, we are ready to present a relaxed GRPDA when $f$ satisfies Assumption
5.1. The relaxation parameters lie in $(0,2/3)$ and the parameter $\psi$,
which is restricted to $(1,\phi]$ in Algorithm 3.1, can now be extended to
$(1,2]$. The relaxed GRPDA is summarized below.

######  Algorithm 5.1 (Relaxed GRPDA)

Step 0.

Let $\tau,\sigma>0$ and $\psi\in(1,2]$ be such that $\tau\sigma L^{2}<\psi$,
$(\eta,\varrho)=(\eta(\sigma),\varrho(\sigma))$ be defined in (67), and
$\\{\rho_{n}\\}_{n\in{\mathbb{N}}}\subseteq(0,\frac{3}{2})$ be such that
$\sum_{n\in\mathbb{N}}\rho_{n}(1-\frac{2\rho_{n}}{3})=+\infty$. Choose
$x_{0}\in\mathbb{R}^{q}$ and $y_{-1}\in\mathbb{R}^{p}$. Set $z_{0}=x_{0}$ and
$n=1$.

Step 1.

Compute

$\displaystyle\left\\{\begin{array}[]{rcl}\widetilde{y}_{n-1}&=&\eta(y_{n-2}+\sigma
Kx_{n-1})+\varrho b,\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
\widetilde{z}_{n}&=&\frac{\psi-1}{\psi}x_{n-1}+{1\over\psi}z_{n-1},\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ \widetilde{x}_{n}&=&{\rm Prox}_{\tau
g}(\widetilde{z}_{n}-\tau K^{\top}\widetilde{y}_{n-1}).\end{array}\right.$

Step 2.

Take a relaxation step

$\displaystyle\left\\{\begin{array}[]{rcl}y_{n-1}&=&y_{n-2}+\rho_{n}(\widetilde{y}_{n-1}-y_{n-2}),\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
z_{n}&=&z_{n-1}+\rho_{n}(\widetilde{z}_{n}-z_{n-1}),\vskip 3.0pt plus 1.0pt
minus 1.0pt\\\
x_{n}&=&x_{n-1}+\rho_{n}(\widetilde{x}_{n}-x_{n-1}).\end{array}\right.$

Step 3.

Set $n\leftarrow n+1$ and return to Step 1.

The convergence of Algorithm 5.1 is established in the following theorem. The
key of its proof is to observe, based on Theorem 5.1 and (Bauschke2011book, ,
Proposition 4.32), that the operator $\mathbb{S}\circ\mathbb{T}$ is
$2/3$-averaged111An operator $P$ is $\alpha$-averaged for some
$\alpha\in(0,1)$ if there exists a nonexpansive operator $Q$ such that
$P=(1-\alpha)I+\alpha Q$., and thus the convergence result follows from
(Bauschke2011book, , Proposition 5.15).

###### Theorem 5.2

Let $f$ satisfy Assumption 5.1 and
$\\{(z_{n},x_{n},y_{n-1})\\}_{n\in{\mathbb{N}}}$ be the sequence generated by
Algorithm 5.1 from any initial point
$(x_{0},y_{-1})\in\mathbb{R}^{q}\times\mathbb{R}^{p}$ and $z_{0}=x_{0}$. Then,
$\\{(x_{n},y_{n})\\}_{n\in{\mathbb{N}}}$ converges to a solution of (2), i.e.,
an element in ${\mathcal{S}}$.

###### Proof

Let $\mathbb{S}$ and $\mathbb{T}$ be defined as in (81),
$\mathbb{G}:=\mathbb{S}\circ\mathbb{T}$, and $\xi_{n}=(z_{n};x_{n};y_{n-1})$
for $n\geq 0$. Then, the sequence
$\\{(z_{n},x_{n},y_{n-1})\\}_{n\in{\mathbb{N}}}$ generated by Algorithm 5.1
from $(z_{0},x_{0},y_{-1})$ satisfies
$\xi_{n}=\xi_{n-1}+\rho_{n}(\mathbb{G}(\xi_{n-1})-\xi_{n-1})$ for $n\geq 1$.
The key of the proof is to show that $\mathbb{G}$ is $2/3$-averaged, which we
argue below.

Apparently, $\eta$ defined in (67) lies in $(0,1]$. It then follows from
Theorem 5.1 that $\mathbb{T}$ defined in (81) is firmly nonexpansive under the
condition of Algorithm 5.1, i.e., $\tau,\sigma>0$ and $\tau\sigma
L^{2}<\psi\in(1,2]$. On the other hand, it is well known that the proximal
operator ${\rm Prox}_{\tau g}$ is firmly nonexpansive for any $\tau>0$. Then,
by the definition in (81), $\mathbb{S}$ is also firmly nonexpansive.
Consequently, it follows from (Bauschke2011book, , Proposition 4.32) that
$\mathbb{G}=\mathbb{S}\circ\mathbb{T}$ is $2/3$-averaged.

Denote by $\mbox{Fix}(\mathbb{G})$ the set of fixed points of $\mathbb{G}$,
i.e.,
$\mbox{Fix}(\mathbb{G}):=\left\\{\xi\in\mathbb{R}^{q}\times\mathbb{R}^{q}\times\mathbb{R}^{p}:\xi=\mathbb{G}(\xi)\right\\}$.
By the construction of $\mathbb{S}$ and $\mathbb{T}$, it is elementary to show
that

$\mbox{Fix}(\mathbb{G})=\\{\xi=(z;\,x;\,y)\mid
z=x,\;(x,y)\in{\mathcal{S}}\\},$

where ${\mathcal{S}}$ is defined in (4). By Assumption 1.1, ${\mathcal{S}}$,
and thus $\mbox{Fix}(\mathbb{G})$, is nonempty. Consequently, it follows from
(Bauschke2011book, , Proposition 5.15) that the sequence
$\\{\xi_{n}\\}_{n\in{\mathbb{N}}}$ is Fejér monotone with respect to
$\text{Fix}(\mathbb{G})$ and converges to a point in $\text{Fix}(\mathbb{G})$.
This completes the proof.

##  6 Numerical Experiments

In this section, we present numerical results on LASSO
Tibshirani1996Regression , nonnegative least-squares and minimax matrix game
problems to demonstrate the performance of Algorithm 3.1 (GRPDA), the
accelerated GRPDA given in Algorithm 4.1 (A-GRPDA), and the relaxed GRPDA
described in Algorithm 5.1 (R-GRPDA). All the experiments were performed
within Python 3.8 running on a 64-bit Windows PC with an Intel(R) Core(TM)
i5-4590<EMAIL_ADDRESS>GHz and 8GB of RAM. All the results presented in this section
are reproducible by specifying the $seed$ of the random number generator in
our code, which is available at https://github.com/cxk9369010/Golden-Ratio-
PDA.

In this section, we let $\|\cdot\|=\|\cdot\|_{2}$ be the $\ell_{2}$-norm
induced by the dot inner product. For LASSO and nonnegative least-squares
problems, the component function $f$ has the form $f(\cdot)={1\over
2}\|\cdot-b\|^{2}$. Thus, Assumption 5.1 is satisfied and the relaxed GRPDA
described in Algorithm 5.1 can be applied. Furthermore, GRPDA with fixed step
sizes (Algorithm 3.1) can adopt the larger value $\psi=2$ as well (equivalent
to setting $\psi=2$ and $\rho_{n}\equiv 1$ in Algorithm 5.1). Since $f^{*}$
(rather than $g$, as required by Algorithm 4.1) is strongly convex, we can
switch $(g,K,x,q)$ with $(f^{*},-K^{T},y,p)$ and reduce (43) to (2), so that
$g$ is strongly convex. As such, the accelerated GRPDA given in Algorithm 4.1
is also applicable.

The algorithmic parameters are specified as follows. For GRPDA, we used fixed
step sizes $\tau=\frac{\sqrt{\psi}}{\sqrt{\beta}L}$ and $\sigma=\beta\tau$.
For the minimax matrix game problem, we set $\psi=1.618$, while for LASSO and
nonnegative least-squares problems, which satisfy Assumption 5.1, we set
$\psi=2$. For A-GRPDA, we set $\psi=1.5$ and $\beta_{0}=1$, while for R-GRPDA
we set $\tau=\frac{\sqrt{2}}{\sqrt{\beta}L}$, $\sigma=\beta\tau$ and
$\rho_{n}\equiv 1.49$. The following state-of-the-art algorithms are compared:

  * •

PDA (as given in (11)) with $\tau=\frac{1}{\sqrt{\beta}L}$ and
$\sigma=\sqrt{\beta}\frac{1}{L}$ and $\delta=1$.

  * •

PGM (proximal gradient method, e.g., Beck2009A ) with fixed step size
$\alpha>0$. For LASSO and nonnegative least-squares problems, it holds that
$f(Kx)={1\over 2}\|Kx-b\|^{2}$, and the PGM appears as $x_{n}={\rm
Prox}_{\alpha g}(x_{n-1}-\alpha K^{\top}(Kx_{n-1}-b))$. We set
$\alpha=1/\|K\|^{2}$.

  * •

FISTA (fast iterative shrinkage thresholding algorithm Beck2009A ) with fixed
step size $\alpha>0$. For LASSO and nonnegative least-squares problems, FISTA
appears as

$\displaystyle\left\\{\begin{array}[]{rcl}t_{n}&=&\Big{(}1+\sqrt{1+4t_{n-1}^{2}}\Big{)}/2,\\\
y_{n-1}&=&x_{n-1}+{t_{n-1}-1\over t_{n}}(x_{n-1}-x_{n-2}),\\\ x_{n}&=&{\rm
Prox}_{\alpha g}(y_{n-1}-\alpha K^{\top}(Ky_{n-1}-b)),\end{array}\right.$

where $t_{0}=1$ and $x_{-1}=x_{0}$. The same as PGM, we set
$\alpha=1/\|K\|^{2}$.

  * •

GRAAL (golden ratio algorithm Malitsky2019Golden , or as given in (15)) with
$\tau=\frac{\phi}{2L}$, where $\phi=1.618$.

The choice of $\beta$ for GRPDA, R-GRPDA and PDA will be specified later. For
all algorithms, we used the same initial points. For the three tested
problems, the evaluation of the proximal point mappings is relatively cheap
compared with matrix-vector multiplications. As such, all the compared
algorithms have the same dominant per iteration cost, i.e., two matrix-vector
multiplications, we only compare the their performance in terms of iteration
numbers.

######  Problem 6.1 (LASSO)

The LASSO problem has the form
$\min_{x}F(x):=\frac{1}{2}\|Kx-b\|^{2}+\mu\|x\|_{1}$, where
$K\in\mathbb{R}^{p\times q}$ and $b\in\mathbb{R}^{p}$ are given, and
$x\in\mathbb{R}^{q}$ is an unknown signal. In the context of compressive
sensing, $p$ is much smaller than $q$.

Apparently, the above LASSO problem corresponds to (1) with
$f(u)=\frac{1}{2}\|u-b\|^{2}$ and $g(x)=\mu\|x\|_{1}$. Thus, Assumption 5.1 is
satisfied and $\psi=2$ can be adopted in GRPDA. The $seed$ of the random
number generator was set to $1$. We generated $x^{*}\in\mathbb{R}^{q}$
randomly. Specifically, $s$ nonzero components of $x^{*}$ were determined
uniformly at random, and their values were drawn from the uniform distribution
in $[-10,10]$. The matrix $K\in{\mathbb{R}}^{p\times q}$ is constructed as in
Malitsky2018A  by one of the following ways:

  * (i)

All entries of $K$ were generated independently from ${\mathcal{N}}(0,1)$, the
normal distribution with mean $0$ and standard deviation $1$.

  * (ii)

First, we generated a matrix $A\in\mathbb{R}^{p\times q}$, whose entries are
independently drawn from ${\mathcal{N}}(0,1)$. Then, for a scalar $v\in(0,1)$
we constructed the matrix $K$ column by column as follows:
$K_{1}=A_{1}/\sqrt{1-v^{2}}$ and $K_{j}=vK_{j-1}+A_{j}$, $j=2,\ldots,q$. Here
$K_{j}$ and $A_{j}$ represent the $j$th column of $K$ and $A$, respectively.
As $v$ increases, $K$ becomes more ill-conditioned. In this experiment, we
tested $v=0.5$ and $0.9$.

The additive noise $\nu\in{\mathbb{R}}^{p}$ was generated from
${\mathcal{N}}(0,0.1)$. Finally, we set $b=Kx^{*}+\nu$.

The step size ratio $\beta=\sigma/\tau$ was set to be $400$ as in
Malitsky2018A , which was determined based on numerical experience. All the
algorithms were initialized at $x_{0}=0$ and $y_{0}=-b$. For the LASSO
problem, we first compared the performance of GRPDA with varying $\psi$. The
results are given in Figure 1, where the decreasing behavior of function value
errors $F(x_{n})-F^{*}$ as the algorithm proceeded is presented. Here
$F^{*}=\inf_{x}F(x)$ was approximately computed by running our algorithms for
sufficiently many iterations. It can be seen from Figure 1 that GRPDA with
larger $\psi$ converges faster for all the tested three cases. Thus, in our
comparison with other algorithms, we set $\psi=2$. Besides, instead of
terminating the algorithms with some stopping criteria, we ran all algorithms
for a fixed number of iterations and examined their convergence behavior, see
Figure 2. The parameters $p$, $q$, $s$ and $v$ are given in the captions of
the figures.

(a) Case (i). (b) Case (ii) with $v=0.5$. (c) Case (ii) with $v=0.9$.

Figure 1: Evolution of $F(x_{n})-F^{*}$ for the LASSO problem with different
$\psi$ and $(p,q,s)=(200,1000,10)$.

It can be seen from Figure 2 that primal-dual type methods outperform the PGM,
FISTA and GRAAL. For all the three tested cases, GRPDA performs better than
PDA, R-GRPDA is even better and A-GRPDA is the best. In particular, the
proposed Gauss-Seidel type golden ratio algorithms, including GRPDA, A-GRPDA
and R-GRPDA, perform much better than the Jacobian type golden ratio algorithm
GRAAL as given in (15).

(a) Case (i) with $(p,q,s)=(200,1000,10)$. (b) Case (i) with
$(p,q,s)=(1000,2000,100)$. (c) Case (ii) with $(p,q,s,v)=(1000,5000,100,0.5)$.
(d) Case (ii) with $(p,q,s,v)=(1000,5000,100,0.9)$.

Figure 2: Numerical results for the LASSO problem.

######  Problem 6.2 (Nonnegative least-squares problem)

The nonnegative least-squares problem aims to find $x\in\mathbb{R}^{q}_{+}$
(the nonnegative orthant) such that $\|Kx-b\|$ is minimized for given
$K\in\mathbb{R}^{p\times q}$ and $b\in\mathbb{R}^{p}$, i.e.,
$\min_{x\in\mathbb{R}^{q}}F(x):=\frac{1}{2}\|Kx-b\|^{2}+\iota_{\mathbb{R}^{q}_{+}}(x)$.

Similar to the LASSO problem, the nonnegative least-squares problem can also
be solved by A-GRPDA and R-GRPDA since the data fitting term
$\frac{1}{2}\|Kx-b\|^{2}$ remains the same. The only difference is that the
$\ell_{1}$-norm $\mu\|\cdot\|_{1}$ in LASSO is now replaced by the indicator
function $\iota_{\mathbb{R}^{q}_{+}}(\cdot)$. Since the proximal mapping of
the indicator function $\iota_{\mathbb{R}^{q}_{+}}(\cdot)$ is just the
projection on $\mathbb{R}^{q}_{+}$, all the compared algorithms keep
feasibility of the nonnegativity constraint. As such,
$F(x_{n})=\frac{1}{2}\|Kx_{n}-b\|^{2}$ since $x_{n}\geq 0$ is always
satisfied. For this experiment, we consider two types of data described below.

  * (i)

Real data from the Matrix Market
library222https://math.nist.gov/MatrixMarket/data/Harwell-
Boeing/lsq/lsq.html.. Two instances were tested, i.e., “illc1033” and
“illc1850”, where $K\in\mathbb{R}^{p\times q}$ is sparse and has sizes
$(p,q)=(1033,320)$ and $(1850,712)$, respectively. These matrices were used in
the testing of the famous LSQR algorithm PaS82acm  and are much more ill-
conditioned than the other two instances “well1033” and “well1850” available
at the library. The entries of $b\in{\mathbb{R}}^{p}$ were drawn independently
from ${\mathcal{N}}(0,1)$.

  * (ii)

Random matrix. $K\in\mathbb{R}^{p\times q}$ has approximately $dpq$ nonzeros,
with $d\in(0,1)$ as in Malitsky2018A , and $x^{*}$ is a sparse vector, whose
$s$ nonzero entries were drawn uniformly from $[0,100]$. We set $b=Kx^{*}$,
which gives $F^{*}=0$. The following two cases were tested.

    1. (a)

$p=1000$, $q=2000$, $d=0.5$ and $s=100$. The nonzero entries of $K$ were
generated independently from the uniform distribution in $[0,1]$.

    2. (b)

$p=10000$, $q=20000$, $d=0.01$ and $s=500$. The nonzero entries of $K$ were
generated independently from the normal distribution $\mathcal{N}(0,1)$.

For GRPDA and R-GRPDA, we set $\beta=1$, except for the random data with
$p=1000$, for which we set $\beta=25$ as in Malitsky2018A . The same as the
LASSO case, all the algorithms were initialized at $x_{0}=0$ and $y_{0}=-b$.

(a) illc1033. (b) illc1850.

Figure 3: Numerical results for nonnegative least-squares problem: Real data.

(a) $p=1000$ (b) $p=10000$

Figure 4: Numerical results for nonnegative least-squares problem: Random
data.

Similar to the LASSO case, the decreasing behavior of $F(x_{n})-F^{*}$ as the
algorithms proceeded is presented in Figure 3 for real data and Figure 4 for
random data. It can be observed from Figure 3 that A-GRPDA performs the best,
followed by FISTA. In comparison, for random data R-GRPDA performs the best,
followed by GRPDA and PDA, as shown in Figure 4. Other algorithms seem to be
less efficient. Nonetheless, all the compared algorithms perform favorably and
have attained fairly high precision in a reasonable number of iterations.

######  Problem 6.3 (Minimax matrix game)

Let $\Delta_{q}=\\{x\in\mathbb{R}^{q}:x\geq 0,\;\sum_{i=1}^{q}x_{i}=1\\}$ be
the standard unit simplex in $\mathbb{R}^{q}$, and $K\in\mathbb{R}^{p\times
q}$. The minimax matrix game problem is given by

$\displaystyle\min_{x\in\mathbb{R}^{q}}\max_{y\in\mathbb{R}^{p}}\iota_{\Delta_{q}}(x)+\left\langle
Kx,y\right\rangle-\iota_{\Delta_{p}}(y).$ (95)

Apparently, (95) is a special case of (2) with $g=\iota_{\Delta_{q}}$ and
$f^{*}=\iota_{\Delta_{p}}$. We note that the projection onto $\Delta_{q}$ can
be computed efficiently, see, e.g., (Beck2017book, , Corollary 6.29).
Therefore, any algorithms depending on the proximal operators of
$\iota_{\Delta_{q}}$ can be implemented efficiently. In our implementation, we
used the algorithm from Duchi2011Diagonal  to compute the projection onto the
standard unit simplex. To compare different algorithms, we used the primal-
dual gap function as defined in (19), which can be easily computed for a
feasible pair $(x,y)\in\Delta_{q}\times\Delta_{p}$ by
$G(x,y):=\max_{i}(Kx)_{i}-\min_{j}(K^{\top}y)_{j}$. Here subscript $i$ (or
$j$) denotes the $i$th (or $j$th) component of the underlying vector.

For this minimax matrix game problem, the only relevant algorithms discussed
at the beginning of this section are PDA, GRPDA and GRAAL. The rest are not
applicable. In this experiment, we set $\beta=1$ for PDA and GRPDA. The
initial point for all algorithms was set to be
$x_{0}=\frac{1}{q}(1,\ldots,1)^{\top}$ and
$y_{0}=\frac{1}{p}(1,\ldots,1)^{\top}$. We tested two types of random matrix
$K\in{\mathbb{R}}^{p\times q}$ with $seed=50$, i.e., (i) $(p,q)=(100,100)$ and
all entries of $K$ were generated independently from the uniform distribution
in $[-1,1]$, and (ii) $(p,q)=(100,500)$ and all entries of $K$ were generated
independently from the normal distribution ${\mathcal{N}}(0,1)$.

(a) Case (i). (b) Case (ii).

Figure 5: Numerical results for minimax matrix game problem.

The decreasing behavior of the primal-dual gap function as the algorithms
proceeded is given in Figure 5, from which it can be seen that in both tests
GRPDA performs comparably with PDA and both GRPDA and PDA outperform GRAAL.
Recall that in this experiment $\beta$ was set to be $1$ and thus the primal
step size $\tau$ and the dual step size $\eta$ are equal. We attribute the
faster convergence of GRPDA and PDA, as compared to GRAAL, to the Guass-Seidel
nature of GRPDA (23) and PDA (11), as compared to the Jacobian type iteration
GRAAL (15).

##  7 Conclusions

In this paper, motivated by the recent work Malitsky2019Golden , we proposed,
analyzed and tested a golden ratio primal-dual algorithm (GRPDA), which can be
viewed as a new remedy to the classical Arrow-Hurwicz method. Different from
the PDA in Chambolle2011A , which uses an extrapolation or inertial step, the
proposed GRPDA uses a convex combination technique originally introduced in
Malitsky2019Golden . Global iterate convergence and ${\mathcal{O}}(1/N)$
ergodic rate of convergence, measured by the primal dual gap function, are
established. When either $g$ or $f^{*}$ is strongly convex, we managed to
modify the algorithm so that it enjoys a faster ${\mathcal{O}}(1/N^{2})$
ergodic convergence rate. For two widely used special cases, i.e., regularized
least-squares problem and linear equality constrained problem, we further show
via spectral analysis that the algorithmic parameter $\psi$ can be enlarged
from $(1,\phi]$ to $(1,2]$, which allows larger primal and dual step sizes.
Moreover, in the fixed point perspective the iterative mapping is
$2/3$-averaged and thus a relaxation step can be taken with parameter
$\rho\in(0,2/3)$. Our preliminary numerical results on LASSO, nonnegative
least-squares and minimax matrix game problems show that the proposed
algorithms perform favorably. In particular, GRPDA with fixed step sizes is
comparable with PDA in general, and the relaxed and accelerated variants can
achieve superior performance.

Some interesting issues remain to be investigated. For example, how to modify
the proposed algorithms so that they are suitable to the cases when the
operator norm $\|K\|$ is hard to evaluate, how to adapt the proposed
algorithms to more general settings, including the general linearly
constrained separable convex optimization problem
$\min_{x,y}\\{f(x)+g(y):Ax+By=b\\}$, or when the operator $K$ is nonlinear.
Moreover, for large scale finite sum problem, which are abundant in machine
learning, it is interesting to design and analyze their stochastic and/or
incremental counterparts. We leave these issues for further investigations.

## References

  * [1] H. H. Bauschke and P. L. Combettes.  Convex analysis and monotone operator theory in Hilbert spaces.  CMS Books in Mathematics/Ouvrages de Mathématiques de la SMC. Springer, Cham, second edition, 2011.
  * [2] A. Beck.  First-Order Methods in Optimization.  MOS-SIAM Series on Optimization. SIAM-Society for Industrial and Applied Mathematics, 2017.
  * [3] A. Beck and M. Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.
  * [4] D. P. Bertsekas and E. M. Gafni.  Projection methods for variational inequalities with application to the traffic assignment problem.  Math. Programming Stud., 17:139–159, 1982.
  * [5] T. Bouwmans, N. S. Aybat, and E. H. Zahzah.  Handbook of “Robust low-rank and sparse matrix decomposition: applications in image and video processing”, volume 45.  2016\.
  * [6] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations & Trends in Machine Learning, 3(1):1–122, 2010.
  * [7] A. Chambolle, M. J. Ehrhardt, P. Richtarik, and C.-B. Schonlieb.  Stochastic primal-dual hybrid gradient algorithm with arbitrary sampling and imaging application.  SIAM Journal on Optimization, 28(4):2783–2808, 2018.
  * [8] A. Chambolle and T. Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of Mathematical Imaging and Vision, 40(1):120–145, 2011\.
  * [9] A. Chambolle and T. Pock.  On the ergodic convergence rates of a first-order primal-dual algorithm.  Mathematical Programming, 159(1–2):253–287, SEP 2016\.
  * [10] S. Chen, M. A. Saunders, and D. L. Donoho.  Atomic decomposition by basis pursuit.  SIAM Review, 43(1):129–159, 2001.
  * [11] D. L. Donoho.  Compressed sensing.  IEEE Trans. Inform. Theory, 52(4):1289–1306, 2006.
  * [12] J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra.  Efficient projections onto the $l_{1}$-ball for learning in high dimensions.  In The 25th international conference on Machine learning, pages 272–279, 2008.
  * [13] E. Esser, X. Zhang, and T. F. Chan.  A general framework for a class of first order primal-dual algorithms for convex optimization in imaging science.  SIAM Journal on Imaging Sciences, 3(4):1015–1046, 2010\.
  * [14] M. Fazel, T. K. Pong, D. F. Sun, and P. Tseng.  Hankel matrix rank minimization with applications in system identification and realization.  SIAM Journal on Matrix Analysis and Applications, 34:946–977, 2013\.
  * [15] D. Gabay and B. Mercier.  A dual algorithm for the solution of nonlinear variational problems via finite-element approximations.  Computers and Mathematics with Applications, 2:17–40, 1976.
  * [16] R. Glowinski and A. Marrocco.  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.  R.A.I.R.O., R2, 9(R-2):41–76, 1975.
  * [17] S. Hayden and O. Stanley.  A low patch-rank interpretation of texture.  SIAM Journal on Imaging Sciences, 6(1):226–262, 2013.
  * [18] B. He and X. Yuan.  Convergence analysis of primal-dual algorithms for a saddle-point problem: from contraction perspective.  SIAM Journal on Imaging Sciences, 5(1):119–149, 2012\.
  * [19] E. Jonathan and P. Dimitri.  On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators.  Mathematical Programming, 55(1-3):293–318, 1992.
  * [20] P. L. Lions and B. Mercier.  Splitting algorithms for the sum of two nonlinear operators.  SIAM Journal on Numerical Analysis, 16(6):964–979, 1979.
  * [21] Y. Liu, Y. Xu, and W. Yin.  Acceleration of primal-dual methods by preconditioning and simple subproblem procedures.  arXiv:1811.08937v2, 2018.
  * [22] Y. Malitsky.  Golden ratio algorithms for variational inequalities.  Mathematical Programming, To appear.
  * [23] Y. Malitsky and T. Pock.  A first-order primal-dual algorithm with linesearch.  SIAM Journal on Optimization, 28(1):411–432, 2018.
  * [24] R. D. C. Monteiro and B. F. Svaiter.  Complexity of variants of Tseng’s modified FB splitting and Korpelevich’s methods for hemivariational inequalities with applications to saddle-point and convex optimization problems.  SIAM J. Optim., 21(4):1688–1720, 2011.
  * [25] A. Nedic and A. Ozdaglar.  Subgradient methods for saddle-point problems.  Journal of Optimization Theory & Applications, 142(1):205–228, JUL 2009.
  * [26] D. Needell and R. Ward.  Near-optimal compressed sensing guarantees for anisotropic and isotropic total variation minimization.  IEEE Transactions on Image Processing, 22(10):3941–3949, 2013.
  * [27] A. Nemirovski.  Prox-Method with Rate of Convergence O(1/t) for Variational Inequalities with Lipschitz Continuous Monotone Operators and Smooth Convex-Concave Saddle Point Problems.  SIAM J. Optim., 15(1):229–251, 2006.
  * [28] C. C. Paige and M. A. Saunders.  Lsqr: An algorithm for sparselinear equations and sparse least squares.  ACM Trans. Math. Soft., 8:43–71, 1982.
  * [29] T. Pock and A. Chambolle.  Diagonal preconditioning for first order primal-dual algorithms in convex optimization.  IEEE International Conference on Computer Vision, pages 1762–1769, 2011.
  * [30] R. T. Rockafellar.  Convex analysis.  Princeton University Press, 1970.
  * [31] R. Shefi and M. Teboulle.  Rate of convergence analysis of decomposition methods based on the proximal method of multipliers for convex minimization.  SIAM J. Optim., 24(1):269–297, 2014.
  * [32] R. Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society, 58(1):267–288, 1996.
  * [33] H. Uzawa.  Iterative methods for concave programming.  Studies in Linear and Nonlinear Programming (K. J. Arrow, L. Hurwicz and H. Uzawa, eds). Stanford University Press, Stanford, CA, 1958.
  * [34] J. Yang and Y. Zhang.  Alternating direction algorithms for $\ell_{1}$-problems in compressive sensing.  SIAM Journal on Scientific Computing, 33(1):250–278, 2011.
  * [35] M. Zhu and T. F. Chan.  An efficient primal-dual hybrid gradient algorithm for total variation image restoration.  CAM Report 08-34, UCLA, Los Angeles, CA, 2008.