# Accelerated forward-backward and
Douglas-Rachford splitting dynamics

Ibrahim K. Ozaslan and<EMAIL_ADDRESS>Mihailo R. Jovanović<EMAIL_ADDRESS>Ming Hsieh Department of Electrical and Computer Engineering, University of
Southern California, Los Angeles, CA 90089

###### Abstract

We examine convergence properties of continuous-time variants of accelerated
Forward-Backward (FB) and Douglas-Rachford (DR) splitting algorithms for
nonsmooth composite optimization problems. When the objective function is
given by the sum of a quadratic and a nonsmooth term, we establish accelerated
sublinear and exponential convergence rates for convex and strongly convex
problems, respectively. Moreover, for FB splitting dynamics, we demonstrate
that accelerated exponential convergence rate carries over to general strongly
convex problems. In our Lyapunov-based analysis we exploit the variable-metric
gradient interpretations of FB and DR splittings to obtain smooth Lyapunov
functions that allow us to establish accelerated convergence rates. We provide
computational experiments to demonstrate the merits and the effectiveness of
our analysis.

###### keywords:

Nesterov’s acceleration; proximal gradient method; forward-backward splitting;
Douglas-Rachford splitting; envelope functions; primal-dual gradient flow
dynamics.

##  1 Introduction

The convergence properties of gradient descent (GD),

$x^{k+1}\,=\,x^{k}\,-\,\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\alpha}\nabla f(x^{k})$ (1)

can be improved by adding a momentum term,

$x^{k+1}\,=\,x^{k}\,+\,\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\gamma}(x^{k}\,-\,x^{k-1})\,+\,\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\alpha}\nabla f(x^{k}\,+\,\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\beta}(x^{k}\,-\,x^{k-1})).$ (2)

where $k$ is the iteration index, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\alpha}$ is the stepsize, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\beta}$ is the extrapolation parameter, and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\gamma}$ is the damping coefficient. In particular, the Polyak’s Heavy-
Ball Method (HBM) is obtained by setting $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\beta}=0$ in (2) [22] and Nesterov’s Accelerated Algorithm (NAA) results
from $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\gamma}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\beta}$ in (2) [14]. For a convex objective function $f$:
$\mathbb{R}^{n}\to\mathbb{R}$ with an $L$-Lipschitz continuous gradient
$\nabla f$, NAA with $(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\alpha},\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\beta})=(1/L,(k-2)/(k+3))$ improves the sublinear convergence rate
$O(1/k)$ of GD to $O(1/k^{2})$. Furthermore, for an $m$-strongly convex $f$,
the linear convergence rate $O(\mathrm{e}^{-k/\kappa})$ of GD [15, Section
2.1.5] is improved to $O(\mathrm{e}^{-k/\sqrt{\kappa}})$ with
$(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\alpha},\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\beta})=(1/L,(\sqrt{\kappa}-1)/(\sqrt{\kappa}+1))$. In both cases, these
rates are optimal in the sense that there exists a smooth convex function $f$
for which the rates cannot be improved [15, Theorems 2.1.7 and 2.1.13].

A wide-spread use of accelerated algorithms has inspired efforts to understand
the underlying mechanisms. In particular, a growing body of literature views
optimization algorithms as continuous-time dynamical systems that allow
familiar methods to be recovered upon proper discretization [7, 13, 6, 1]. For
example, the continuous-time equivalent of GD is given by the gradient flow
dynamics,

$\dot{x}\,=\,-\alpha\nabla f(x)$ (3)

and celebrated Arrow-Hurwicz-Uzawa primal-dual gradient flow dynamics can be
used to effectively solve a class of constrained optimization problems [17]. A
dynamical systems perspective of NAA offered in [24] was further exploited by
leveraging a variational formulation based on Lagrangian and Hamiltonian
frameworks from classical mechanics [25, 26] and contraction theory [5].
Furthermore, NAA was obtained using the semi-implicit Euler discretization of
the differential equation [13],

$\ddot{x}\,+\,\gamma\dot{x}\,+\,\alpha\nabla f(x\,+\,\beta\dot{x})\,=\,0$ (4)

where $\alpha$, $\beta$, and $\gamma$ are possibly time-varying algorithmic
parameters that differ for convex and strongly convex problems. The following
Lyapunov function,

$V(x,\dot{x})\,=\,af(x)\,+\,(1/2)\|bx\,+\,\dot{x}\|_{2}^{2}$ (5)

with positive parameters $a$ and $b$ was also utilized in [13] to show that
(4) achieves the optimal convergence rate of discrete-time algorithms for both
convex and strongly convex problems with proper choices of
($\alpha,\beta,\gamma$). Moreover, for $\beta+\gamma=1$, the discretization of
(4) preserves the accelerated convergence rate. This provides a system-
theoretic characterization of acceleration and bridges the gap between
discrete and continuous-time analyses.

Our contributions. We generalize accelerated dynamics (4) to nonsmooth
composite optimization problems in which the objective function is given by
the sum of a smooth and a nonsmooth term. Our approach is motivated by [10]
which obtains the Forward-Backward (FB) and Douglas-Rachford (DR) splitting
dynamics by replacing $\nabla f$ in (3) with the generalized gradient map
associated with the nonsmooth objective function. Based on the relation
between gradient flow dynamics (3) and its accelerated variant (4), we
introduce accelerated FB and DR splitting dynamics by replacing gradient
$\nabla f$ in (4) with the generalized gradient map. We also utilize a
Lyapunov-based approach to establish accelerated sublinear and linear
convergence rates that match the optimal discrete-time rates for convex and
strongly convex problems, respectively. Our Lyapunov function is similar to
(5) with two important differences: (i) the Euclidean distance is replaced
with a suitable weighted norm; and (ii) the FB or DR envelopes are used
instead of the objective function. The latter choice is motivated by the
variable-metric gradient method interpretation of associated splittings [21,
20]. We also note that, in our analysis, the condition $\beta+\gamma=1$ is
satisfied, which is an integral requirement for discretization analysis of
[13] to be applicable to the proposed accelerated dynamics.

Although accelerated rates for the FB and DR splittings are known in discrete-
time, the continuous-time analysis has not been done before. A differential
inclusion model of accelerated FB splitting was introduced in [4, 2], but no
analysis was provided. Furthermore, the continuous-time analysis of
accelerated FB splitting resorts to subgradients and it only captures
sublinear convergence rates [24, Thm. 24].

It is also worth mentioning that the FB and DR splittings are special cases of
the primal-dual gradient flow dynamics associated with the proximal augmented
Lagrangian [6, 10, 16]. Thus, insights about acceleration of these splitting
methods can also pave the way to accelerated primal-dual algorithms, which is
of critical importance to the design of fast gradient-based methods for
constrained optimization problems. We note that (4) was used as a motivation
in [3, 27, 11] to accelerate primal-dual gradient flow dynamics. However, in
contrast to [13], the algorithmic parameters are not bounded in time in these
references. Thus, it is unclear if the primal-dual dynamics proposed in [3,
27, 11] allow rate-preserving discretization that leads to efficient first-
order iterative schemes.

Paper structure. In Section 2, we introduce the problem setup and provide
background material. In Section 3, we summarize our main findings, and in
Section 4 we provide the proofs. In Section 5, we utilize computational
experiments to demonstrate the merits of our analyses, and in Section 6, we
conclude the paper with remarks.

##  2 Problem formulation and preliminaries

In this section, we formulate the problem and provide background on the
proximal operators; Moreau, FB, and DR envelopes; as well as FB and DR
splittings.

We consider convex optimization problem

$\operatorname*{minimize}\limits_{x}~{}F(x)\;\mathrel{\mathop{:}}=\;f(x)\,+\,g(x)$
(6)

where $f$: $\mathbb{R}^{n}\to\mathbb{R}$ is, unless noted otherwise, a
continuously differentiable convex function with an $L$-Lipschitz continuous
gradient $\nabla f$ and $g$: $\mathbb{R}^{n}\to\mathbb{R}$ is a possibly non-
differentiable closed proper convex function.

###  2.1 Proximal operator and Moreau envelope

The proximal operator of a proper lower semi-continuous convex function $g$
for a positive penalty parameter $\mu$ is the mapping defined as
$\operatorname{\mathbf{prox}}_{\mu g}(v)=(I+\mu\partial g)^{-1}(v),$ where
$\partial g$ is a subdifferential of $g$ and $I$ is the identity map [19].
Alternatively, $\operatorname{\mathbf{prox}}_{\mu g}$ can be obtained as

$\operatorname{\mathbf{prox}}_{\mu
g}(v)\,=\,\operatorname*{argmin}\limits_{z}\left(g(z)\,+\,\tfrac{1}{2\mu}\|z\,-\,v\|_{2}^{2}\right).$

The value function of this optimization problem determines the Moreau envelope
associated with $g$,

