A Control-Theoretic Approach to Analysis and Parameter Selection of Douglas-Rachford Splitting

TL;DR

This paper introduces a control-theoretic framework for analyzing and selecting parameters in Douglas-Rachford splitting and ADMM algorithms, providing convergence guarantees and simplifying rate proofs.

Contribution

It offers a unified, dynamical system-based analysis method for convergence and parameter tuning of Douglas-Rachford splitting, applicable under various convexity and smoothness conditions.

Findings

Derived a dimensionally independent matrix inequality for convergence analysis.

Provided performance guarantees and optimal parameter settings.

Simplified proof of O(1/k) convergence rate without strong convexity.

Abstract

Douglas-Rachford splitting and its equivalent dual formulation ADMM are widely used iterative methods in composite optimization problems arising in control and machine learning applications. The performance of these algorithms depends on the choice of step size parameters, for which the optimal values are known in some specific cases, and otherwise are set heuristically. We provide a new unified method of convergence analysis and parameter selection by interpreting the algorithm as a linear dynamical system with nonlinear feedback. This approach allows us to derive a dimensionally independent matrix inequality whose feasibility is sufficient for the algorithm to converge at a specified rate. By analyzing this inequality, we are able to give performance guarantees and parameter settings of the algorithm under a variety of assumptions regarding the convexity and smoothness of the objective function. In particular, our framework enables us to obtain a new and simple proof of the O(1/k) convergence rate of the algorithm when the objective function is not strongly convex.

Equations

$${\mathrm{minimize}}_{x\in {\mathbb{R}}^d}\{F(x)=f(x)+g(x)\},$$

$$y_k$$

$$z_k$$

$$x_{k+1}$$

$${\mathrm{prox}}_{\alpha f}(x)={\text{argmin}}_y\left\{f(y)+\frac{1}{2\alpha }\Vert x-y{\Vert }_2^2\right\}.$$

$$y=x-\alpha \partial f(y).$$

$${\left[\begin{array}{c}x-y\\ \varphi (x)-\varphi (y)\end{array}\right]}^{\top }Q\left[\begin{array}{c}x-y\\ \varphi (x)-\varphi (y)\end{array}\right]\ge 0.$$

$$Q(m,L)=\left[\begin{array}{cc}-\frac{mL}{m+L}& 1/2\\ 1/2& -\frac{1}{m+L}\end{array}\right]\otimes I_d,$$

$$Q_p(m,L,\alpha )=\left[\begin{array}{cc}0& I_d\\ \alpha I_d& -I_d\end{array}\right]Q(m,L)\left[\begin{array}{cc}0& \alpha I_d\\ I_d& -I_d\end{array}\right],$$

$$x_{k+1}=x_k+{\lambda }_k\varphi (x_k),$$

$$\varphi (x_k):={\mathrm{prox}}_{\alpha g}(2{\mathrm{prox}}_{\alpha f}(x_k)-x_k)-{\mathrm{prox}}_{\alpha g}(x_k).$$

$$y_k$$

$$z_k$$

$$x_{k+1}$$

$$x_{\star }=y_{\star }+\alpha \partial f(y_{\star }),y_{\star }=z_{\star },\partial f(y_{\star })+\partial g(z_{\star })=0.$$

$$\varphi (x_k)=z_k-y_k=-\alpha (\partial f(y_k)+\partial g(z_k)).$$

$$V_k=\Vert x_k-x_{\star }{\Vert }_2^2+\sum _{i=0}^{k-1}{\theta }_i\Vert \partial f(y_i)+\partial g(z_i){\Vert }_2^2,$$

$$\underset{i=0,\dots ,k-1}{\min }\Vert \partial f(y_i)+\partial g(z_i){\Vert }_2^2\le \frac{1}{{\Theta }_k}\Vert x_0-x_{\star }{\Vert }_2^2.$$

$$\Vert x_k-x_{\star }{\Vert }_2^2+\sum _{i=0}^{k-1}{\theta }_i\Vert \partial f(y_i)+\partial g(z_i){\Vert }_2^2\le \Vert x_0-x_{\star }{\Vert }_2^2.$$

