A Control-Theoretic Approach to Analysis and Parameter Selection of Douglas-Rachford Splitting
Jacob H. Seidman Mahyar Fazlyab Victor M. Preciado George J. Pappas




[
    Loren N. Pfeiffer
    December 30, 2023
=====================

empty
empty





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.


Optimization algorithms, Lyapunov methods



§ INTRODUCTION


























In this paper, we consider problems of the form

    minimize_x ∈ℝ^d {F(x)=f(x)+g(x)},

where f,g:^d →∪{+∞} are convex, closed, and proper (c.c.p.).  Douglas-Rachford splitting (DRS) solves problem (<ref>) with the following iterations:
 

    y_k    = _α f(x_k),  
    
    z_k    =_α g(2y_k-x_k),  
    
    x_k+1   = x_k + λ_k (z_k-y_k),


where  is the proximal operator (see Definition <ref>) and α and λ_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 (<ref>).  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 α and λ_k using semidefinite programs (SDPs).

The algorithm was first proposed in <cit.>, and has since found application in general separable optimization problems <cit.>.  Its dual formulation, ADMM, has been particularly useful in distributed optimization problems <cit.>.  Since the iterates of ADMM can be written as applying DRS to the dual problem <cit.>, 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 <cit.>, <cit.>.  These techniques have led to proofs of a O(1/k) convergence rate for the non-strongly convex case <cit.>, and linear convergence when f is smooth and strongly convex <cit.>, <cit.>.  In the most general case, a condition for convergence is that λ_k ∈ (0,2) with ∑_k=0^∞λ_k(2-λ_k) = ∞ <cit.>, though there exist cases where the algorithm converges with λ_k > 2 <cit.>.

Recently, there has been interest in automating the analysis and design of optimization algorithms via SDPs, <cit.>. In particular, 



through the method of integral quadratic constraints proposed in <cit.>, the authors of <cit.> 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 <cit.> provide evidence that as the relaxation parameter approaches 2 from below, the linear convergence rate is close to being optimal and in <cit.> are able to analytically solve the SDP to give a convergence rate.  The work in <cit.> gives an optimal choice for the relaxation parameter when f is quadratic.  Furthermore, <cit.> gives a set of assumptions in which a bound on the linear convergence rate is minimized by setting λ_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 <cit.> 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 λ_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 λ_k.




§ PRELIMINARIES


We denote the set of real numbers by ℝ, the set of real n-dimensional vectors by ℝ^n, the set of real m× n-dimensional matrices by ℝ^m× n, and the n-dimensional identity matrix and zero matrix by I_n and 0_n, respectively.  For a function f ℝ^d →ℝ∪{+∞}, we denote by f={x ∈ℝ^d  f(x)<∞} the effective domain of f. The subdifferential of a function f at a point x is ∂ f(x) := {g | f(y) - f(x) ≥ g^⊤(y - x), ∀ y ∈(f)}.  By abuse of notation we will also refer to a subgradient, that is an element of the subdifferential by ∂ f(x) as well.  The indicator function of a set C is given by 1_C(x) = 0 if x ∈ C and 1_C(x) = ∞ if x ∉ C.  For two matrices A ∈ℝ^m× n and B∈ℝ^p× q their Kronecker product is A ⊗ B.  
 
We say a differentiable function f ℝ^d →ℝ is L_f-smooth on 𝒮⊆  f if ∇ f(x)-∇ f(y)_2 ≤ L_f x-y_2 for some L_f>0 and all x,y ∈𝒮. This also implies for all x,y ∈𝒮, f(y) ≤ f(x) + ∇ f(x)^⊤ (y-x) +(L_f/2) y-x_2^2.  A differentiable function f ℝ^d →ℝ is m_f-strongly convex on 𝒮⊆  f if m_f x-y_2^2 ≤ (x-y)^⊤ (∇ f(x)-∇ f(y)) for some m_f>0 and all x,y ∈𝒮.  The class of functions which are L_f-smooth and m_f-strongly convex is denoted by ℱ(m_f,L_f). 

 
	Given a c.c.p. function f:^d →∪{+∞} and α > 0, the proximal operator _α f:^d →^d is defined as
	
    _α f(x) = argmin_y {f(y) + 1/2αx - y_2^2}.

	The point y = _α f(x) also is given by the implicit solution to the subgradient equation
	
    y = x - α∂ f(y).




We say that a nonlinear function ϕ:^d →^d satisfies the incremental quadratic constraint <cit.> (or point-wise integral quadratic constraint <cit.>) defined by Q∈^2d × 2d if for all x,y,

    [       x-y; ϕ(x)-ϕ(y) ]^⊤ Q [       x-y; ϕ(x)-ϕ(y) ]≥ 0.



For 0≤ m, L < ∞ define

	
    Q(m,L) = [ -m L/m+L      1/2;      1/2   -1/m+L ]⊗ I_d,


and define Q(m,∞) as lim_L →∞ Q(m,L).  It was noted in <cit.> that a differentiable function f belongs to the class ℱ(m_f,L_f) on 𝒮 if and only if the gradient ∇ f satisfies the incremental QC defined by Q(m_f,L_f).  When L_f = ∞, the subgradient ∂ f satisfies the QC defined by Q(m_f,∞).  If we define

    Q_p(m,L,α) = [     0   I_d; α I_d  -I_d ]Q(m,L) [     0 α I_d;   I_d  -I_d ],

then the proximal operator of a function f ∈ℱ(m,L), _α f, satisfies the incremental QC defined by Q_p(m,L,α) <cit.>.

























§ ANALYSIS OF DOUGLAS-RACHFORD SPLITTING VIA MATRIX INEQUALITIES




 §.§ Douglas-Rachford Splitting as a Dynamical System

We can write the updates in (<ref>) as a linear system with state x_k and feedback nonlinearity ϕ(x_k),

    x_k+1 = x_k + λ_k ϕ(x_k),

where

    ϕ(x_k) := _α g(2_α f(x_k)-x_k)-_α g(x_k).


Our main technique is to describe the nonlinearity ϕ with incremental QCs representing the  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 ∈ℱ(0,∞) and g ∈ℱ(0,∞),


  * Case 2: f ∈ℱ(0,L_f) and g ∈ℱ(0,∞), with 0 < L_f < ∞,


  * Case 3: f ∈ℱ(m_f,L_f) and g ∈ℱ(0,∞), with 0 < m_f ≤ L_f < ∞.


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.




 §.§ Characterization of Fixed Points


From relation (<ref>), the iterates (<ref>) can be rewritten as

 

    y_k    = x_k -α∂ f(y_k),  
    
    z_k    =2y_k-x_k - α∂ g(z_k), 
    
    x_k+1   = x_k + λ_k (z_k-y_k) .



The fixed points of (<ref>) satisfy


    x_⋆ = y_⋆+α∂ f(y_⋆),  y_⋆=z_⋆,  ∂ f(y_⋆) +∂ g(z_⋆)=0.


Since y_⋆ = z_⋆, the rightmost equality is exactly the optimality condition for (<ref>). 


We will also make use of the following relation, obtained from adding the equations in (<ref>) and the definition of ϕ,


    ϕ(x_k) = z_k - y_k = -α(∂ f(y_k) + ∂ g(z_k)).


From this, we can interpret the feedback nonlinearity ϕ as the optimality residual of problem (<ref>), which is driven to zero by the linear system in the feedback interconnection.




 §.§ Convergence Certificates via Matrix Inequalities




  §.§.§ Case 1: Non-strongly convex and non-smooth case


We first assume that f,g ∈ℱ(0,∞).  We propose the following family of Lyapunov functions parameterized by a sequence {θ_i}_i=0^∞ with θ_i > 0,



    V_k =x_k-x_⋆_2^2+ ∑_i=0^k-1θ_i ∂ f(y_i) + ∂ g(z_i)_2^2,


for all k > 0 and V_0 = x_0 - x_⋆^2.  For notational convenience define the partial sums Θ_k = ∑_i=0^k-1θ_i.  The presence of the running sum of subgradients is reminiscent of the Popov criterion <cit.>.  It can also be interpreted as the running weighted sum of fixed point residuals (see (<ref>)).  The next lemma shows how this Lyapunov function can ensure a convergence rate in terms of the growth of Θ_k.
 
	Consider the algorithm in (<ref>). Suppose there exists a sequence {θ_i}_i=0^∞ with θ_i>0 such that V_k+1≤ V_k for all k≥ 0. Then 
	
    min_i=0,…,k-1∂ f(y_i) + ∂ g(z_i)_2^2 ≤1/Θ_kx_0-x_⋆_2^2.

	



Since V_k+1≤ V_k for all k, in particular we have that V_k ≤ V_0, or


    x_k-x_⋆_2^2+∑_i=0^k-1θ_i∂ f(y_i) +∂ g(z_i)_2^2 ≤x_0-x_⋆_2^2.


Removing the first term on the left, and dividing through by Θ_k gives

    ∑_i=0^k-1θ_i/Θ_k∂ f(y_i) + ∂ g(z_i)_2^2 ≤x_0-x_ ⋆_2^2Θ_k.

	The result follows from the fact that the left side is a weighted average, as θ_i > 0 and ∑_i=0^k-1θ_i/Θ_k = 1.


In the following theorem, we derive a matrix inequality in terms of α,λ, and {θ_i}_i=0^∞ as a sufficient condition to guarantee V_k+1≤ V_k, which in turn implies (<ref>).


























	Let m_f = 0, L_f = ∞, and consider the following matrix inequality
	
    W_k^(0) + σ_k^(1) Q^(1) + σ_k^(2) Q^(2)≼ 0,

	where 



    W_k^(0)   = [                  0               -λ_k                λ_k;               -λ_k     λ_k^2 +θ_k/α^2 -(λ_k^2 + θ_k/α^2);                λ_k  -(λ_k^2 +θ_k/α^2)     λ_k^2 +θ_k/α^2 ]⊗I_d , 
    
    Q^(1)   = [     0   I_d; α I_d - I_d;     0     0 ]Q(m_f,L_f)  [     0 α I_d     0;   I_d - I_d     0 ], 
    
    Q^(2)   =[     0  -I_d;     0 2 I_d; α I_d - I_d ]Q(0,∞) [     0     0 α I_d; - I_d 2 I_d - I_d ].



	If  α, λ_k,θ_k > 0 and σ_k^(1), σ_k^(2)≥ 0 are chosen so that (<ref>) is satisfied for all k ≥ 0, then for all f ∈ℱ(0,∞) and g ∈ℱ(0,∞) the iterates in (<ref>) satisfy 


    min_i=0,…,k-1∂ f(y_i) + ∂ g(z_i)_2^2 ≤1/Θ_kx_0-x_⋆_2^2.



We first see that V_k+1 - V_k can be written as a quadratic form.  Define the error signal


    e_k :=[ (x_k-x_⋆)^⊤ (y_k-y_⋆)^⊤ (z_k-z_⋆)^⊤ ]^⊤.


Using the updates in (<ref>), the fact that z_⋆ = y_⋆ (see (<ref>)), and the relation (<ref>), it can be verified that


    V_k+1 - V_k = e_k^⊤ W_k^(0) e_k,

where W_k^(0) is given by (<ref>). Next, note that

    e_k^⊤ Q^(1) e_k = x_k - x_⋆
     y_k - y_⋆^⊤ Q_p(m_f,L_f,α) x_k - x_⋆
     y_k - y_⋆,
  
where Q_p(m_f,L_f,α) is defined in (<ref>).  Since y_k = _α f(x_k) and y_⋆ =_α f(x_⋆), this is exactly the incremental QC that the _α f operator satisfies.  Thus, we have for all k, e_k^⊤ Q^(1) e_k ≥ 0. We also note that

    [     0     0 α I_d; - I_d 2 I_d - I_d ]e_k  = 0    αI_d 
     I_d    -I_d(2y_k - x_k) - (2y_⋆ -x_⋆) 
     z_k - z_⋆.


As z_k = _α g(2y_k - x_k) and z_⋆ = _α g(2y_⋆ - x_⋆), we similarly conclude that e_k^⊤ Q^(2) e_k ≥ 0 is implied from the incremental QC that _α g satisfies. Returning to (<ref>), if we multiply from the left and right by e_k^⊤ and e_k respectively, we obtain

    e_k^⊤ W_k^(0) e_k + σ_k^(1) e_k^⊤ Q^(1) e_k + σ_k^(2) e_k^⊤ Q^(2) e_k ≤ 0.

Since σ_k^(1), σ_k^(2)≥ 0 and we have shown that e_k^⊤ Q^(1) e_k ≥ 0 and e_k^⊤ Q^(2) e_k ≥ 0, it must be that e_k^⊤ W_k^(0) e_k ≤ 0.  Hence, V_k+1- V_k ≤ 0, and the result now follows from Lemma <ref>.





  §.§.§ Case 2: Non-strongly convex and smooth f

If f ∈ℱ(0,L_f) with 0< L_f< ∞, 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 

    V_k = x_k-x_⋆_2^2 +  ∑_i=0^k-1θ_i[F(z_i)-F(z_⋆)],

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 <ref>.

 
	Consider the algorithm in (<ref>). Suppose there exists a sequence {θ_i}_i=0^∞ with θ_i>0 such that V_k+1≤ V_k for all k ≥0. Then 
	
    min_i=0,…,k-1 [F(z_i) - F(z_⋆)] ≤1Θ_kx_0 - x_⋆_2^2.



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

 
Let 0 = m_f < L_f < ∞ and consider the following matrix inequality
	
    W_k^(1) + σ_k^(1) Q^(1) + σ_k^(2) Q^(2)≼ 0,

	where 
	
    W_k^(1) =[                       0                    -λ_k                     λ_k;                    -λ_k       θ_k L_f/2 + λ_k^2 θ_k/2(1/α - L_f) -λ_k^2;                     λ_k  θ_k/2(1/α- L_f) -λ_k^2 θ_k(L_f/2 - 1/α) +λ_k^2 ]⊗ I_d,

	and Q^(1) and Q^(2) are defined in (<ref>).  If α, λ_k,θ_k > 0 and σ_k^(1), σ_k^(2)≥ 0 are chosen so that (<ref>) is satisfied for all k ≥ 0, then for all f ∈ℱ(0,L_f) with 0<L_f<∞ and g ∈ℱ(0,∞) the iterates in (<ref>) satisfy


    min_i=0,…,k-1 [F(z_i) - F(z_⋆)] ≤1/Θ_kx_0-x_⋆_2^2.




We begin by bounding the difference of the Lyapunov function defined in (<ref>), V_k+1 - V_k, by a quadratic form in the error signal e_k (see (<ref>)).   From the convexity and smoothness of f, we can write

    f(z_k) - f(y_k)   ≤∇ f(y_k)^⊤(z_k - y_k) +L_f/2z_k - y_k_2^2, 
    
    f(y_k) - f(z_⋆)   ≤∇ f(y_k)^⊤(y_k -y_⋆),

where we have used that z_⋆ = y_⋆.  From the convexity of g

    g(z_k) -g(z_⋆) ≤∂ g(z_k)^⊤(z_k - z_⋆).

Adding these three inequalities together and using the relation (<ref>) allows us to conclude

    F(z_k) - F(z_⋆)    ≤L_f/2y_k-y_⋆_2^2 + (L_f/2 - 1/α)z_k - z_⋆_2^2
         + (1/α - L_f)(y_k-y_⋆)^⊤(z_k - z_⋆).

Using the recursion for x_k+1, we then find that 

    V_k+1 - V_k ≤ e_k^⊤ W_k^(1) e_k.


The proof now proceeds identically as in the proof of Theorem <ref> up to the statement that (<ref>) implies e_k^⊤ W_k^(1) e_k ≤ 0.  Then by (<ref>), we have that V_k+1 - V_k ≤ 0, and the result follows from Lemma <ref>.





  §.§.§ Case 3: Strongly convex and smooth f

We now assume that f ∈ℱ(m_f,L_f) and g ∈ℱ(0, ∞), with 0 < m_f ≤ L_f < ∞.  For this scenario let



    V_k = x_k - x_∗_2^2.


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

 
Consider the algorithm in (<ref>). Suppose there exists ρ∈ (0,1) such the Lyapunov function V_k defined by (<ref>) satisfies V_k+1≤ρ^2 V_k for all k ≥ 0.  Then 

    x_k-x_∗_2^2≤ρ^2kx_0 - x_∗_2^2.





The proof follows immediately V_k+1≤ρ^2 V_k, the definition (<ref>) of V_k, and induction.



We again see that the difference V_k+1 - ρ^2 V_k can be written as a quadratic form acting on the error signal e_k as defined in (<ref>).  Using the definition for x_k+1 in terms of the previous iterates, we can write


    V_k+1 - ρ^2 V_k = e_k^⊤ Q_ke_k,

where Q_k is given by


    Q_k = [ 1- ρ^2   -λ_k    λ_k;   -λ_k  λ_k^2 -λ_k^2;    λ_k -λ_k^2  λ_k^2 ]⊗ I_d.


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



	Let 0 < m_f ≤ L_f < ∞ and consider the following matrix inequality
	
    Q_k + σ_k^(1) Q^(1) + σ_k^(2) Q^(2)≼ 0,

	where Q_k is given in (<ref>) and Q^(1) and Q^(2) are given in (<ref>).  If α, λ_k > 0, σ_k^(1), σ_k^(2)≥ 0, and ρ∈ (0,1) are chosen so that (<ref>) is satisfied for all k ≥ 0, then for all f ∈ℱ(m_f,L_f) and g ∈ℱ(0,∞) with 0<m_f≤ L_f<∞, the iterates in (<ref>) satisfy the following linear convergence rate 


    x_k - x_∗_2^2 ≤ρ^2kx_0 - x_∗_2^2.




We proceed identically as in the proof of Theorems <ref> and <ref>.  If (<ref>) is satisfied, then e_k^⊤ Q_k e_k ≤ 0 for all k.  This is equivalent to V_k+1 -  ρ^2 V_k ≤ 0, by which (<ref>) follows from Lemma <ref>.




§ 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.



  §.§.§ Case 1: Non-strongly convex and non-smooth case


We now select algorithm parameters that satisfy the matrix inequality in (<ref>).  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.

 
If f ∈ℱ(0,∞) and g∈ℱ(0,∞), then for any choice of λ_k = λ∈ (0,2) and α>0, if we set σ_k^(1) = σ_k^(2) = σ_k, with


    σ_k:=2λ_k/α,   θ_k := α^2λ_k(2-λ_k),


then σ_k^(1), σ_k^(2), α, λ_k, and θ_k satisfy the matrix inequality (<ref>).


Making these substitutions gives W_k^(0) + σ_k^(1)Q^(1) + σ_k^(2)Q^(2) = 0_n ≼ 0.







The convergence rate bound provided by the parameter choices in Theorem <ref> guarantees convergence only if lim_k →∞Θ_k = ∞, which in this case means ∑_i=0^∞λ_k(2-λ_k) = ∞.  This is consistent with the conditions on the relaxation parameter found in <cit.>.


 
After setting θ_k as in (<ref>), we can maximize Θ_k by setting λ_k = 1 which results in Θ_k = α^2 k and the following result

    min_i=0,…,k-1∂ f(y_i) + ∂ g(z_i)_2^2 ≤1/α^2 kx_0-x_⋆_2^2.




Using the relation (<ref>) we can rewrite (<ref>) as

    min_i=0,…,k-1z_i - y_i_2^2 ≤1/kx_0-x_⋆_2^2.

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




  §.§.§ Case 2: Non-strongly convex and smooth case


When f ∈ℱ(0,L_f), g∈ℱ(0,∞) with 0<L_f<∞, we have the following result on feasibility of the matrix inequality (<ref>).

 
For any α > 0 and λ_k = λ with 0 < λ < 2, if we set σ_k^(1)=σ_k^(2) =σ_k,


    σ_k    :=2λ_k/α(√((2-λ_k/α L_f)^2 +1) - (2-λ_k/α L_f)), 
    θ_k    := 2λ_k α(1 + (2-λ_k/α L_f) - √((2-λ_k/α L_f)^2 +1)).


Then σ_k^(1), σ_k^(2), α, λ_k, and θ_k satisfy the matrix inequality (<ref>).


This can be verified by substituting the expressions for σ_k^(1), σ_k^(2), and θ_k into the minors of W_k^(1) + σ_k^(1)Q^(2) + σ_k^(2)Q^(2) and seeing that Sylverster's criterion is satisfied <cit.>.



Note that for a constant selection of θ_k = θ, Θ_k = θ k.  If (<ref>) holds then

    min_i=0,…,k-1 [F(z_i) - F(z_⋆)] ≤1θ kx_0 - x_⋆_2^2.



For moderate values of α L_f, we can take a second order Taylor expansion of the rightmost term in (<ref>) and maximize the resulting expression with respect to λ_k.  This suggests that we should set λ_k to

    λ_k = 2/3(2 - α L_f + √(1 - α L_f + α^2 L_f^2)).






  §.§.§ Case 3: Strongly convex and smooth case

When f ∈ℱ(m_f, L_f), g∈ℱ(0,∞), with 0 < m_f ≤ L_f < ∞, we can modify the matrix inequality (<ref>) to get a linear dependence on the relaxation parameter.  If we define Λ_k := [    0 -λ_k  λ_k ]^⊤⊗ I_d and 

    M_k := (Q_k - Λ_k Λ_k^⊤) +σ_k^(1) Q^(1) + σ_k^(2) Q^(2),

then (<ref>) is equivalent to


    M_k - Λ_k [-1] Λ_k^⊤≼ 0.


As Λ_k Λ_k^⊤≽ 0, if (<ref>) is satisfied then it must be the case that M_k ≼ 0.  We now recognize that M_k- Λ_k [-1] Λ_k^⊤ is the Schur Complement of the bottom right entry in the matrix


    Σ_k := [   M_k   Λ_k; Λ_k^⊤    -1 ].


By the properties of the Schur complement <cit.>, we can conclude that (<ref>) is satisfied if and only if Σ_k ≼ 0.  The advantage of using Σ_k ≼ 0 instead of (<ref>), is that now both the convergence rate ρ^2 and the relaxation parameter λ_k appear linearly.  If we set σ_k^(1) = σ^(1), σ_k^(2) = σ^(2), and λ_k = λ for all k, the corresponding bound in (<ref>) can be optimized by solving the following SDP, where λ is now a decision variable,


    minimize  ρ^2,  subject to  Σ_k ≼ 0,
 

where the decision variables are ρ^2, λ > 0 and σ^(1),σ^(2)≥ 0.  The optimal ρ^2 from solving this program over a range of step sizes α and condition numbers κ_f = L_f/m_f is shown in Figure <ref>.  We see that with increasing κ_f, the optimal choice of α decreases.  Across the range of values of κ_f and α, the SDP (<ref>) returns λ=2 as the optimal relaxation parameter.


We may repeat the same derivation for the inequalities (<ref>) and (<ref>) to linearize their dependence on λ_k as well.


















§ 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 <cit.>),

   x,z∈^nminimize         1_{y∈^n |Ay=b}(x) + z_1 

   subject to         x - z = 0, 
with data A ∈^300 × 10000 and b ∈ℝ^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 λ_k = λ, including λ^⋆ = 1 (see Remark <ref>), and fixed α = 1.

For the smooth cases we consider a LASSO problem,

    minimize_x  1/2Ax-b_2^2 + γx_1.

with A ∈^300 × 200, b ∈^300 and γ = 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 α with λ_k = 1 fixed.  The value α^⋆ is found by performing a grid search over possible values and choosing that which maximizes the rate given by Theorem <ref>.  For the strongly convex case A is set to be full rank and plot the convergence of DRS over a range of α with λ_k = 2 fixed.  Again, α^⋆ is the value of α which gives the best rate as provided by Theorem <ref>.  The results are presented in Figure <ref>.
























§ 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 α.  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 ∈ℱ(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.









































IEEEtran