$\begin{array}[]{rcl}\\!\\!\\!\mathcal{M}_{\mu
g}(v)&=&g(\operatorname{\mathbf{prox}}_{\mu
g}(v))\,+\,\frac{1}{2\mu}\|\operatorname{\mathbf{prox}}_{\mu
g}(v)\,-\,v\|_{2}^{2}\end{array}$ (7)

which is, even for a non-differentiable $g$, a continuously differentiable
function [19],

$\nabla\mathcal{M}_{\mu
g}(v)\,=\,\tfrac{1}{\mu}(v\,-\,\operatorname{\mathbf{prox}}_{\mu
g}(v))\,\in\,\partial g(\operatorname{\mathbf{prox}}_{\mu g}(v)).$ (8)

###  2.2 The FB and DR envelopes

Using Lipschitz continuity of $\nabla f$, we have the following upper bound on
the objective function $F$ in (6) [21],

$\begin{array}[]{l}\\!\\!F(y)\leq\;J(x,y)\;\mathrel{\mathop{:}}=\\\
\\!\\!f(x)\,+\,g(y)\,+\,\tfrac{1}{2\mu}\|y\,-\,(x\,-\,\mu\nabla
f(x))\|_{2}^{2}\,-\,\frac{\mu}{2}\|\nabla f(x)\|_{2}^{2}\end{array}$

which holds for any $x,y\in\mathbb{R}^{n}$ and $\mu\in(0,1/L)$. Minimizing
$J(x,y)$ over $y$ yields the FB envelope of $F$ [21]

$\begin{array}[]{rrl}\mathcal{F}_{\mu}(x)&\mathrel{\mathop{:}}=&\operatorname*{minimize}\limits_{y}~{}J(x,y)\\\
&=&f(x)\,+\,\mathcal{M}_{\mu g}(x\,-\,\mu\nabla
f(x))\,-\,\tfrac{\mu}{2}\|\nabla f(x)\|_{2}^{2}.\end{array}$

For a particular choice of a distance function (instead of the Euclidean
norm), the FB envelope can also be interpreted as a generalized Moreau
envelope of $F$ [12]. The DR envelope associated with (6) is another useful
value function resulting from evaluation of the FB envelope
$\mathcal{F}_{\mu}(x)$ at $x=\operatorname{\mathbf{prox}}_{\mu f}(z)$ with
$z\in\mathbb{R}^{n}$ [20, 8], i.e.,

$\mathcal{D}_{\mu}(z)\,\mathrel{\mathop{:}}=\,\mathcal{F}_{\mu}(\operatorname{\mathbf{prox}}_{\mu
f}(z)).$

When $f$ is a twice continuously differentiable function with a bounded
Hessian $\nabla^{2}f(x)$, the FB envelope is also continuously differentiable
with gradient

$\nabla\mathcal{F}_{\mu}(x)\,=\,(I\,-\,\mu\nabla^{2}f(x))G_{\mu}(x)$ (9)

where $G_{\mu}$ is the generalized gradient map defined as

$G_{\mu}(x)\,=\,\tfrac{1}{\mu}(x\,-\,\operatorname{\mathbf{prox}}_{\mu
g}(x\,-\,\mu\nabla f(x))).$ (10)

Consequently, the DR envelope is also continuously differentiable and its
gradient is determined by

$\nabla\mathcal{D}_{\mu}(z)\,=\,\nabla\operatorname{\mathbf{prox}}_{\mu
f}(z)\nabla\mathcal{F}_{\mu}(\operatorname{\mathbf{prox}}_{\mu f}(z)).$ (11)

Here, $\nabla\operatorname{\mathbf{prox}}_{\mu f}(z)$ can be obtained by
differentiating (8)

$\nabla\operatorname{\mathbf{prox}}_{\mu
f}(z)\;\mathrel{\mathop{:}}=\;(I\,+\,\mu\nabla^{2}f(\operatorname{\mathbf{prox}}_{\mu
f}(z)))^{-1}$ (12)

Using (9) and (11), $\nabla\mathcal{D}_{\mu}$ can also be expressed in terms
of the generalized gradient map as

$\nabla\mathcal{D}_{\mu}(z)\,=\,(2\nabla\operatorname{\mathbf{prox}}_{\mu
f}(z)\,-\,I)G_{\mu}(\operatorname{\mathbf{prox}}_{\mu f}(z)).$ (13)

We note that for $\mu\in(0,~{}1/L)$,

$\begin{array}[]{c}\min~{}F\,=\,\min~{}\mathcal{F}_{\mu}\,=\,\min~{}\mathcal{D}_{\mu}\\\
\operatorname*{argmin}\;F\,=\,\operatorname*{argmin}\;\mathcal{F}_{\mu}\,=\,\operatorname{\mathbf{prox}}_{\mu
f}(\operatorname*{argmin}\;\mathcal{D}_{\mu})\end{array}$ (14)

which together with the smoothness of the envelope functions motivates
studying acceleration of dynamics (16) and (20) through the FB and DR
envelopes [21].

###  2.3 The FB and DR splittings

Proximal gradient algorithm [19],

$x^{k+1}\,=\,x^{k}\,-\,\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\alpha}G_{\mu}(x^{k})$ (15)

replaces $\nabla f$ in (1) with the generalized gradient map $G_{\mu}$,
thereby extending gradient descent to composite optimization problem (6).
Iterations (15) are also known as the FB splitting and their continuous-time
equivalent is determined by the proximal gradient flow dynamics [10],

$\dot{x}\,=\,-\alpha G_{\mu}(x).$ (16)

Furthermore, evaluation of the generalized gradient map at
$x=\operatorname{\mathbf{prox}}_{\mu f}(z)$ facilitates extension of FB
splitting (15) to composite optimization problem (6) in which both $f$ and $g$
are allowed to be non-differentiable,

$z^{k+1}\;\in\;z^{k}\,-\,\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\alpha}G_{\mu}(\operatorname{\mathbf{prox}}_{\mu f}(z^{k})).$ (17)

Using (10) and setting $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\alpha}=\mu$ in (17) yields,

$\begin{array}[]{rcl}z^{k+1}&\in&z^{k}\,-\,\operatorname{\mathbf{prox}}_{\mu
f}(z)\;+\\\ &&\operatorname{\mathbf{prox}}_{\mu
g}(\operatorname{\mathbf{prox}}_{\mu f}(z)\,-\,\mu\partial
f(\operatorname{\mathbf{prox}}_{\mu f}(z))).\end{array}$ (18)

Moreover, (8) can be used to rewrite (18) as,

$z^{k+1}\,=\,z^{k}\,-\,\operatorname{\mathbf{prox}}_{\mu
f}(z)\;+\;\operatorname{\mathbf{prox}}_{\mu
g}(2\operatorname{\mathbf{prox}}_{\mu f}(z)-z)$ (19)

which is known as DR splitting [23]. While (6) can be solved via (19) even for
a non-differentiable $f$, when $f$ is continuously differentiable
$G_{\mu}(\operatorname{\mathbf{prox}}_{\mu f}(z))$ becomes a single-valued
function and (17) can be used to obtain the DR splitting dynamics [10],

$\dot{z}\,=\,-\alpha G_{\mu}(\operatorname{\mathbf{prox}}_{\mu f}(z)).$ (20)

Lemma 1 characterizes stability properties of proximal gradient flow (16) and
DR splitting dynamics (20) under the following assumption.

######  Assumption 1.

In (6), $f$ is an $m$-strongly convex function with an $L$-Lipschitz
continuous gradient $\nabla f$ and $g$ is a proper closed convex function.

######  Lemma 1 (Thms. 3 and 7, [10]).

Let Assumption 1 hold with $m>0$ and let $\mu=1/(2L)$. Then, the proximal
gradient flow (16) and DR splitting dynamics (20) are globally exponentially
stable with rate $\rho=\alpha m$.

##  3 Main results

In this section, we examine convergence properties of continuous-time variants
of accelerated proximal gradient and DR splitting algorithms for nonsmooth
composite optimization problem (6). In Theorems 1 and 2, we show that the
corresponding dynamical systems provide acceleration when the smooth part of
the objective function $f$ is either convex or strongly convex quadratic
function, respectively. Furthermore, in Theorem 3, we demonstrate that
properties of accelerated proximal gradient flow dynamics carry over from
quadratic to general strongly convex problems. To improve readability, we
relegate the proofs of theorems to Sections 4.

Inspired by (4), which provides accelerated convergence for smooth problems,
we propose the following differential equation to accelerate proximal gradient
flow dynamics (16),

rcl ¨x + γ˙x \+ αG_μ(x + β˙x) & = 0. Similarly, we introduce rcl ¨z + γ˙z \+
αG_μ(prox_μf(z + β˙z)) & = 0
x = prox_μf(z)