$$\sum _{i=0}^{k-1}\frac{{\theta }_i}{{\Theta }_k}\Vert \partial f(y_i)+\partial g(z_i){\Vert }_2^2\le \frac{\Vert x_0-x_{\star }{\Vert }_2^2}{{\Theta }_k}.$$

$$W_k^{(0)}+{\sigma }_k^{(1)}{Q}^{(1)}+{\sigma }_k^{(2)}{Q}^{(2)}⪯0,$$

$$W_k^{(0)}$$

$$Q^{(1)}$$

$$Q^{(2)}$$

$$\underset{i=0,\dots ,k-1}{\min }\Vert \partial f(y_i)+\partial g(z_i){\Vert }_2^2\le \frac{1}{{\Theta }_k}\Vert x_0-x_{\star }{\Vert }_2^2.$$

$$e_k:={\left[\begin{array}{ccc}(x_k-x_{\star }{)}^{\top }& (y_k-y_{\star }{)}^{\top }& (z_k-z_{\star }{)}^{\top }\end{array}\right]}^{\top }.$$

$$V_{k+1}-V_k=e_k^{\top }{W}_k^{(0)}{e}_k,$$

$$e_k^{\top }{Q}^{(1)}{e}_k={\left[\begin{array}{c}{x}_k-x_{\star }\\ y_k-y_{\star }\end{array}\right]}^{\top }{Q}_p(m_f,L_f,\alpha )\left[\begin{array}{c}{x}_k-x_{\star }\\ y_k-y_{\star }\end{array}\right],$$

$$\left[\begin{array}{ccc}0& 0& \alpha I_d\\ -I_d& 2I_d& -I_d\end{array}\right]e_k=\left[\begin{array}{cc}0& \alpha I_d\\ I_d& -I_d\end{array}\right]\left[\begin{array}{c}(2y_k-x_k)-(2y_{\star }-x_{\star })\\ z_k-z_{\star }\end{array}\right].$$

$$e_k^{\top }{W}_k^{(0)}{e}_k+{\sigma }_k^{(1)}{e}_k^{\top }{Q}^{(1)}{e}_k+{\sigma }_k^{(2)}{e}_k^{\top }{Q}^{(2)}{e}_k\le 0.$$

$$V_k=\Vert x_k-x_{\star }{\Vert }_2^2+\sum _{i=0}^{k-1}{\theta }_i[F(z_i)-F(z_{\star })],$$

$$\underset{i=0,\dots ,k-1}{\min }[F(z_i)-F(z_{\star })]\le \frac{1}{{\Theta }_k}\Vert x_0-x_{\star }{\Vert }_2^2.$$

$$W_k^{(1)}+{\sigma }_k^{(1)}{Q}^{(1)}+{\sigma }_k^{(2)}{Q}^{(2)}⪯0,$$

Full text

I Introduction

In this paper, we consider problems of the form

[TABLE]

where $f,g:\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\}$ are convex, closed, and proper (c.c.p.). Douglas-Rachford splitting (DRS) solves problem (1) with the following iterations:

[TABLE]

where $\mathrm{prox}$ is the proximal operator (see Definition 1) and $\alpha$ and $\lambda_{k}$ are known as the proximal step size and relaxation parameter, respectively. For a proper selection of these parameters, the limiting values of both $y_{k}$ and $z_{k}$ will be a solution to (1). The goal of this work is to provide convergence rates for DRS over various assumptions on $f$ and $g$, and optimize these rates with respect to the algorithm parameters $\alpha$ and $\lambda_{k}$ using semidefinite programs (SDPs).

The algorithm was first proposed in [1], and has since found application in general separable optimization problems [2]. Its dual formulation, ADMM, has been particularly useful in distributed optimization problems [3]. Since the iterates of ADMM can be written as applying DRS to the dual problem [4, 5], convergence results for one algorithm are valid for the other as well when strong duality holds.