as an accelerated variant of DR splitting dynamics (20), where $G_{\mu}$ is
the generalized gradient map associated with problem (6); see (10) for
definition.

Our analysis of (21) and (21) utilizes Lyapunov functions obtained by
replacing $f$ in (5) with smooth envelope functions $\mathcal{F}_{\mu}$ or
$\mathcal{D}_{\mu}$, respectively. This choice is motivated by the variable-
metric gradient method interpretations of FB and DR splittings. The lack of
convexity of the FB and DR envelopes when $f$ in (6) is not a convex quadratic
function is the main challenge in convergence analysis. When $f$ is either a
convex or a strongly convex quadratic function, we show that (21) and (21)
converge to a solution of problem (6) at an accelerated rate. We then remove
this quadratic restriction on $f$ and show that acceleration carries over to
general strongly convex problems for FB splitting dynamics (21).

###  3.1 Convex quadratic $f$

We first restrict our attention to a nonsmooth optimization problem (6) with
quadratic $f$. In Theorem 1, we establish accelerated sublinear convergence
for convex problems and, in Theorem 2, we show that the convergence rate can
be improved to exponential for strongly convex problems. This extends results
from the literature from a discrete-time to a continuous-time setting.

######  Assumption 2.

In (6), $f$ is a convex quadratic function with a strong convexity constant
$m\geq 0$ and a Lipschitz continuity constant $L$, and $g$ is a proper closed
convex function.

######  Theorem 1.

Let Assumption 2 hold and let $\mu\in(0,1/L)$. If the algorithmic parameters
are given by $(\alpha>0,\gamma(t)=3/(t+3),\beta(t)=1-\gamma(t))$ then the
solutions to (21) and (21) for $t\geq 0$ satisfy

$F(p_{\mu}(x(t)))-F(x^{\star})\,\leq\,\tfrac{c_{1}}{(t\,+\,3)^{2}}(\|x(0)-x^{\star}\|_{2}^{2}\,+\,\|\dot{x}(0)\|_{2}^{2})$

where $p_{\mu}(x)\mathrel{\mathop{:}}=\operatorname{\mathbf{prox}}_{\mu
g}(x-\mu\nabla f(x))$, $x^{\star}\in\operatorname*{argmin}F(x)$, and $c_{1}>0$
is a constant that depends on the problem parameters.

Theorem 1 establishes accelerated sublinear convergence rate for convex
problems with quadratic smooth part of the objective function. This matches
the optimal rate $O(1/k^{2})$ of discrete-time gradient-based methods [14].
Furthermore, in Theorem 2 we show that the use of constant algorithmic
parameters for strongly convex problems with quadratic $f$ and known $m>0$
leads to an accelerated exponential convergence rate. In proof of Theorem 2,
we utilize Lemma 2 which provides estimates of the smoothness and strong
convexity constants of the FB and DR envelopes under Assumption 2.

######  Lemma 2 (Thm. 2.3 [21] and Proposition 4.6 [8]).

Let Assumption 2 hold and let $\mu\in(0,1/L)$. The FB and DR envelopes
$\mathcal{F}_{\mu}$ and $\mathcal{D}_{\mu}$ are convex functions with the
following smoothness and strong convexity constants

$\begin{array}[]{rcl}\tilde{L}&=&\left\\{\begin{array}[]{ll}2(1\,-\,\mu
m)/\mu&~{}~{}\mbox{\rm for}~{}{\mathcal{F}_{\mu}}\\\\[2.84544pt] (1\,-\,\mu
m)/\big{(}\mu(1\,+\,\mu m)^{2}\big{)}&~{}~{}\mbox{\rm
for}~{}{\mathcal{D}_{\mu}}\end{array}\right.\\\\[14.22636pt]
\tilde{m}&=&\left\\{\begin{array}[]{ll}\min\left\\{(1\,-\,\mu m)m,(1\,-\,\mu
L)L\right\\}&~{}~{}\mbox{\rm for}~{}{\mathcal{F}_{\mu}}\\\\[2.84544pt]
\min\left\\{\dfrac{(1\,-\,\mu m)m}{(1\,+\,\mu m)^{2}},\dfrac{(1\,-\,\mu
L)L}{(1\,+\,\mu L)^{2}}\right\\}&~{}~{}\mbox{\rm
for}~{}{\mathcal{D}_{\mu}}\end{array}\right.\end{array}$

If the condition number $\kappa=L/m$ of $f$ satisfies $\kappa\geq 2$, then
Lemma 2 shows that for $\mu=1/(kL)$ with $k\geq 2$, the condition numbers of
the FB and DR envelopes,
$\tilde{\kappa}\mathrel{\mathop{:}}=\tilde{L}/\widetilde{m}$, are both upper
bounded by $2k\kappa$.

######  Theorem 2.

Let Assumption 2 hold with $m>0$, let $\mu\in(0,1/L)$, and let $\tilde{L}$ and
$\widetilde{m}$ be the smoothness and strong convexity constants of either the
FB or the DR envelope. If the algorithmic parameters in (21) are given by
$(\alpha>0,\gamma=2\sqrt{\alpha\widetilde{m}}\big{/}(\sqrt{\alpha\widetilde{m}}+1),\beta=1-\gamma)$
then the solutions to (21) and (21) for $t\geq 0$ satisfy

$\|x(t)\,-\,x^{\star}\|_{2}^{2}\;\leq\;c_{2}(\|x(0)\,-\,x^{\star}\|_{2}^{2}\,+\,\|\dot{x}(0)\|_{2}^{2})\,\mathrm{e}^{-\rho
t}$

where $\rho=\sqrt{\alpha\widetilde{m}}-\alpha\widetilde{m}/2$ is the
convergence rate and $c_{2}>0$ is a constant that depends on the problem
parameters.

###  3.2 Strongly convex $f$

We next show that an additional restriction on $\mu$ allows us to obtain an
accelerated exponential convergence rate for (21) even when the smooth part of
the objective function $f$ in (6) is not quadratic but strongly convex.

######  Theorem 3.

Let Assumption 1 hold and let $\mu\in\left(0,\sqrt{\gamma\beta}/(2L)\right]$.
If the algorithmic parameters in (21) are given by
$(\alpha>0,\gamma=2\sqrt{\alpha m}\big{/}(\sqrt{\alpha m}+1),\beta=1-\gamma)$
then the solutions to (21) for $t\geq 0$ satisfy

$\|x(t)\,-\,x^{\star}\|_{2}^{2}\;\leq\;c_{3}(\|x(0)\,-\,x^{\star}\|_{2}^{2}\,+\,\|\dot{x}(0)\|_{2}^{2})\mathrm{e}^{-\rho
t}$

where $\rho=\sqrt{\alpha m}-\alpha m/2$ is the convergence rate, and $c_{3}>0$
is a constant that depends on the problem parameters.

Setting $\alpha=1/\tilde{L}$ and $\alpha=1/L$ in Theorems 2 and 3,
respectively, gives an exponential convergence rate, $\rho\sim
O(1/\sqrt{\kappa})$. This rate matches the best-achievable rate of discrete-
time gradient-based methods for strongly convex problems [15]. Moreover, for
$\alpha=1/L$, the upper bound on the penalty parameter $\mu$ in Theorem 3 is
of $O(1/(L\sqrt[4]{\kappa}))$ and is smaller than the upper bound on $\mu$ in
Theorem 2 by a factor of $O(1/\sqrt[4]{\kappa})$.

######  Remark 1.

For the rate-preserving discretization proposed in [13] to be applicable to
accelerated dynamics (21) and (21), we require $\beta+\gamma=1$ in Theorems
1-3. For general convex problems, this additional constraint on algorithmic
parameters together with the lack of convexity of envelope functions makes the
analysis challenging. For (21), we obtain acceleration by introducing
additional constraints on the penalty parameter $\mu$ and two novel
inequalities that substitute the strong convexity inequalities; see Lemma 3
for details.

##  4 Proof of main results

We next present a system-theoretic framework that can be specialized to
particular problem setups addressed in Theorems 1-3. To simplify notation, we
assume without loss of generality that in (21),
$\min_{x}\mathcal{F}_{\mu}(x)=0$ and $x^{\star}=0$ is a minimizer; similarly,
for (21), $\min_{z}\mathcal{D}_{\mu}(z)=0$ and $z^{\star}=0$ is a minimizer.

System representation. We rewrite (21) as a feedback of a linear system with a
static nonlinear block,

$\dot{\psi}\,=\,A\psi\,+\,Bu,\quad y\,=\,C\psi,\quad u\,=\,\Delta(y)$ (22a)
where $\psi=[~{}x^{T}~{}\dot{x}^{T}~{}]^{T}$ and
$\psi=[~{}z^{T}~{}\dot{z}^{T}~{}]^{T}$ denote the respective state vectors for
(21) and (21). On the other hand, the system matrices are given by
$\left[\begin{array}[]{c|c}A&B\\\ \hline\cr
C&0\end{array}\right]\,=\,\left[\begin{array}[]{cc|c}0&I&0\\\ 0&-\gamma
I&-\alpha I\\\ \hline\cr I&\beta I&0\end{array}\right].$ (22b) Here, the
partition is done conformably with the partition of the state vector
$\psi=[~{}\psi_{1}^{T}~{}\psi_{2}^{T}~{}]^{T}$ and the corresponding nonlinear
terms are $\Delta(y)\,=\,\left\\{\begin{array}[]{ll}G_{\mu}(y),&\text{FB
splitting}\\\\[-1.42271pt] G_{\mu}(\operatorname{\mathbf{prox}}_{\mu
f}(y)),&\text{DR splitting.}\end{array}\right.$ (22c)

###  4.1 Theorems 1 and 2: convex quadratic $f$

We first prove our results for convex quadratic $f$ in (6),
$f(x)=(1/2)x^{T}Qx+q^{T}x$, where $Q\in\mathbb{R}^{n\times n}$ is a positive
semidefinite matrix and $q\in\mathbb{R}^{n}$. In this case, the Lipschitz
continuity and strong convexity parameters $L$ and $m$ are given by the
largest and smallest eigenvalues of $Q$.

Lyapunov function. Inspired by (5) which was used to establish accelerated
convergence rates for (4) in [13], we propose the Lyapunov function candidate
for (21),

$V(t;\psi)\,=\,\alpha\mathcal{E}(\psi_{1})\,+\,(1/2)\|\theta(t)\psi_{1}+\psi_{2}\|_{H}^{2}$
(23)

where the envelope $\mathcal{E}$ (either $\mathcal{F}_{\mu}$ or
$\mathcal{D}_{\mu}$) and the associated weight matrix $H$ are determined by
the selection of the nonlinear term $\Delta$ in (22c). Our aim is to show that
$\dot{V}(t)+\theta(t)V(t)\leq 0$ for $t\geq t_{0}$ which in conjunction with
the Grönwall inequality implies $V(t)\leq
V(t_{0})\exp(-\int_{t_{0}}^{t}\theta(\tau)\mathrm{d}\tau)$. The inequalities
presented in the theorems are then obtained by leveraging various properties
of associated envelopes and the proximal operators. Hence, $\theta(t)$ in (23)
determines the convergence rate. As we elaborate in the forthcoming
subsections, in Theorem 1 we set $\theta(t)=2/(t+3)$ and in Theorem 2 we set
$\theta(t)=\rho$ to establish sublinear and exponential convergence rates,
respectively.

Characterization of FB and DR envelopes. To upper bound $\theta V$, we first
derive an upper bound on the envelope function $\mathcal{E}$. For a (strongly)
convex quadratic $f$, the FB and DR envelopes $\mathcal{E}$ are also
(strongly) convex; see Lemma 2. Thus, for any $x,x^{\prime}\in\mathbb{R}^{n}$,

$\mathcal{E}(x)\,-\,\mathcal{E}(x^{\prime})\,\leq\,\langle\nabla\mathcal{E}(x),x\,-\,x^{\prime}\rangle\,-\,\dfrac{\widetilde{m}}{2}\,\|x\,-\,x^{\prime}\|_{2}^{2}$
(24)

where $\widetilde{m}$ is the parameter of strong convexity of $\mathcal{E}$.
Moreover, (9) and (13) yield $\nabla\mathcal{E}(x)=H\Delta(x)$. For the FB
envelope $\mathcal{E}=\mathcal{F}_{\mu}$, $\Delta(x)=G_{\mu}(x)$, and

$H\,=\,I\,-\,\mu\nabla^{2}f(x)\,=\,I\,-\,\mu Q.$

On the other hand, for the DR envelope $\mathcal{E}=\mathcal{D}_{\mu}$,
$\Delta(x)=G_{\mu}(\operatorname{\mathbf{prox}}_{\mu f}(x))$ and the matrix
$H$ is determined by

$H\,=\,2\nabla\operatorname{\mathbf{prox}}_{\mu f}(x)\,-\,I\,=\,(I\,-\,\mu
Q)(I\,+\,\mu Q)^{-1}.$

In both cases, $0\prec H\preceq I$ for $\mu\in(0,1/L)$, which together with
(24) implies that for any $x,x^{\prime}\in\mathbb{R}^{n}$

$\mathcal{E}(x)\,-\,\mathcal{E}(x^{\prime})\,\leq\,\langle\Delta(x),x\,-\,x^{\prime}\rangle_{H}\,-\,\dfrac{\widetilde{m}}{2}\,\|x\,-\,x^{\prime}\|_{H}^{2}$
(25)

where $\|x\|_{H}^{2}\mathrel{\mathop{:}}=\langle x,x\rangle_{H}=x^{T}Hx$.

An upper bound on the Lyapunov function. We first rewrite the Lyapunov
function candidate using system representation (22a) as

$V(t;\psi)\,=\,\alpha\mathcal{E}(C_{1}\psi)\,+\,(1/2)\|R(t)\psi\|_{H}^{2}$
(26)

where $C_{1}=[I~{}0]$, $R(t)=[\theta(t)I~{}I]$. We define $C_{2}=[0~{}\beta
I]$, $v=\Delta(C_{1}\psi)$, and recall our assumption that
$\psi_{1}^{\star}=C_{1}\psi^{\star}=0$ and $\mathcal{E}(\psi_{1}^{\star})=0$.
Using inequality (25) twice, for $(x,x^{\prime})=(C\psi,C_{1}\psi^{\star})$
and $(x,x^{\prime})=(C_{1}\psi,C\psi)$, along with $u=\Delta(C\psi)$ and
$C=C_{1}+C_{2}$ yield,

$\begin{array}[]{rcl}\mathcal{E}(C\psi)&\leq&\langle
u,C\psi\rangle_{H}\,-\,(\widetilde{m}/2)\|C\psi\|_{H}^{2}\\\
\mathcal{E}(C_{1}\psi)\,-\,\mathcal{E}(C\psi)&\leq&-\langle
v,C_{2}\psi\rangle_{H}\,-\,(\widetilde{m}/2)\|C_{2}\psi\|_{H}^{2}.\end{array}$

Finally, using the definition of $V$ in (26) and summing the last two
inequalities gives an upper bound on $\theta V$,

$\begin{array}[]{rcl}\\!\\!\\!\\!\\!\theta V&\leq&\alpha\theta\left(\langle
u,C\psi\rangle_{H}\,-\,\langle v,C_{2}\psi\rangle_{H}\right)\\\
&&-\,\dfrac{\alpha\widetilde{m}\theta}{2}\left(\|C\psi\|_{H}^{2}\,+\,\|C_{2}\psi\|_{H}^{2}\right)\,+\,\dfrac{\theta}{2}\|R\psi\|_{H}^{2}.\end{array}$
(27)

We next examine the derivative of $V$ along the solutions of (22) and use (27)
to obtain an upper bound on $\dot{V}+\theta V$.

Lyapunov analysis. The derivative of (26) along the solutions of (22) is
determined by

$\dot{V}=\alpha\langle v,C_{1}(A\psi+Bu)\rangle_{H}+\langle
R\psi,\dot{R}\psi+R(A\psi+Bu)\rangle_{H}$ (28)

where $\dot{R}=[(\mathrm{d}\theta(t)/\mathrm{d}t)I~{}0]$. Finally, combining
(28) with (27) leads to

$\dot{V}\,+\,\theta
V\,\leq\,\frac{1}{2}\left[\begin{array}[]{c}\psi\\\\[-21.33955pt] \\\
u\\\\[-21.33955pt] \\\
v\end{array}\right]^{T}\left[\begin{array}[]{ccc}\Pi&\Gamma&\Lambda\\\\[-21.33955pt]
\\\ \Gamma^{T}&0&0\\\\[-21.33955pt] \\\
\Lambda^{T}&0&0\end{array}\right]\left[\begin{array}[]{c}\psi\\\\[-21.33955pt]
\\\ u\\\\[-21.33955pt] \\\ v\end{array}\right]$ (29)

where matrices $\Pi$, $\Gamma$, and $\Lambda$ differ for convex and strongly
convex cases. Now, we prove that the matrix on the right-hand-side of (29) is
negative semidefinite for the parameters given in Theorems 1 and 2.

####  4.1.1 Theorem 1: convex quadratic $f$

For a (non-strongly) convex quadratic function $f$, $\widetilde{m}=0$, the
upper bound on $\theta V$ given in (27) simplifies to

$\theta V\,\leq\,\alpha\theta\left(\langle u,C\psi\rangle_{H}\,+\,\langle
v,-C_{2}\psi\rangle_{H}\right)\,+\,(\theta/2)\|R\psi\|_{H}^{2}.$

We let $P\mathrel{\mathop{:}}=R^{T}HR$ and utilize this upper bound and (28)
to write $\Pi$ in (29) as

$\begin{array}[]{rcl}\Pi&=&A^{T}P\,+\,PA\,+\,\theta P\,+\,\dot{P}\\\
&=&\left[\begin{array}[]{cc}2\dot{\theta}\theta\,+\,\theta^{3}&\dot{\theta}\,+\,2\theta^{2}\,-\,\theta\gamma\\\\[-18.49411pt]
\\\
\dot{\theta}\,+\,2\theta^{2}\,-\,\theta\gamma&3\theta\,-\,2\gamma\end{array}\right]\otimes
H\end{array}$ (30)

where $\dot{P}=\dot{R}^{T}HR+R^{T}H\dot{R}$ and $\otimes$ is the Kronecker
product. A necessary condition for negative semidefiniteness of the matrix on
the right-hand-side of (29) is the negative semidefiniteness of $\Pi$. One way
to simplify the conditions that guarantee $\Pi\preceq 0$ is to set the off-
diagonal blocks in (30) to zero, i.e.,

$\dot{\theta}\,+\,2\theta^{2}\,-\,\theta\gamma\,=\,0.$ (31)

When (31) is satisfied, $\Pi\preceq 0$ if its diagonal blocks are negative
semidefinite. The first diagonal block requires $\theta\leq 2/(t+r)$ and the
largest $\theta$ that satisfies this condition is $\theta(t)=2/(t+r)$ for some
$r>0$. Substituting this value to (31) yields $\gamma(t)=3/(t+r)$. For this
choice of parameters, $\Pi\equiv 0$. While any choice of $r>0$ works here,
different constants $c_{1}$ in Theorem 1 are obtained for different values of
$r$. In what follows, we set $r=3$ to ensure positivity of $\beta(t)$ in
Theorem 1 for any $t\geq 0$.

Since the off-diagonal blocks in (29) are given by

$\begin{array}[]{rclll}\Gamma&=&PB\,+\,\alpha\theta
C^{T}H&=&\left[\begin{array}[]{c}0\\\
-\alpha(1\,-\,\theta\beta)H\end{array}\right]\\\
\Lambda&=&\dfrac{\alpha}{\beta}(1\,-\,\theta\beta)C_{2}^{T}H&=&\left[\begin{array}[]{c}0\\\
\alpha(1\,-\,\theta\beta)H\end{array}\right]\,=\,-\Gamma.\end{array}$ (32)

for the above selection of $\theta$ and $\gamma$, (29) simplifies to

$\begin{array}[]{rcl}\dot{V}\,+\,\theta
V&\leq&-(\alpha/\beta)(1\,-\,\theta\beta)\langle
u\,-\,v,C_{2}\psi\rangle_{H}.\end{array}$ (33)

Moreover, monotonicity of the gradient of the envelope function implies
$\langle u-v,C_{2}\psi\rangle_{H}\geq 0$. Thus, the right-hand-side of (33) is
upper bounded by zero and the use of Grönwall inequality yields $V(t)\leq
V(0)/(t+3)^{2}$, $t\geq 0$. The inequity in Theorem 1 is obtained as follows.
We first use the smoothness of the envelope function together with our
assumption that $\min_{x}\mathcal{F}_{\mu}(x)=0$ and
$x^{\star}\in\operatorname*{argmin}_{x}\mathcal{F}_{\mu}(x)$ to obtain

$\mathcal{E}(\psi_{1}(0))\,\leq\,\tfrac{\tilde{L}}{2}\|\psi_{1}(0)\|^{2}\,\leq\,\tfrac{\tilde{L}}{2}\|\psi(0)\|^{2}$

where $\tilde{L}$ is given in Lemma 2. Then, application of the triangle
inequality, $H\preceq I$, and the fact that $\theta(t)\leq 1$ for $t\geq 0$
yield the upper bound on the Lyapunov function,

$\displaystyle V(\psi(0))$
$\displaystyle\,=\,\tfrac{\alpha\tilde{L}}{2}\|\psi(0)\|^{2}\,+\,\tfrac{1}{2}\|\theta(t)\psi_{1}(0)\,+\,\psi_{2}(0)\|_{H}^{2}$
$\displaystyle\,\leq\,(\tfrac{\alpha\tilde{L}}{2}\,+\,1)\|\psi(0)\|^{2}.$ (34)

We note that for both (21) and (21),
$\mathcal{E}(\psi_{1})=\mathcal{F}_{\mu}(x)$. Hence, a lower bound on the
Lyapunov function can be obtained using the identity [21, Theorem 2.2],
$F(p_{\mu}(x))\leq\mathcal{F}_{\mu}(x)$ for any $x\in\mathbb{R}^{n}$. This
lower bound in conjunction with (34) leads to the inequality in Theorem 1 for
(21) with constant $c_{1}=1/\tilde{L}+1/\alpha$.

In the case of (21) for which $\psi=[z^{T}~{}\dot{z}^{T}]^{T}$, while the
lower bound on the Lyapunov function stays the same, the upper bound (34) can
be expressed in terms of $x$ and $\dot{x}$ using the following relations [20],

$\begin{array}[]{rcl}x&=&\operatorname{\mathbf{prox}}_{\mu
f}(z)\,=\,(I\,+\,\mu Q)^{-1}z,\\\\[2.84544pt]
\dot{x}&=&\nabla\operatorname{\mathbf{prox}}_{\mu f}(z)\dot{z}\,=\,(I\,+\,\mu
Q)^{-1}\dot{z}\end{array}$

which lead to the inequality in Theorem 1 with constant $c_{1}=(1+\mu
L)(1/\tilde{L}+1/\alpha)$.

####  4.1.2 Theorem 2: strongly convex quadratic $f$

For a strongly convex $f$, the use of constant $\alpha$, $\beta$, and $\gamma$
in (21) leads to time-invariant dynamics (22a). In this case, the Lyapunov
function candidate does not explicitly depend on time, i.e., $R$ in (26)
becomes a constant matrix. We again let $P=R^{T}HR$ and utilize (27) and (28)
to write $\Pi$ in (29) as

$\begin{array}[]{rcl}\Pi&=&A^{T}P\,+\,PA\,+\,\theta
P\,-\,\alpha\widetilde{m}\theta(C^{T}HC\,+\,C_{2}^{T}HC_{2}).\end{array}$

Since the off-diagonal blocks $\Gamma$ and $\Lambda$ in (29) remain the same
as in (32), we have

$\langle\psi,\Gamma u\,+\,\Lambda
v\rangle\,=\,-(\alpha/\beta)(1\,-\,\theta\beta)\langle
C_{2}\psi,u-v\rangle_{H}.$ (35)

Moreover, strong monotonicity of gradient of the envelope function along with
$H\preceq I$ implies $\langle
u-v,C_{2}\psi\rangle_{H}\geq\widetilde{m}\|C_{2}\psi\|_{H}^{2}$ which together
with (35) simplifies (29) to

$\dot{V}\,+\,\theta V\,\leq\,(1/2)\psi^{T}(\Pi\,-\,\Upsilon)\psi$ (36)

where

$\Upsilon\,=\,(2\alpha\widetilde{m}/\beta)(1\,-\,\theta\beta)C_{2}^{T}HC_{2}.$

Substitution of the system matrices in (22b) gives

$\Pi-\Upsilon=\left[\begin{array}[]{cc}\theta(\theta^{2}-\alpha\widetilde{m})&\theta(2\theta-\gamma-\alpha\widetilde{m}\beta)\\\
\theta(2\theta-\gamma-\alpha\widetilde{m}\beta)&3\theta-2\gamma-2\alpha\widetilde{m}\beta\end{array}\right]\otimes
H.$

The conditions that guarantee $\Pi-\Upsilon\preceq 0$ can be simplified by
setting off-diagonal blocks in $\Pi-\Upsilon$ to zero, i.e.,

$\theta\,=\,(\gamma\,+\,\alpha\widetilde{m}\beta)/2.$ (37)

When (37) is satisfied, the second diagonal block in $\Pi-\Upsilon$ is clearly
negative definite. Since $\beta=1-\gamma$, the negative semidefiniteness of
the first diagonal block requires
$\gamma\leq(2\sqrt{\alpha\widetilde{m}}-\alpha\widetilde{m})/(1-\alpha\widetilde{m})$
which is trivially satisfied by the choice of
$\gamma=(2\sqrt{\alpha\widetilde{m}}-2\alpha\widetilde{m})/(1-\alpha\widetilde{m})$
given in the theorem. Hence, the right-hand-side of (36) is negative
semidefinite, leading to $V(t)\leq V(0)\mathrm{e}^{-\theta t}$ for $t\geq 0$.

The strong convexity and smoothness of the envelope functions yield quadratic
lower and upper bounds on the Lyapunov function, which completes the proof for
(21). For (21), the proof is completed as in Section 4.1.1, i.e., utilizing
linearity of the proximal operator together with
$\dot{x}=\nabla\operatorname{\mathbf{prox}}_{\mu f}(z)\dot{z}$ and (12).

###  4.2 Theorem 3: strongly convex $f$

To prove our results for general strongly convex $f$, we use system
representation given by (22) with constant parameters but with a different
Lyapunov function and a different characterization of the FB envelope.

Lyapunov function. For general strongly convex $f$, since the weight matrix
$H(t)=I-\mu\nabla^{2}f(x(t))$ is time-dependent, in lieu of (26) we propose
the following time-invariant Lyapunov function candidate

$V(\psi)\,=\,\alpha\mathcal{F}_{\mu}(C\psi)\,+\,(1/2)\psi^{T}P\psi$ (38)

where $C$ is given in (22b), $P=R^{T}R$, and $R$ is given in (26). In contrast
to (26), the envelope function in (38) is evaluated at $y=C\psi$ and the
quadratic form does not use the time-dependent weight matrix $H$.

Characterization of the FB envelope. For general strongly convex problems, the
FB envelope is not necessarily convex and the upper bounds (24) and (25) no
longer hold. In Lemmas 3, we derive two novel inequalities that substitute
(24); see Appendix for the proof.

######  Lemma 3.

Let Assumption 1 hold and let $\mu\in(0,1/L)$. For every
$x,x^{\prime}\in\mathbb{R}^{n}$ and $x^{\star}\in\operatorname*{argmin}F(x)$,
FB envelope satisfies

$\displaystyle\begin{array}[]{l}\\!\\!\\!\\!\\!\mathcal{F}_{\mu}(x)\,-\,F(x^{\prime})\,\leq\\\\[-1.42271pt]
\\!\\!\\!\\!\\!\langle
G_{\mu}(x),x\,-\,x^{\prime}\rangle\,-\,\frac{m}{2}\|x\,-\,x^{\prime}\|_{2}^{2}\,-\,\frac{\mu}{2}\|G_{\mu}(x)\|_{2}^{2}.\end{array}$
(39c)
$\displaystyle\mathcal{F}_{\mu}(x)\,-\,F(x^{\star})\,\geq\,\tfrac{m^{2}(1\,-\,\mu
L)}{2L}\|x\,-\,x^{\star}\|_{2}^{2}.$ (39d)

We note that in contrast to strong convexity inequality (24), in Lemma 3, the
gradient of $\mathcal{F}_{\mu}$ is not involved and we characterize the
difference between envelope $\mathcal{F}_{\mu}$ and objective function $F$.
Still, Lyapunov function (38) can be upper bounded via Lemma 2 based on the
equivalence given in (14).

An upper bound on the Lyapunov function. The following upper bound on (38) can
be directly obtained using Lemma 3, (14), and system representation (22a),

$\begin{array}[]{rcl}\theta V&\leq&\alpha\theta\langle
u,C\psi\rangle-\frac{\alpha
m\theta}{2}\|C\psi\|_{2}^{2}-\frac{\alpha\mu\theta}{2}\|u\|_{2}^{2}+\tfrac{\theta}{2}\psi^{T}P\psi.\end{array}$
(40)

Lyapunov analysis. The derivative of (38) along the solutions of (21) is
determined by

$\dot{V}\,=\,\alpha\langle Hu,CA\psi\,+\,CBu\rangle\,+\,\langle
P\psi,A\psi\,+\,Bu\rangle.$

Summing $\dot{V}$ with the upper bound given in (40) yields

$\dot{V}\,+\,\theta V\,\leq\,\frac{1}{2}\left[\begin{array}[]{c}\psi\\\
u\end{array}\right]^{T}\left[\begin{array}[]{cc}\Pi&\Gamma+\Lambda H\\\
\Gamma^{T}+H\Lambda^{T}&-2\alpha^{2}\beta
H\end{array}\right]\left[\begin{array}[]{c}\psi\\\ u\end{array}\right]$

where

$\begin{array}[]{rcl}\Pi&=&A^{T}P\,+\,PA\,+\,\theta P\,-\,\alpha m\theta
C^{T}C\\\ \Gamma&=&PB\,+\,\alpha\theta C^{T}\\\ \Lambda&=&\alpha
A^{T}C^{T}.\end{array}$ (41)

Since $H\succ 0$, the Schur complement implies that $\dot{V}+\theta V\leq 0$
if and only if

$2\alpha^{2}\beta\Pi+(\Gamma\Lambda^{T}+\Lambda\Gamma^{T})+\Gamma
H^{-1}\Gamma^{T}+\Lambda H\Lambda^{T}\preceq 0.$ (42)

Substitution of system matrices (22b) into (41) gives

$\displaystyle\Gamma\Lambda^{T}\,+\,\Lambda\Gamma^{T}\,=\,\left[\begin{array}[]{cc}0&0\\\
0&-2\alpha^{2}(1\,-\,\theta\beta)(1\,-\,\gamma\beta)I\end{array}\right]$ (43c)
$\displaystyle\Gamma H^{-1}\Gamma^{T}\,+\,\Lambda H\Lambda^{T}\,=\,$
$\displaystyle\left[\begin{array}[]{cc}0&~{}0\\\
0&~{}\alpha^{2}\\!\left((1-\theta\beta)^{2}H^{-1}+(1-\
\gamma\beta)^{2}H\right)\end{array}\right].$ (43f)

The eigenvalue decomposition of $H$ can be used to show

$\displaystyle\Gamma H^{-1}\Gamma^{T}\,+\,\Lambda
H\Lambda^{T}\,\preceq\,\left[\begin{array}[]{cc}0&0\\\ 0&\alpha^{2}\eta
I\end{array}\right]$ (44c)
$\displaystyle\\!\\!\eta\mathrel{\mathop{:}}=\operatorname*{maximize}\limits_{\sigma\,\in\,[m,L]}~{}\dfrac{(1-\theta\beta)^{2}}{1-\mu\sigma}\,+\,(1-\gamma\beta)^{2}(1-\mu\sigma).$
(44d)

If $(1-\theta\beta)\geq(1-\mu\sigma)(1-\gamma\beta)$ for all $\sigma\in[m,L]$,
then the objective function in (44d) is increasing for $\sigma\in[m,L]$ and
the maximum occurs at $\sigma=L$,

$\eta\,=\,(1\,-\,\theta\beta)^{2}/(1\,-\,\mu
L)\,+\,(1\,-\,\gamma\beta)^{2}(1\,-\,\mu L).$ (45)

Combining (43c) with (44c) and (45) yields

$\begin{array}[]{l}(\Gamma\Lambda^{T}\,+\,\Lambda\Gamma^{T})\,+\,\Gamma
H^{-1}\Gamma^{T}\,+\,\Lambda H\Lambda^{T}\,\preceq\\\
\left[\begin{array}[]{cc}0&0\\\ 0&\frac{\alpha^{2}}{1-\mu
L}\big{(}(1-\theta\beta)-(1-\mu
L)(1-\gamma\beta)\big{)}^{2}I\end{array}\right].\end{array}$

Moreover, since the entries of the matrix $\Pi$ are given by

$\Pi\,=\,\left[\begin{array}[]{cc}\theta(\theta^{2}\,-\,\alpha
m)&\theta(2\theta\,-\,\gamma\,-\,\alpha m\beta)\\\
\theta(2\theta\,-\,\gamma\,-\,\alpha m\beta)&3\theta\,-\,2\gamma\,-\,\alpha
m\theta\beta^{2}\end{array}\right]\otimes I$

setting the off-diagonal blocks in $\Pi$ to zero yields

$\theta\,=\,(\gamma\,+\,\alpha m\beta)/2.$ (46)

When (46) is satisfied, verifying $\dot{V}+\theta V\leq 0$, or equivalently
(42), amounts to checking these three conditions,

  1. (i)

$(1-\theta\beta)\geq(1-\mu\sigma)(1-\gamma\beta)$ for all $\sigma\in[m,L]$

  2. (ii)

$\theta^{2}\leq\alpha m$

  3. (iii)

$\frac{((1-\theta\beta)-(1-\mu L)(1-\gamma\beta))^{2}}{2\beta(1-\mu
L)}\,\leq\,\alpha m\theta\beta^{2}+2\gamma-3\theta$.

While condition (i) is necessary for upper bound (45), conditions (ii) and
(iii) ensure negative semidefiniteness of the diagonal blocks in (42).

For $\gamma$ given in Theorem 3 and $\beta=1-\gamma$, we obtain
$\theta=\sqrt{\alpha m}-\alpha m/2$ using (46). Thus, $\theta\leq\gamma$ and
conditions (i) and (ii) are trivially satisfied. For large values of
$\kappa=L/m$, the numerator of the left-hand-side in condition (iii) is
dominated by $\mu L$. On the other hand, the largest term in the product of
the denominator and the right-hand-side is $\gamma\beta$, which suggests that
condition (iii) holds if $\mu\sim O(\sqrt{\gamma\beta}/L)$. Due to space
limitations, we use computational tools to verify this hypothesis. With the
above choices of $(\gamma,\beta,\theta)$ and $\mu L=\sqrt{\gamma\beta}/2$,
condition (iii) can be written as $wh(w)\leq 0$ where $w=\sqrt{\alpha
m}\in[0,1]$ and $h$: $[0,1]\to\mathbb{R}$. In Fig. 2, we show that $h$ is
negative for all $w\in[0,1]$. Moreover, since (i) implies that $wh(w)$ is an
increasing function of $\mu L$, condition (iii) is satisfied for all $\mu
L\in(0,\sqrt{\gamma\beta}/2]$.

|  $h(w)$
---
|
---
$w$
Figure 1: Plot of the function $h$ that demonstrates $h<0$. Thus, $wh(w)\leq
0$ for $w\in[0,1]$ and condition (iii) is satisfied.

As a result, using Grönwall inequality, we conclude that $V(t)\leq
V(0)\mathrm{e}^{-\theta t}$ for $t\geq 0$ along the solutions of (21). A
quadratic upper bound on the Lyapunov function can be obtain by using the
smoothness of the FB envelope as in the proof of Theorem 1. Moreover, since
$C^{T}C+P$ is invertible, inequality (39d) in Lemma 3 yields a quadratic lower
bound on $V(t)$. These upper and lower quadratic bounds completes the proof.

##  5 Computational experiments

In our computational experiments we choose zero initial conditions and
$\alpha=2\mu=1/L$ unless stated otherwise.

###  5.1 Example 1: $\ell_{1}$-regularized least squares

Let us examine an $\ell_{1}$-regularized least squares problem,

$\operatorname*{minimize}\limits_{x}~{}(1/2)\|Ex\,-\,q\|_{2}^{2}\,+\,\lambda\|x\|_{1}$
(47)

where the problem data $E\in\mathbb{R}^{100\times 2000}$,
$q\in\mathbb{R}^{100}$, and $\lambda>0$ are generated as described in [18].
Since $E^{T}E$ is a singular matrix, this problem is not strongly convex.
Thus, the time-varying parameters given in Theorem 1 are used to illustrate
improved sublinear convergence rate of accelerated proximal gradient flow and
DR splitting dynamics (21) and (21) in Fig. 2.

###  5.2 Example 2: Box-constrained quadratic program

For a box-constrained quadratic program,

$\begin{array}[]{l}\operatorname*{minimize}\limits_{x}~{}(1/2)x^{T}Qx\,+\,q^{T}x\,+\,I_{\mathcal{C}}(x)\end{array}$
(48)

where $I_{\mathcal{C}}(x)$ is the indicator function of the convex set
$\mathcal{C}=\\{x~{}|~{}l\leq x\leq u\\}$, we generate problem data with
positive definite $Q\in\mathbb{R}^{500\times 500}$, and vectors
$q,l,u\in\mathbb{R}^{500}$ as described in [9] and set the condition number
$\kappa=10^{5}$. Linear convergence of the underlying algorithms for constant
parameters given in Theorem 2 is illustrated in Fig. 2.

###  5.3 Example 3: $\ell_{1}$-regularized logistic regression

We examine the $\ell_{1}$-regularized logistic regression problem

$\operatorname*{minimize}\limits_{x}\sum\limits_{i\,=\,1}^{n}\\!\\!\left(-y_{i}a_{i}^{T}x+\log(1+\mathrm{e}^{a_{i}^{T}x})\right)+\lambda\|x\|_{1}+\tfrac{\varrho}{2}\|x\|_{2}^{2}$
(49)

where the problem data $A=[a_{1}\;\cdots~{}a_{n}]^{T}\in\mathbb{R}^{200\times
1000}$ as well as the scalars $\lambda$ and $y_{i}$ are generated as described
in [24]. In Fig. 2, we validate linear convergence of accelerated proximal
gradient flow dynamics (21) with constant parameters given in Theorem 3 and
$\mu=1/(L\sqrt[4]{\kappa})\approx 1/(26L)$ for this strongly convex problem
with $\varrho=0.1$, and the condition number $\kappa\approx 5\cdot 10^{5}$.

|  |  |  |  |
---|---|---|---|---|---
|  $F(x(t))-F(x^{\star})$
---
|
---
$t$
|  $\|x(t)\,-\,x^{\star}\|^{2}_{2}/\|x^{\star}\|^{2}_{2}$
---
|
---
$t$
|  $\|x(t)\,-\,x^{\star}\|^{2}_{2}/\|x^{\star}\|^{2}_{2}$
---
|
---
$t$

Figure 2: (a) Sublinear convergence of accelerated dynamics (21) and (21) with
the time-varying parameters given in Theorem 1 for $\ell_{1}$-regularized
least squares problem (47). (b) Linear convergence of accelerated dynamics
(21) and (21) with the constant parameters given in Theorem 2 for box-
constrained quadratic program (48). (c) Linear convergence of accelerated
dynamics (21) eith the constant parameters given in Theorem 3 and
$\mu=1/(L\sqrt[4]{\kappa})\approx 1/(26L)$ for the $\ell_{1}$-regularized
logistic regression problem (49) with $\varrho=0.1$, and condition number
$\kappa\approx 5\cdot 10^{5}$.

##  6 Concluding remarks

We have introduced continuous-time accelerated FB and DR splitting dynamics
for composite optimization problems where the objective function is given by
the sum of smooth and nonsmooth terms. Motivated by the variable-metric
gradient method interpretation of FB and DR splittings, we have leveraged
smoothness of associated envelopes and utilized a Lyapunov-based approach to
establish convergence rates for both convex and strongly convex problems.

##  Appendix A Proof of Lemma 3

Our proof is based on an equivalent definition of the FB envelope [21, Eq.
(2.8b)]

$\\!\mathcal{F}_{\mu}(x)=f(x)+g(p(x))-\mu\langle\nabla
f(x),G_{\mu}(x)\rangle+\tfrac{\mu}{2}\|G_{\mu}(x)\|_{2}^{2}$ (50)

and the subgradient inequality,

$g(x)\geq g(x^{\prime})+\langle r,x-x^{\prime}\rangle,~{}\forall
x,x^{\prime}\in\mathbb{R}^{n},~{}\forall r\in\partial g(x^{\prime}).$ (51)

Proof of (39c). Let
$p(x)\mathrel{\mathop{:}}=\operatorname{\mathbf{prox}}_{\mu g}(x-\mu\nabla
f(x))$. The identity (8) together with (10) implies that

$G_{\mu}(x)-\nabla f(x)\in\partial g(p(x)).$

Hence, the subgradient inequality (51) for an arbitrary
$x^{\prime}\in\mathbb{R}^{n}$ yields

$g(x^{\prime})\;\geq\;g(p(x))\,+\,\langle G(x)\,-\,\nabla
f(x),x^{\prime}\,-\,p(x)\rangle$

which in conjunction with (50) results in

$\begin{array}[]{rcl}\mathcal{F}_{\mu}(x)&\leq&f(x)\,+\,g(x^{\prime})\,-\,\langle
G(x)\,-\,\nabla f(x),x^{\prime}\,-\,p(x)\rangle\\\ &&\,-\,\mu\langle\nabla
f(x),G_{\mu}(x)\rangle\,+\,\frac{\mu}{2}\|G_{\mu}(x)\|_{2}^{2}.\end{array}$

Rearranging terms and then upper bounding $f(x)$ via strong convexity,

$f(x)\,\leq\,f(x^{\prime})\,+\,\langle\nabla
f(x),x\,-\,x^{\prime}\rangle\,-\,\tfrac{m}{2}\|x\,-\,x^{\prime}\|^{2}$

gives

$\begin{array}[]{rcl}\mathcal{F}_{\mu}(x)&\leq&f(x^{\prime})\,+\,g(x^{\prime})\,-\,\frac{m}{2}\|x\,-\,x^{\prime}\|_{2}^{2}\\\
&&+\,\langle
G(x),p(x)\,-\,x^{\prime}\rangle\,+\,\frac{\mu}{2}\|G_{\mu}(x)\|_{2}^{2}\\\
&=&F(x^{\prime})\,-\,\frac{m}{2}\|x\,-\,x^{\prime}\|_{2}^{2}\,+\,\langle
G_{\mu}(x),x\,-\,x^{\prime}\rangle\\\ &&\,-\,\langle
G_{\mu}(x),x\,-\,p(x)\rangle\,+\,\frac{\mu}{2}\|G_{\mu}(x)\|_{2}^{2}\end{array}$

where the equality is obtained by adding and subtracting $x$ from the second
argument of the inner product. Substituting (10) into the second inner product
completes the proof.

Proof of (39d). Let
$p(x)\mathrel{\mathop{:}}=\operatorname{\mathbf{prox}}_{\mu g}(x-\mu\nabla
f(x))$, $\tilde{x}\mathrel{\mathop{:}}=x-x^{\star}$, and
$\tilde{p}\mathrel{\mathop{:}}=p(x)-p(x^{\star})$. Substitution of (10) into
(50) gives

$\begin{array}[]{rcl}\mathcal{F}_{\mu}(x)\,-\,F(x^{\star})&=&f(x)\,-\,f(x^{\star})\,+\,g(p(x))\,-\,g(x^{\star})\\\
&&\,-\,\langle\nabla
f(x),\tilde{x}\,-\,\tilde{p}\rangle\,+\,\frac{\mu}{2}\|G_{\mu}(x)\|^{2}.\end{array}$
(52)

The optimality condition for problem (6) implies that $-\nabla
f(x^{\star})\in\partial g(x^{\star})$ and $x^{\star}=p(x^{\star})$. Using
these two identities in the subgradient inequality (51) yields

$g(p(x))\,\geq\,g(x^{\star})\,-\,\langle\nabla f(x^{\star}),\tilde{p}\rangle$

which together with (52) results in

$\begin{array}[]{rcl}\mathcal{F}_{\mu}(x)\,-\,F(x^{\star})&\geq&f(x)\,-\,f(x^{\star})\,-\,\langle\nabla
f(x),\tilde{x}\rangle\\\ &&\\!\\!\\!\\!+\langle\nabla f(x)\,-\,\nabla
f(x^{\star}),\tilde{p}\rangle\,+\,\frac{\mu}{2}\|G_{\mu}(x)\|^{2}.\end{array}$

Substituion of the lower bound given by the cocoercivity of gradient $\nabla
f(x)$,

$f(x)\,-\,f(x^{\star})\;\geq\;\langle\nabla
f(x^{\star}),\tilde{x}\rangle\,+\,\tfrac{1}{2L}\|\nabla f(x)\,-\,\nabla
f(x^{\star})\|_{2}^{2}$

into the right-hand-side yields

$\begin{array}[]{rcl}\mathcal{F}_{\mu}(x)\,-\,F(x^{\star})&\geq&\langle\nabla
f(x^{\star})\,-\,\nabla f(x),\tilde{x}\,-\,\tilde{p}\rangle\\\
&&\\!\\!\\!\\!+\tfrac{1}{2L}\|\nabla f(x)\,-\,\nabla
f(x^{\star})\|_{2}^{2}\,+\,\frac{\mu}{2}\|G_{\mu}(x)\|^{2}\\\
&=&\tfrac{\mu}{2}\|G_{\mu}(x)-(\nabla f(x)\,-\,\nabla f(x^{\star}))\|^{2}\\\
&&+\,\tfrac{1}{2L}(1\,-\,\mu L)\|\nabla f(x)\,-\,\nabla f(x^{\star})\|^{2}\\\
&\geq&\tfrac{m^{2}}{2L}(1\,-\,\mu L)\|x\,-\,x^{\star}\|_{2}^{2}\end{array}$

where the equality follows from (10) and the completion of squares; the second
inequality is obtained by ignoring the first (nonnegative) term and using
$m$-strong convexity of $f$ on the second one.

## References

  * [1] A. Allibhoy and J. Cortés.  Control barrier function-based design of gradient flows for constrained nonlinear programming.  IEEE Trans. Automat. Control, 2023.
  * [2] H. Attouch and Z. Chbani.  Fast inertial dynamics and FISTA algorithms in convex optimization. Perturbation aspects.  arXiv preprint arXiv:1507.01367, 2015.
  * [3] H. Attouch, Z. Chbani, J. Fadili, and H. Riahi.  Fast convergence of dynamical ADMM via time scaling of damped inertial dynamics.  J. Optim. Theory Appl., 193(1-3):704–736, 2022.
  * [4] H. Attouch, Z. Chbani, J. Peypouquet, and P. Redont.  Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity.  Math. Program., 168:123–175, 2018.
  * [5] P. Cisneros-Velarde and F. Bullo.  A contraction theory approach to optimization algorithms from acceleration flows.  In Proc. Int. Conf. Artif. Intell. Stat., pages 1321–1335, 2022\.
  * [6] N. K. Dhingra, S. Z. Khong, and M. R. Jovanović.  The proximal augmented Lagrangian method for nonsmooth composite optimization.  IEEE Trans. Automat. Control, 64(7):2861–2868, 2019.
  * [7] G. Franca, D. P. Robinson, and R. Vidal.  Gradient flows and proximal splitting methods: A unified view on accelerated and stochastic optimization.  Phys. Rev. E, 103(5):53304, 2021.
  * [8] P. Giselsson and M. Falt.  Envelope functions: Unifications and further properties.  J. Optim. Theory Appl., 178:673–698, 2018.
  * [9] C. C. Gonzaga, E. W. Karas, and D. R. Rossetto.  An optimal algorithm for constrained differentiable convex optimization.  SIAM J. Control Optim., 23(4):1939–1955, 2013.
  * [10] S. Hassan-Moghaddam and M. R. Jovanović.  Proximal gradient flow and Douglas-Rachford splitting dynamics: Global exponential stability via integral quadratic constraints.  Automatica J. IFAC, 123:109311, 2021.
  * [11] X. He, R. Hu, and Y. Fang.  “Second-order primal”+“first-order dual” dynamical systems with time scaling for linear equality constrained convex optimization problems.  IEEE Trans. Automat. Control, 67(8):4377–4383, 2022.
  * [12] T. Liu and T. K. Pong.  Further properties of the forward–backward envelope with applications to difference-of-convex programming.  Comput. Optim. Appl., 67:489–520, 2017.
  * [13] M. Muehlebach and M. I. Jordan.  A dynamical systems perspective on Nesterov acceleration.  In Proc. Int. Conf. Mach. Learn., pages 4656–4662, 2019.
  * [14] Y. Nesterov.  A method of solving a convex programming problem with convergence rate ${O}(1/k^{2})$.  In Dokl. Akad. Nauk, volume 269, pages 543–547, 1983.
  * [15] Y. Nesterov.  Introductory lectures on convex optimization: A basic course, volume 87.  Springer, Berlin, 2003.
  * [16] I. K. Ozaslan and M. R. Jovanović.  Exponential convergence of primal-dual dynamics for multi-block problems under local error bound condition.  In Proc. IEEE Conf. Decis. Control, pages 7579–7584, 2022.
  * [17] I. K. Ozaslan and M. R. Jovanović.  Tight lower bounds on the convergence rate of primal-dual dynamics for equality constrained convex problems.  In Proc. IEEE Conf. Decis. Control, pages 7318–7323, 2023.
  * [18] B. O’Donoghue and E. Candes.  Adaptive restart for accelerated gradient schemes.  Found. Comput. Math., 15:715–732, 2015.
  * [19] N. Parikh and S. Boyd.  Proximal algorithms.  Found. Trends Opt., 1(3):123–231, 2014.
  * [20] P. Patrinos, L. Stella, and A. Bemporad.  Douglas-Rachford splitting: Complexity estimates and accelerated variants.  In Proc. IEEE Conf. Decision Control, pages 4234–4239, 2014.
  * [21] P. Patrinos, L. Stella, and A. Bemporad.  Forward-backward truncated Newton methods for convex composite optimization.  arXiv:1402.6655, 2014.
  * [22] B. T. Polyak.  Some methods of speeding up the convergence of iteration methods.  Comput. Math. Phys., 4(5):1–17, 1964.
  * [23] E. K Ryu and S. Boyd.  Primer on monotone operator methods.  Appl. comput. math, 15(1):3–43, 2016.
  * [24] W. Su, S. Boyd, and E. Candes.  A differential equation for modeling Nesterov’s accelerated gradient method: theory and insights.  J. Mach. Learn. Res., 17(153):1–43, 2016.
  * [25] A. Wibisono, A. C. Wilson, and M. I. Jordan.  A variational perspective on accelerated methods in optimization.  Proc. Natl. Acad. Sci. USA, 113(47):E7351–E7358, 2016.
  * [26] A. C Wilson, B. Recht, and M. I. Jordan.  A Lyapunov analysis of accelerated methods in optimization.  J. Mach. Learn. Res., 22(1):5040–5073, 2021.
  * [27] X. Zeng, J. Lei, and J. Chen.  Dynamical primal-dual Nesterov accelerated method and its application to network optimization.  IEEE Trans. Automat. Control, 68(3):1760–1767, 2023.