The convergence of DRS has previously been analyzed using monotone operator theory and variational inequalities, see [6], [7]. These techniques have led to proofs of a $O(1/k)$ convergence rate for the non-strongly convex case [8], and linear convergence when $f$ is smooth and strongly convex [9], [10]. In the most general case, a condition for convergence is that $\lambda_{k}\in(0,2)$ with $\sum_{k=0}^{\infty}\lambda_{k}(2-\lambda_{k})=\infty$ [11], though there exist cases where the algorithm converges with $\lambda_{k}>2$ [12].

Recently, there has been interest in automating the analysis and design of optimization algorithms via SDPs, [13, 14, 15, 16, 17, 18, 19]. In particular, through the method of integral quadratic constraints proposed in [14], the authors of [12] derive an SDP for choosing the parameters of ADMM in the case of smooth and strongly convex $f$. Using a similar framework, the authors of [20] provide evidence that as the relaxation parameter approaches $2$ from below, the linear convergence rate is close to being optimal and in [21] are able to analytically solve the SDP to give a convergence rate. The work in [22] gives an optimal choice for the relaxation parameter when $f$ is quadratic. Furthermore, [23] gives a set of assumptions in which a bound on the linear convergence rate is minimized by setting $\lambda_{k}=2$.

Our Contribution: By viewing DRS as a linear system with non-linear feedback, we derive a dimensionally independent matrix inequality which gives convergence guarantees via Lyapunov functions. Whereas such an approach was previously applied in [12] to the case of smooth and strongly convex $f$, our framework is novel in that it encompasses varying assumptions on the smoothness and convexity of $f$. By changing a single term in the Lyapunov function for each scenario, we are able to relate the satisfaction of a matrix inequality to the convergence of the algorithm. In particular, we give a new and simple proof of $O(1/k)$ convergence in the non-strongly convex case. These symbolic results can then be used to select step sizes that optimize the derived rates.

In the strongly convex case, the corresponding matrix inequality is sufficient to guarantee a linear convergence rate. We are able to modify the matrix inequality to linearize the dependence on $\lambda_{k}$, allowing us to numerically optimize its value for the convergence rate directly. While previous work derived SDP’s which can verify the performance of the algorithm for a given parameter setting, to the best of our knowledge this is the first time such a method immediately gives an optimal relaxation parameter when solved numerically, as opposed to having to search over a range of values for $\lambda_{k}$.

II Preliminaries

We denote the set of real numbers by $\mathbb{R}$, the set of real $n$-dimensional vectors by $\mathbb{R}^{n}$, the set of real $m\times n$-dimensional matrices by $\mathbb{R}^{m\times n}$, and the $n$-dimensional identity matrix and zero matrix by $I_{n}$ and $0_{n}$, respectively. For a function $f\colon\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\}$, we denote by $\mathrm{dom}\,f=\{x\in\mathbb{R}^{d}\colon f(x)<\infty\}$ the effective domain of $f$. The subdifferential of a function $f$ at a point $x$ is $\partial f(x):=\{g\;|\;f(y)-f(x)\geq g^{\top}(y-x),\forall y\in\mathrm{dom}(f)\}$. By abuse of notation we will also refer to a subgradient, that is an element of the subdifferential by $\partial f(x)$ as well. The indicator function of a set $C$ is given by $\mathrm{1}_{C}(x)=0$ if $x\in C$ and $\mathrm{1}_{C}(x)=\infty$ if $x\notin C$. For two matrices $A\in\mathbb{R}^{m\times n}$ and $B\in\mathbb{R}^{p\times q}$ their Kronecker product is $A\otimes B$.

We say a differentiable function $f\colon\mathbb{R}^{d}\to\mathbb{R}$ is $L_{f}$-smooth on $\mathcal{S}\subseteq\mathrm{dom}\,f$ if $\|\nabla f(x)-\nabla f(y)\|_{2}\leq L_{f}\|x-y\|_{2}$ for some $L_{f}>0$ and all $x,y\in\mathcal{S}$. This also implies for all $x,y\in\mathcal{S}$, $f(y)\leq f(x)+\nabla f(x)^{\top}(y-x)+(L_{f}/2)\|y-x\|_{2}^{2}$. A differentiable function $f\colon\mathbb{R}^{d}\to\mathbb{R}$ is $m_{f}$-strongly convex on $\mathcal{S}\subseteq\mathrm{dom}\,f$ if $m_{f}\|x-y\|_{2}^{2}\leq(x-y)^{\top}(\nabla f(x)-\nabla f(y))$ for some $m_{f}>0$ and all $x,y\in\mathcal{S}$. The class of functions which are $L_{f}$-smooth and $m_{f}$-strongly convex is denoted by $\mathcal{F}(m_{f},L_{f})$.

Definition 1** **(Proximal Operator)

Given a c.c.p. function $f:\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\}$ and $\alpha>0$, the proximal operator $\mathrm{prox}_{\alpha f}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is defined as

[TABLE]

The point $y=\mathrm{prox}_{\alpha f}(x)$ also is given by the implicit solution to the subgradient equation

[TABLE]

We say that a nonlinear function $\phi:\mathbb{R}^{d}\to\mathbb{R}^{d}$ satisfies the incremental quadratic constraint [24] (or point-wise integral quadratic constraint [14]) defined by $Q\in\mathbb{R}^{2d\times 2d}$ if for all $x,y$,

[TABLE]

For $0\leq m,L<\infty$ define

[TABLE]

and define $Q(m,\infty)$ as $\lim_{L\to\infty}Q(m,L)$. It was noted in [25, 14] that a differentiable function $f$ belongs to the class $\mathcal{F}(m_{f},L_{f})$ on $\mathcal{S}$ if and only if the gradient $\nabla f$ satisfies the incremental QC defined by $Q(m_{f},L_{f})$. When $L_{f}=\infty$, the subgradient $\partial f$ satisfies the QC defined by $Q(m_{f},\infty)$. If we define

[TABLE]

then the proximal operator of a function $f\in\mathcal{F}(m,L)$, $\mathrm{prox}_{\alpha f}$, satisfies the incremental QC defined by $Q_{p}(m,L,\alpha)$ [15].

III Analysis of Douglas-Rachford Splitting via Matrix Inequalities

III-A Douglas-Rachford Splitting as a Dynamical System

We can write the updates in (2) as a linear system with state $x_{k}$ and feedback nonlinearity $\phi(x_{k})$,

[TABLE]

where

[TABLE]

Our main technique is to describe the nonlinearity $\phi$ with incremental QCs representing the $\mathrm{prox}$ operator. This allows us to derive a matrix inequality as a sufficient condition for closed-loop stability of the system via a Lyapunov function argument. We perform this derivation in the following three cases:

• •

Case 1: $f\in\mathcal{F}(0,\infty)$ and $g\in\mathcal{F}(0,\infty)$,

• •

Case 2: $f\in\mathcal{F}(0,L_{f})$ and $g\in\mathcal{F}(0,\infty)$, with $0<L_{f}<\infty$,

• •

Case 3: $f\in\mathcal{F}(m_{f},L_{f})$ and $g\in\mathcal{F}(0,\infty)$, with $0<m_{f}\leq L_{f}<\infty$.

We will see that for each case only one term in the Lyapunov function needs to be modified to obtain the convergence result. We then use the matrix inequality condition for each case to obtain information about optimal choices of the algorithm parameters both symbolically and numerically.

III-B Characterization of Fixed Points

From relation (4), the iterates (2) can be rewritten as

[TABLE]

The fixed points of (10) satisfy

[TABLE]

Since $y_{\star}=z_{\star}$, the rightmost equality is exactly the optimality condition for (1).

We will also make use of the following relation, obtained from adding the equations in (10) and the definition of $\phi$,

[TABLE]

From this, we can interpret the feedback nonlinearity $\phi$ as the optimality residual of problem (1), which is driven to zero by the linear system in the feedback interconnection.

III-C Convergence Certificates via Matrix Inequalities

III-C1 Case 1: Non-strongly convex and non-smooth case

We first assume that $f,g\in\mathcal{F}(0,\infty)$. We propose the following family of Lyapunov functions parameterized by a sequence $\{\theta_{i}\}_{i=0}^{\infty}$ with $\theta_{i}>0$,

[TABLE]

for all $k>0$ and $V_{0}=\|x_{0}-x_{\star}\|^{2}$. For notational convenience define the partial sums $\Theta_{k}=\sum_{i=0}^{k-1}\theta_{i}$. The presence of the running sum of subgradients is reminiscent of the Popov criterion [26]. It can also be interpreted as the running weighted sum of fixed point residuals (see (11)). The next lemma shows how this Lyapunov function can ensure a convergence rate in terms of the growth of $\Theta_{k}$.

Lemma 1

Consider the algorithm in (2). Suppose there exists a sequence $\{\theta_{i}\}_{i=0}^{\infty}$ with $\theta_{i}>0$ such that $V_{k+1}\leq V_{k}$ for all $k\geq 0$. Then

[TABLE]

Proof:

Since $V_{k+1}\leq V_{k}$ for all $k$, in particular we have that $V_{k}\leq V_{0}$, or

[TABLE]

Removing the first term on the left, and dividing through by $\Theta_{k}$ gives

[TABLE]

The result follows from the fact that the left side is a weighted average, as $\theta_{i}>0$ and $\sum_{i=0}^{k-1}\theta_{i}/\Theta_{k}=1$. ∎

In the following theorem, we derive a matrix inequality in terms of $\alpha,\lambda$, and $\{\theta_{i}\}_{i=0}^{\infty}$ as a sufficient condition to guarantee $V_{k+1}\leq V_{k}$, which in turn implies (14).

Theorem 1

Let $m_{f}=0$, $L_{f}=\infty$, and consider the following matrix inequality

[TABLE]

where

[TABLE]

If $\alpha,\lambda_{k},\theta_{k}>0$ and $\sigma_{k}^{(1)},\sigma_{k}^{(2)}\geq 0$ are chosen so that (17) is satisfied for all $k\geq 0$, then for all $f\in\mathcal{F}(0,\infty)$ and $g\in\mathcal{F}(0,\infty)$ the iterates in (2) satisfy

[TABLE]

Proof:

We first see that $V_{k+1}-V_{k}$ can be written as a quadratic form. Define the error signal

[TABLE]

Using the updates in (10), the fact that $z_{\star}=y_{\star}$ (see (11)), and the relation (12), it can be verified that

[TABLE]

where $W_{k}^{(0)}$ is given by (18a). Next, note that

[TABLE]

where $Q_{p}(m_{f},L_{f},\alpha)$ is defined in (7). Since $y_{k}=\mathrm{prox}_{\alpha f}(x_{k})$ and $y_{\star}=\mathrm{prox}_{\alpha f}(x_{\star})$, this is exactly the incremental QC that the $\mathrm{prox}_{\alpha f}$ operator satisfies. Thus, we have for all $k$, $e_{k}^{\top}Q^{(1)}e_{k}\geq 0$. We also note that

[TABLE]

As $z_{k}=\mathrm{prox}_{\alpha g}(2y_{k}-x_{k})$ and $z_{\star}=\mathrm{prox}_{\alpha g}(2y_{\star}-x_{\star})$, we similarly conclude that $e_{k}^{\top}Q^{(2)}e_{k}\geq 0$ is implied from the incremental QC that $\mathrm{prox}_{\alpha g}$ satisfies. Returning to (17), if we multiply from the left and right by $e_{k}^{\top}$ and $e_{k}$ respectively, we obtain

[TABLE]

Since $\sigma_{k}^{(1)},\sigma_{k}^{(2)}\geq 0$ and we have shown that $e_{k}^{\top}Q^{(1)}e_{k}\geq 0$ and $e_{k}^{\top}Q^{(2)}e_{k}\geq 0$, it must be that $e_{k}^{\top}W_{k}^{(0)}e_{k}\leq 0$. Hence, $V_{k+1}\!-\!V_{k}\leq 0$, and the result now follows from Lemma 1. ∎

III-C2 Case 2: Non-strongly convex and smooth $f$

If $f\in\mathcal{F}(0,L_{f})$ with $0<L_{f}<\infty$, we may leverage the smoothness of $f$ to refine the result of the previous section. In particular, we can use the inequality for $L_{f}$-smooth functions (see Preliminaries) to relate the behavior of the subgradients to the objective values, whereas in the previous section this inequality was not available. For $k>0$ let

[TABLE]

with $V_{0}$ defined as in the previous case.

This Lyapunov function leads to the following Lemma, the proof of which is identical to that of Lemma 1.

Lemma 2

Consider the algorithm in (2). Suppose there exists a sequence $\{\theta_{i}\}_{i=0}^{\infty}$ with $\theta_{i}>0$ such that $V_{k+1}\leq V_{k}$ for all $k\geq 0$. Then

[TABLE]

This allows us to prove the following theorem for when $f$ is non-strongly convex and smooth.

Theorem 2

Let $0=m_{f}<L_{f}<\infty$ and consider the following matrix inequality

[TABLE]

where

[TABLE]

and $Q^{(1)}$ and $Q^{(2)}$ are defined in (18). If $\alpha,\lambda_{k},\theta_{k}>0$ and $\sigma_{k}^{(1)},\sigma_{k}^{(2)}\geq 0$ are chosen so that (25) is satisfied for all $k\geq 0$, then for all $f\in\mathcal{F}(0,L_{f})$ with $0<L_{f}<\infty$ and $g\in\mathcal{F}(0,\infty)$ the iterates in (2) satisfy

[TABLE]

Proof:

We begin by bounding the difference of the Lyapunov function defined in (24), $V_{k+1}-V_{k}$, by a quadratic form in the error signal $e_{k}$ (see (20)). From the convexity and smoothness of $f$, we can write

[TABLE]

where we have used that $z_{\star}=y_{\star}$. From the convexity of $g$

[TABLE]

Adding these three inequalities together and using the relation (12) allows us to conclude

[TABLE]

Using the recursion for $x_{k+1}$, we then find that

[TABLE]

The proof now proceeds identically as in the proof of Theorem 1 up to the statement that (25) implies $e_{k}^{\top}W_{k}^{(1)}e_{k}\leq 0$. Then by (32), we have that $V_{k+1}-V_{k}\leq 0$, and the result follows from Lemma 2. ∎

III-C3 Case 3: Strongly convex and smooth $f$

We now assume that $f\in\mathcal{F}(m_{f},L_{f})$ and $g\in\mathcal{F}(0,\infty)$, with $0<m_{f}\leq L_{f}<\infty$. For this scenario let

[TABLE]

The following lemma characterizes when we can extract a linear convergence rate from this Lyapunov function.

Lemma 3

Consider the algorithm in (2). Suppose there exists $\rho\in(0,1)$ such the Lyapunov function $V_{k}$ defined by (33) satisfies $V_{k+1}\leq\rho^{2}V_{k}$ for all $k\geq 0$. Then

[TABLE]

Proof:

The proof follows immediately $V_{k+1}\leq\rho^{2}V_{k}$, the definition (33) of $V_{k}$, and induction. ∎

We again see that the difference $V_{k+1}-\rho^{2}V_{k}$ can be written as a quadratic form acting on the error signal $e_{k}$ as defined in (20). Using the definition for $x_{k+1}$ in terms of the previous iterates, we can write

[TABLE]

where $Q_{k}$ is given by

[TABLE]

From (35), we can use the same reasoning developed in the non-strongly convex case to arrive at the following theorem.

Theorem 3

Let $0<m_{f}\leq L_{f}<\infty$ and consider the following matrix inequality

[TABLE]

where $Q_{k}$ is given in (36) and $Q^{(1)}$ and $Q^{(2)}$ are given in (18). If $\alpha,\lambda_{k}>0$, $\sigma_{k}^{(1)},\sigma_{k}^{(2)}\geq 0$, and $\rho\in(0,1)$ are chosen so that (37) is satisfied for all $k\geq 0$, then for all $f\in\mathcal{F}(m_{f},L_{f})$ and $g\in\mathcal{F}(0,\infty)$ with $0<m_{f}\leq L_{f}<\infty$, the iterates in (2) satisfy the following linear convergence rate

[TABLE]

Proof:

We proceed identically as in the proof of Theorems 1 and 2. If (37) is satisfied, then $e_{k}^{\top}Q_{k}e_{k}\leq 0$ for all $k$. This is equivalent to $V_{k+1}-\rho^{2}V_{k}\leq 0$, by which (38) follows from Lemma 3. ∎

IV Optimizing the Bound and Relaxation Parameter

For each of the cases presented in the previous section we use the associated matrix inequality to optimize our bounds on the convergence rate. In the case of non-strongly convex $f$ we provide analytic convergence rates, while for strongly convex $f$ we optimize our bounds numerically.

IV-1 Case 1: Non-strongly convex and non-smooth case

We now select algorithm parameters that satisfy the matrix inequality in (17). In doing so, we arrive at a new and simple proof of the $O(1/k)$ convergence of DRS in the non-strongly convex and non-smooth case.

Theorem 4

If $f\in\mathcal{F}(0,\infty)$ and $g\in\mathcal{F}(0,\infty)$, then for any choice of $\lambda_{k}=\lambda\in(0,2)$ and $\alpha>0$, if we set $\sigma_{k}^{(1)}=\sigma_{k}^{(2)}=\sigma_{k}$, with

[TABLE]

then $\sigma_{k}^{(1)},\sigma_{k}^{(2)},\alpha,\lambda_{k},$ and $\theta_{k}$ satisfy the matrix inequality (17).

Proof:

Making these substitutions gives $W_{k}^{(0)}+\sigma_{k}^{(1)}Q^{(1)}+\sigma_{k}^{(2)}Q^{(2)}=0_{n}\preceq 0.$ ∎

Remark 1

The convergence rate bound provided by the parameter choices in Theorem 4 guarantees convergence only if $\lim_{k\to\infty}\Theta_{k}=\infty$, which in this case means $\sum_{i=0}^{\infty}\lambda_{k}(2-\lambda_{k})=\infty$. This is consistent with the conditions on the relaxation parameter found in [11].

Remark 2

After setting $\theta_{k}$ as in (39), we can maximize $\Theta_{k}$ by setting $\lambda_{k}=1$ which results in $\Theta_{k}=\alpha^{2}k$ and the following result

[TABLE]

Remark 3

Using the relation (12) we can rewrite (40) as

[TABLE]

Thus we see that Theorem 4 also gives a $O(1/k)$ rate toward the iterates being a fixed point of the algorithm.

IV-2 Case 2: Non-strongly convex and smooth case

When $f\in\mathcal{F}(0,L_{f}),g\in\mathcal{F}(0,\infty)$ with $0<L_{f}<\infty$, we have the following result on feasibility of the matrix inequality (25).

Theorem 5

For any $\alpha>0$ and $\lambda_{k}=\lambda$ with $0<\lambda<2$, if we set $\sigma_{k}^{(1)}=\sigma_{k}^{(2)}=\sigma_{k}$,

[TABLE]

Then $\sigma_{k}^{(1)},\sigma_{k}^{(2)},\alpha,\lambda_{k}$, and $\theta_{k}$ satisfy the matrix inequality (25).

Proof:

This can be verified by substituting the expressions for $\sigma_{k}^{(1)}$, $\sigma_{k}^{(2)}$, and $\theta_{k}$ into the minors of $W_{k}^{(1)}+\sigma_{k}^{(1)}Q^{(2)}+\sigma_{k}^{(2)}Q^{(2)}$ and seeing that Sylverster’s criterion is satisfied [27]. ∎

Remark 4

Note that for a constant selection of $\theta_{k}=\theta$, $\Theta_{k}=\theta k$. If (25) holds then

[TABLE]

Remark 5

For moderate values of $\alpha L_{f}$, we can take a second order Taylor expansion of the rightmost term in (42b) and maximize the resulting expression with respect to $\lambda_{k}$. This suggests that we should set $\lambda_{k}$ to

[TABLE]

IV-3 Case 3: Strongly convex and smooth case

When $f\in\mathcal{F}(m_{f},L_{f})$, $g\in\mathcal{F}(0,\infty)$, with $0<m_{f}\leq L_{f}<\infty$, we can modify the matrix inequality (37) to get a linear dependence on the relaxation parameter. If we define $\Lambda_{k}:=\begin{bmatrix}0&-\lambda_{k}&\lambda_{k}\end{bmatrix}^{\top}\otimes I_{d}$ and

[TABLE]

then (37) is equivalent to

[TABLE]

As $\Lambda_{k}\Lambda_{k}^{\top}\succeq 0$, if (46) is satisfied then it must be the case that $M_{k}\preceq 0$. We now recognize that $M_{k}-\Lambda_{k}[-1]\Lambda_{k}^{\top}$ is the Schur Complement of the bottom right entry in the matrix

[TABLE]

By the properties of the Schur complement [28], we can conclude that (37) is satisfied if and only if $\Sigma_{k}\preceq 0$. The advantage of using $\Sigma_{k}\preceq 0$ instead of (37), is that now both the convergence rate $\rho^{2}$ and the relaxation parameter $\lambda_{k}$ appear linearly. If we set $\sigma_{k}^{(1)}=\sigma^{(1)}$, $\sigma_{k}^{(2)}=\sigma^{(2)}$, and $\lambda_{k}=\lambda$ for all $k$, the corresponding bound in (38) can be optimized by solving the following SDP, where $\lambda$ is now a decision variable,

[TABLE]

where the decision variables are $\rho^{2},\lambda>0$ and $\sigma^{(1)},\sigma^{(2)}\geq 0$. The optimal $\rho^{2}$ from solving this program over a range of step sizes $\alpha$ and condition numbers $\kappa_{f}=L_{f}/m_{f}$ is shown in Figure 1. We see that with increasing $\kappa_{f}$, the optimal choice of $\alpha$ decreases. Across the range of values of $\kappa_{f}$ and $\alpha$, the SDP (48) returns $\lambda=2$ as the optimal relaxation parameter.

Remark 6

We may repeat the same derivation for the inequalities (17) and (25) to linearize their dependence on $\lambda_{k}$ as well.

V Numerical Experiments

We investigate how our theoretical results compare with the experimental performance of DRS in the three scenarios described above. For the non-smooth and non-strongly convex case, we consider a basis pursuit problem (see [3]),

[TABLE]

with data $A\in\mathbb{R}^{300\times 10000}$ and $b\in\mathbb{R}^{300}$. We run DRS on the dual of this problem (ADMM) which is non-smooth and non-strongly convex. We test the convergence over a range of values of $\lambda_{k}=\lambda$, including $\lambda^{\star}=1$ (see Remark 2), and fixed $\alpha=1$.

For the smooth cases we consider a LASSO problem,

[TABLE]

with $A\in\mathbb{R}^{300\times 200}$, $b\in\mathbb{R}^{300}$ and $\gamma=0.1$. For the non-strongly convex case we set $A$ to be rank deficient and plot the convergence of DRS over a range of $\alpha$ with $\lambda_{k}=1$ fixed. The value $\alpha^{\star}$ is found by performing a grid search over possible values and choosing that which maximizes the rate given by Theorem 2. For the strongly convex case $A$ is set to be full rank and plot the convergence of DRS over a range of $\alpha$ with $\lambda_{k}=2$ fixed. Again, $\alpha^{\star}$ is the value of $\alpha$ which gives the best rate as provided by Theorem 3. The results are presented in Figure 2.

VI Conclusion

We presented a unified framework for deriving convergence bounds for DRS and parameter settings that optimize these bounds. Our framework encompasses different assumptions on the smoothness and convexity of $f$. We are able to give simple proofs of convergence and find optimal choices for the relaxation parameter by solving a small convex program for a fixed $\alpha$. It is important to note that the parameter selections optimize our bounds in the sense of the best worst-case convergence rate over the entire class of objective functions with $f\in\mathcal{F}(m_{f},L_{f})$. While there are scenarios where additional structure in the problem might make alternate parameter settings more effective, the settings we see here bound the worst-case performance agnostic of any additional problem structure. For future work, this framework will be extended to encompass accelerated variants of DRS, as well as three or more operator splitting and multi-block ADMM.