A Proximal Alternating Direction Method of Multiplier for Linearly Constrained Nonconvex Minimization
Jiawei Zhang Zhi-Quan Luo

Organic molecular tuning of many-body interaction energies in air-suspended carbon nanotubes
    Y. K. Kato
    December 30, 2023
============================================================================================


  Consider the minimization of a nonconvex differentiable function over a polyhedron. A popular primal-dual first-order method for this problem is to perform a gradient projection iteration for the augmented Lagrangian function and then update the dual multiplier vector using the constraint residual. However, numerical examples show that this approach can exhibit “oscillation" and may not converge. In this paper, we propose a proximal alternating direction method of multipliers for the multi-block version of this problem. A distinctive feature of this method is the introduction of a “smoothed" (i.e., exponentially weighted) sequence of primal iterates, and the inclusion, at each iteration, to the augmented Lagrangian function a quadratic proximal term centered at the current smoothed primal iterate. The resulting proximal augmented Lagrangian function is inexactly minimized (via a gradient projection step) at each iteration while the dual multiplier vector is updated using the residual of the linear constraints. When the primal and dual stepsizes are chosen sufficiently small, we show that suitable “smoothing" can stabilize the “oscillation", and the iterates of the new proximal ADMM algorithm converge to a stationary point under some mild regularity conditions.

The iteration complexity of our algorithm for finding an ϵ-stationary solution is 𝒪(1/ϵ^2), which improves the best known complexity of 𝒪(1/ϵ^3) for the problem under consideration. Furthermore, when the objective function is quadratic, we establish the linear convergence of the algorithm. Our proof is based on a new potential function and a novel use of error bounds.







§ INTRODUCTION

Consider the following linearly constrained optimization problem:

    [                    f(x);       [5pt] Ax=b, x∈ P, ]

where f is differentiable (not necessarily convex)and

    P:={x∈ℝ^n|ℓ_i≤ x_i≤ u_i, i=1,2,...,n }

is a bounded box with ℓ_i<u_i for all i, and the matrix A has dimension ℝ^m× n.
Problems of the form (<ref>) arise in many applications involving big data, including nonnegative matrix factorization <cit.>, phase retrieval <cit.>, distributed matrix factorization <cit.>, polynomial optimization <cit.>, asset allocation <cit.>, zero variance discriminant analysis <cit.>, to name just a few. A popular approach to solve problem (<ref>) is to dualize the linear equality constraint and apply a primal-dual type algorithm to the resulting augmented Lagrangian function. Such approach is particularly attractive when the objective function f(x) has a separable structure since in this case the corresponding primal minimization problem can be decomposed and often times solvable efficiently in parallel, while the dual update can be carried out in closed form.
The local convergence analysis of the classical augmented Lagrangian method was for smooth objective function with smooth equality constraint <cit.>. Global convergence of a primal-dual method for smooth objective and smooth equality constraint was given recently in <cit.>. When the decision variable x consists of many small variable blocks, the popular alternating direction method of multipliers (ADMM) is often the preferred algorithm to solve (<ref>), see <cit.> for a detailed coverage of the method and many applications from a set of diverse fields. It is well known <cit.> that for the two-block (strongly) convex case ADMM can be viewed as a variant of proximal-point method or operator-splitting method, from which one can derive linear convergence of the method. The paper <cit.> proves that the direct extension of ADMM to three-block case may not converges and state some sufficient conditions for the extended algorithm to converge.
The reference <cit.> uses the error bound analysis <cit.> to establish the linear convergence rate of ADMM for a family of convex programming problem with any number of variable blocks with a reduced dual stepsize. However, for nonconvex problems, the convergence of the augmented Lagrangian method or ADMM has not been well understood, despite the fact they have been widely used in applications.
In <cit.>, the convergence of ADMM was established for some special nonconvex problems such as consensus-based sharing problems by using the augmented Lagrangian function as the potential function. This approach was further extended <cit.> to a larger family of nonconvex-nonsmooth problems under some technical assumption such as the prox-regularity of the objective function. The papers <cit.> proved the convergence of some inexact ADMM for certain nonsmoooth, nonconvex problems. However, these references all require at least one block of the variable to be unconstrained and a strong feasibility assumption holds, namely, suppose the linear equality constraint is A_1x_1+A_2x_2=b, then the image of A_1 is contained in the image of A_2 and the variable block x_2 does not have any other constraint.
Recently, <cit.> established the convergence of multi-block ADMM algorithm for the so called multi-affine constraints which are linear in each variable block but otherwise nonconvex. Similar to <cit.>, this paper also has some technical assumptions, including the similar strong feasibility constraint and that the objective function for some block must be strongly convex. Thus, the results from the existing studies <cit.> do not cover the general problem (<ref>).

The penalty method of <cit.> does apply to the linearly constrained nonconvex problem (<ref>). However, the penalty method is usually slower than the method of multipliers which do not use diminishing stepsizes or any unbounded parameters. To our knowledge, the fastest first order algorithm for the linearly constrained problem (<ref>) is the one proposed by <cit.>, which achieves an iteration complexity 𝒪(1/ϵ^3).
The contribution of this paper is as follows. We propose a proximal alternating direction method of multipliers (ADMM) to solve the linearly constrained nonconvex differentiable minimization problem (<ref>). A distinctive feature of the algorithm is to introduce a “smoothed" (i.e., exponentially weighted) sequence of primal iterates, and at each iteration add to the augmented Lagrangian function an extra quadratic proximal term centered at the smoothed primal iterate. The resulting proximal augmented Lagrangian function is inexactly minimized at each iteration while the dual multiplier vector is updated using the residual of the linear constraints. The algorithm is well suited for large scale optimization involving big data, and easily extends to the multiple variable block case, resulting in a variant of the well-known ADMM algorithm. When the primal and dual stepsizes are chosen sufficiently small, we show that the iterates of the proximal ADMM algorithm converge to a stationary point of the nonconvex problem under some mild regularity conditions. Furthermore, we show that this algorithm has an iteration complexity 𝒪(1/ϵ^2) to find an ϵ-stationary solution of problem (<ref>), which improves the best known complexity of 𝒪(1/ϵ^3) <cit.>. Moreover, we present a numerical example showing that the “smoothing" step is necessary for the convergence of the proximal ADMM when the objective function is nonconvex. Furthermore, for a quadratic objective function, we establish the linear convergence of the algorithm.




§ PRELIMINARIES



 §.§ Some notations in this paper

We define some notations used in this paper.


  * ℝ^n is the n-dimension Euclidian space.

  * ℝ^m× n is the space of all m× n real matrix.

  * [·]_+ that  satisfies [x]_+=min_x̅∈ Px-x̅ is the  projection operator to the box P. In particular, the i-th coordinate of [x]_+ is given by

    ([x]_+)_i={[              ℓ_i,          x_i<ℓ_i,;              x_i, x_i ∈ [ℓ_i, u_i],;              u_i,          x_i>u_i. ].


  * (⋯)_+ is the projection operator to the non-negative orthant.

  * ≽ means the coordinate-wise ≥ and ≼ means the coordinate-wise ≤.




 §.§ The set of stationary solutions

We first define the solution of the problem (<ref>) in this subsection.
Due to the linearity of the constraints, there exists a set of Lagrangian multipliers for each stationary point of (<ref>) such that the KKT condition holds. We denote the set of stationary points of (<ref>) by X^*.
Writing down the KKT condition, letting y, μ, ν be the multipliers, we have:

    ∇ f(x^*)+A^Ty^*-μ^*+ν^*   =   0,
    
    Ax^*   =   b,
    ℓ_i≤   x_i^*   ≤    u_i,  i
    μ^*   ≽    0,
    ν^*   ≽    0,
    μ_i^*(ℓ_i-x_i^*)   =   0,   i
    ν_i^*(x_i^*-u_i)   =   0,  i,

where μ_i^* and ν_i^* denote the i-th component of μ^* and ν^* respectively.
Let X^*, Y^* be the sets of all x^* and y^* satisfying the KKT condition. And let W={(x, y)| x, y is a pair of primal-dual solution }.

Note that (<ref>) and (<ref>) are the complementarity conditions. It means that either (ℓ_i-x_i^*) or μ_i^* must be zero for all i and similarly for ν^*.
A stronger condition, which holds generically, is called “strict complementarity condition Γ”.


If for all solutions (x^*, y^*, μ^*, ν^*) of the KKT system, for any i, exactly one of μ_i^* and (ℓ_i-x_i^*) is zero and exact one of ν_i^* and (x_i^*-u_i) is zero, then we say the original problem satisfies the strict complementarity condition.



 §.§ Assumptions




In this subsection, we give our main assumptions, which are valid in many practical problems.


 
 [-10pt]



  (a) The origin is in the relative interior of the set AP-b={Ax-b| x∈ P}.

  (b) The strict complementarity condition holds for (<ref>).

  (c) The objective function f is a differentiable function with Lipschitz continuous gradient

    ∇ f(x)-∇ f(x')≤ Lx-x',  L>0 ∀  x, x'∈ P.



Note that Assumption <ref>(a) is equivalent to the feasibility of (<ref>) for all slightly perturbed b from the range space of A; in particular it does not require the full row rank of A. By using Cauchy-schwarz inequality, assumption <ref>(c) implies the existence of a constant γ (possibly negative) such that

    ⟨∇ f(x)-∇ f(x'),x-x'⟩≥γx-x'^2,   .

Assumption <ref>(b) is reasonable since the strict complementarity is valid generically, as we argue in the proposition below.

Suppose f(x)=g(x)+v^Tx, g is Lipschitz-differentiable and v is a constant vector. If the data vector (v, b) is generated from a continuous probability distribution, then with probability 1, the strict complementarity condition holds for (<ref>). Here a probability distribution is said to be continuous if the probability of a Lesbergue-zero-measure set is zero.


We will use the fact that Lipschitz continuous functions map a zero-measure set to a zero-measure set.
For active sets S_1 and S_2 the KKT condition with respect to {1,2,...,n} is

    ∇ g(x)+A^Ty-μ+ν   =   -v
    
    Ax   =   b,
    
    x_i   =   ℓ_i,  μ_i≥ 0,    i∈ S_1
    
    x_i   =   u_i,  ν_i≥ 0,    i∈ S_2
    
    x_i   ∈   [ℓ_i, u_i], i
    μ_i   =   0,    i∉ S_1
    ν_i   =   0,   i∉ S_2.

We prove that for any S_1, S_2⊆{1,2,...,n}, with probability 0, the strictly complementary condition does not hold.
Since S_1, S_2 have only finitely many choices, we only need to consider fixed S_1 and S_2.
Without loss of generality, assume that x_1=ℓ_1, μ_1=0 and 1∉ S_1.

Consider the Lipschitz continuous map Φ

    Φ(x, y, μ, ν)=(∇ g(x)+A^Ty-μ+ν, Ax)

from the set

    T={(x ,y, μ, ν).   |    x_i=ℓ_i, i∈ S_1,  x_i=u_i, i∈ S_2, 
          ℓ_i≤ x_i≤ u_i, μ_i=0, i ∉ S_1, ν_i=0, i∉ S_2, 
          . μ≽ 0, ν≽ 0, x_1=ℓ_1}

to ℝ^n+m.
Clearly Φ maps from a n+m-1 dimension subset to n+m-dimension space.
Hence, the image is zero-measure in ℝ^n+m.
Consequently, the choice for (v, b) such that the solution exists is of measure zero. Hence, the probability for it is 0 since the probability distribution is continuous.



 §.§ A Proximal Inexact Augmented Lagrangian Multiplier Method

We will state our algorithm in this subsection based on the augmented Lagrangian function.
The Augmented Lagrangian function for (<ref>) is given by

    L(x; y)=f(x)+y^T(Ax-b)+Γ/2Ax-b^2,

where Γ>0 is a constant.
The classical augmented Lagrangian multiplier method minimizes L(x;y) for a fixed y over the box constraint P, and then updates y using the residual of the primal equality constraint Ax=b. Unfortunately, due to the nonconvexity of f, the exact minimization of L(x;y) with respect to x can be difficult. Thus, it is often more practical to minimize L(x;y) inexactly with respect to x. In particular, we recall the following simple inexact augmented Lagrangian multiplier method (which also corresponds to the linearized ADMM algorithm when there is only one primal variable block). 

Though easy to implement numerically, the above inexact augmented Lagrangian multiplier method can behave erratically or even diverge for nonconvex problems (see Figure 1 in Section 6 for a numerical example). To stabilize the convergence behaviour of the inexact augmented Lagrangian multiplier method, we propose a proximal version of the augmented Lagrangian multiplier method. In this new method, we introduce an exponential averaging (or smoothing) scheme to generate an extra sequence {z^t} and insert an extra quadratic proximal term centered at z^t to the augmented Lagrangian function so that the next primal iterate x^t+1 does not deviate too much from the stabilized iterate z^t. More specifically,

let

    K(x, z;y)=L(x; y)+p/2x-z^2,

where p is a positive parameter. Note that the function K is Lipschitz differentiable with modulus L_K=L+p+Γσ^2, and can be made strongly convex in x with modulus γ_K=p+γ>0 if p is chosen to be larger than -γ.
We consider the following proximal inexact augmented Lagrangian multiplier method.









Our main claim is that the introduction of the proximal term can ensure the global convergence of Algorithm <ref> for the nonconvex problem (<ref>).

Suppose Assumption <ref> holds. Moreover, suppose the parameters c and p are selected to satisfing 1/L_K>c>0, p>-γ and that the primal and dual stepsizes β and α to be sufficiently small. Then the dual iterates {y^t} are bounded. Moreover, there holds

    lim_t→∞x^t+1-x^t   =   0,
    lim_t→∞(x^t, X^*)   =   0,
    lim_t →∞(z^t, X^*)   =   0,
    lim_t→∞((x^t, y^t), W)   =   0

and every limit point of the sequence {(x^t, y^t)} generated by Algorithm <ref> is a primal-dual stationary point of (<ref>).

Let

    d(y, z)   =   min_x∈ PK(x, z;y),
    
    x(y, z)   =   _x∈ PK(x, z;y),
    
    M(z)   =   min_x∈ P, Ax=b(f(x)+p/2x-z^2),
    
    x^*(z)   =   _x∈ P, Ax=b(f(x)+p/2x-z^2).

It should be noted that if p>-γ, then f(x)+p/2x-z^2 is strongly convex, so there holds

    K(x,z;y)≥ d(y,z),   M(z)≥ d(y,z),  ∀ y,z,

where the first inequality is due to (<ref>), while the second inequality follows from the strong duality

    M(z)=max_yd(y,z).


Note that the subproblem of the Algorithm <ref> is easy since it involves only projection to the box P.
Compared to Algorithm <ref>, the proximal inexact augmented Lagrangian method (Algorithm <ref>) constructs an auxiliary sequence {z^t} which is a recursive average of the primal sequence {x^t}, and uses it to build a quadratic proximal term in the augmented Lagrangian function L(x^t;y^t). Notice that the recursive averaging step is computationally simple, so Algorithm <ref> has a similar per-computational complexity to Algorithm <ref>. It should be noted that this new quadratic proximal term in the augmented Lagrangian function introduces an extra term p(x^t-z^t) in the gradient of L(x^t;y^t+1). This extra term is an exponentially weighted average of all the previously generated primal iterates {x^0,x^1,...,x^t}. As such, it is different from the well known “momentum term" in the backpropagation training algorithm which is equal to the difference of the previous two primal iterates. Also, this extra term is different from the Nesterov's acceleration scheme which adds to the gradient descent direction a specific (iteration dependent) average of previous two primal iterates.

In the rest of the paper, we fix parameters c<1/L_K and p>-γ.


§ CONVERGENCE ANALYSIS

We will prove Theorem <ref> in this section. The proof consists of a series of lemmas. We will give the proof of the key lemmas and prove the rest in Appendix <ref>, <ref> and <ref>.


 §.§ Key Lemmas

In this subsection, we will give some key lemmas that are needed to establish the main theorem.


  §.§.§ Three Descent Lemmas

In this subsection we give three descent lemmas to estimate the changes in the primal, the dual and the proximal function respectively after one iteration of Algorithm <ref>.
 
For any t, if c<1/L_K,

    K(x^t, z^t; y^t)-K(x^t+1, z^t+1; y^t+1)≥1/2cx^t-x^t+1^2+p/2βz^t-z^t+1^2-αAx^t-b^2.





First, we have the trivial equality:

    K(x^t, z^t; y^t)-K(x^t, z^t; y^t+1)
    =-αAx^t-b^2.

Next, notice that updating x is a standard gradient projection, hence,

    K(x^t, z^t; y^t+1)-K(x^t+1, z^t; y^t+1)
    ≥1/2cx^t-x^t+1^2.

Moreover, recall that in Algorithm <ref>, z^t+1=z^t+β(x^t+1-z^t), we have

    K(x^t+1, z^t; y^t+1)-K(x^t+1, z^t+1; y^t+1)    =   p/2(x^t+1-z^t^2-x^t+1-z^t+1^2) 
       =   p/2(z^t+1-z^t)^T((x^t+1-z^t)+(x^t+1-z^t+1)) 
       =   p/2(2/β-1)z^t-z^t+1^2 
       ≥   p/2βz^t-z^t+1^2,

for β≤ 1.
Combining the above three inequalities yields the desired result.


For any t, we have

    d(y^t+1, z^t+1)-d(y^t, z^t)
    ≥      α(Ax^t-b)^T(Ax(y^t+1, z^t)-b)
          +p/2(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1)).



First, recall that

    K(x, z;y)=f(x)+⟨ y, Ax-b⟩+Γ/2Ax-b^2+p/2x-z^2

so we have

    d(y^t+1, z^t)-d(y^t, z^t)
       =    K(x(y^t+1, z^t), z^t; y^t+1)-K(x(y^t, z^t), z^t; y^t)
       ≥   K(x(y^t+1, z^t), z^t; y^t+1)-K(x(y^t+1, z^t), z^t; y^t)
       =   ⟨ y^t+1-y^t,Ax(y^t+1, z^t)-b⟩,
       =   α(Ax^t-b)^T(Ax(y^t+1, z^t)-b),

where the inequality is because x(y^t, z^t)=min_xK(x, z^t; y^t).
Next, using the same technique, we have

    d(y^t+1,z^t+1)-d(y^t+1, z^t)   =    K(x(y^t+1, z^t+1), z^t+1; y^t+1)-K(x(y^t+1, z^t), z^t; y^t+1)
       ≥    K(x(y^t+1, z^t+1), z^t+1; y^t+1)-K(x(y^t+1, z^t+1), z^t; y^t+1)
       =   p/2(x(y^t+1, z^t+1-z^t+1^2-x(y^t+1, z^t+1)-z^t^2)
       =   p/2(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1)).

Combining these, we get the desired result.


For any t≥ 0, there holds

    M(z^t+1)-M(z^t)≤ p(z^t+1-z^t)^T(z^t-x^*(z^t))+pL̃/2z^t-z^t+1^2,

where L̃=p/p+γ+1.

The proof of Lemma <ref> will be given later since it needs the error bound in Lemma <ref>.

The three terms K(x^t, z^t; y^t), -d(y^t, z^t) and M(z^t) individually do not need to decrease after each iteration; however, some weighted sum of them does! This is the main idea of the proof of Theorem <ref>.
To establish the convergence of Algorithm <ref>, we construct a potential function which decreases sufficiently after each iteration. This potential function is a linear combination of the primal, dual and proximal terms .
Specifically, we will prove that the potential function

    ϕ^t=ϕ(x^t, z^t;y^t)=K(x^t, z^t; y^t)-2d(y^t, z^t)+2M(z^t)

decreases sufficiently after each iteration for sufficiently small α and β, where the functions K, d and P are defined in (<ref>)-(<ref>).
Let f=min_x∈ P, Ax=bf(x). Then M(z)≥f for any z∈ P.
Note that x^t∈ P and z^t∈ P for all t (see the definition of the Algorithm <ref>), so M(z^t) is bounded from below. Moreover, it follows from (<ref>) that

    ϕ^t=(K(x^t, z^t; y^t)-d(y^t, z^t))+(M(z^t)-d(y^t, z^t))+M(z^t)≥ M(z^t)≥f

is also bounded below.
When considering ϕ^t-ϕ^t+1, we have some negative terms, which need to be bounded. To this aim, we need to some so-called ""error bounds”, which bound the iteration sequences by the residuals.


  §.§.§ Error Bounds

In this subsection, we prove some error bounds, which bound the distance of the iterative points to the solution set and bound the difference of x(·, ·), x^*(·) when y and z are perturbed.
The following lemma implies that if the dual residual Ax-b→ 0, the dual variable y is bounded.

Suppose Assumption <ref>(a) holds. If {y^i}⊆ℝ^m be a sequence such that y^i∈ỹ^0+range(A), where ỹ^0 is some fixed vector . If Ax(y^i, z^i)-b→ 0 for some sequence {z^i},
then {y^i} is bounded.


According to Assumption <ref>(a), there exists a positive r>0, such that for any direction d∈ Range(A), we can find an x∈ P satisfying Ax-b=r and Ax-b has the same direction as d.
We claim that if y^i goes to infinity, then Ax(y^i, z^i)-b must be bounded away from 0, i.e., we can not have Ax(y^i, z^i)-b→ 0.
We prove this by contradiction. Assume that y^i→∞ and Ax(y^i, z^i)-b→ 0 for some sequences {y^i} and {z^i}. Since ỹ^0 is fixed, it follows that ỹ^i→∞. By Assumption <ref>(a), there exists a x^i∈ P such that Ax^i-b is of the same direction as -ỹ^i and Ax^i-b=r.
Let

    M=max_x, z∈ P{|f(x)|+p/2x-z^2+Γ/2Ax-b^2+⟨ỹ^0, Ax-b⟩}.

When i is sufficiently large, we have ỹ^i>4M/r and Ax(y^i,z^i)-b<r/2. Therefore, it follows that

    (ỹ^i)^T(Ax^i-b)=-ỹ^i·Ax^i-b=-ỹ^ir<-4M,


    (ỹ^i)^T(Ax^i-b)   =    -rỹ^i
       ≤    -r |(ỹ^i)^T(Ax(y^i, z^i)-b)|/Ax(y^i, z^i)-b
       ≤    -r|(ỹ^i)^T(Ax(y^i, z^i)-b)|/r/2
       ≤   2(ỹ^i)^T(Ax(y^i, z^i)-b),

where the second step follows from Cauchy-Schwartz inequality.
Hence, we have

    K(x^i, z^i; y^i)-K(x(y^i, z^i), z^i; y^i)   ≤    2M+(ỹ^i)^T(Ax^i-b)-(ỹ^i)^T(Ax(y^i, z^i)-b)
       ≤    2M+(ỹ^i)^T(Ax^i-b)-1/2(ỹ^i)^T(Ax^i-b)
       =   2M+1/2(ỹ^i)^T(Ax^i-b)<0,

where the last step is due to (<ref>).
This further implies

    K(x^i, z^i; y^i)<K(x(y^i, z^i), z^i; y^i)

which contradicts the definition of x(y^i, z^i).
Hence, when Ax(y^i, z^i)-b→ 0, we have that y^i is bounded.

Remark Note that each iterate y^t of Algorithm <ref> satisfies the condition of the above lemma with ỹ^0=y^0.


The boundedness of y^t will be used to establish the dual error bound later in Lemma <ref>.
Next result shows that x^*(z) is continuous in z and x(y, z) is Lipschitz-continuous in (y, z).

Suppose p>-γ and z, z'∈ P, then we have

    x^*(z)-x^*(z')≤p/p+γz-z';
          x(y, z)-x(y, z')≤p/p+γz-z';
          x(y, z)-x(y', z)≤p+γ/σy-y'.



We only prove the Lipschitz continuity of x^*(·), the other claims can be proved similarly. Let f(x;z)=f(x)+p/2x-z^2.

    f(x^*(z); z')-f(x^*(z'); z'))
       =   (f(x^*(z); z)-f(x^*(z'); z))-(f(x^*(z'); z')-f(x^*(z'); z))+(f(x^*(z); z')-f(x^*(z); z))
       =   (f(x^*(z); z)-f(x^*(z'); z))-p/2(-2(z'-z)^Tx^*(z')+z'^2-z^2)
          +p/2(-2(z'-z)^Tx^*(z)+z'^2-z^2)
       =   (f(x^*(z); z)-f(x^*(z'); z))+p(z'-z)^T(x^*(z')-x^*(z))
       ≤   -p+γ/2x^*(z)-x^*(z')^2+p(z'-z)^T(x^*(z')-x^*(z)),

where the last inequality is due to the strong convexity of f(x; z) in variable x.
On the other hand, again by the strong convexity, we have

    f(x^*(z); z')-f(x^*(z'); z')≥p+γ/2x^*(z)-x^*(z')^2.

Hence, we have

    -(p+γ)x^*(z)-x^*(z')^2+p(z'-z)^T(x^*(z')-x^*(z))≥ 0,

which by Cauchy-Schwartz inequality further implies

    x^*(z)-x^*(z')≤p/p+γz-z'.

Hence, x^*(·) is Lipschitz continuous with modulus p/(p+γ).


We denote x^+, y^+, z^+ be the updated variables of x, y, z by Algorithm <ref>, namely,

    y^+   =   y+α(Ax-b);
    
    x^+   =   [x-c∇ K(x, z; y^+)]_+;
    
    z^+   =   z+β(x^+-z).

The following simple lemma says that if the algorithm stops, it finds a pair of primal-dual solution.

If

    (x, y, z)=(x^+, y^+, z^+),

Then (x, y)∈ W is a pair of solution.


The proof is just to check the KKT conditions and is omitted since it is trivial.

A couple of corollaries are in order, which are direct from the two lemmas above.

Suppose that Assumption <ref>(a) holds. For any ϵ>0, there exists a δ(ϵ)>0, such that for vectors x,  z∈ P and any y∈ỹ^0+range(A) with a fixed ỹ^0 satisfying

    max{x-x^+, z-x^+, Ax(y^+, z)-b}<δ(ϵ),

we have

    max{(x,X^*), (z, X^*)}<ϵ.


The proof of Corollary <ref> is relegated to Appendix <ref>.

Suppose that Assumption <ref>(a) holds. For any ϵ>0, there exists a δ(ϵ)>0 such that if

    Ax(y^+, z)-b<δ(ϵ)

for some y∈ỹ^0+range(A) with some fixed ỹ^0,
then

    x(y^+, z)-x^*(z)< ϵ,.


The proof of Corollary <ref> can be seen in Appendix <ref>.

Next we develop some primal and dual error bounds. To prove the dual error bound, we need to make use of the Hoffman bound, which can be seen in <cit.>.

Let A∈ℝ^m× n, C∈ℝ^k× n and b∈ℝ^m ,d∈ℝ^k, then the distance from a point x̅∈ℝ^n to the set S={x| Ax≤ b, Cx=d}
is bounded by:

    (x̅, S)^2≤θ^2((Ax̅-b)_+^2+Cx̅-d^2),

where (·)_+ means the projection to the nonegative orthant and θ is a positive constant depending on A and C only.



Suppose p>-γ, Γ>0 are fixed. Then there exist positive constants σ_1,...,σ_4>0 (independent of y^t+1 and z^t) such that the following error bounds hold:

    x^t+1-x^t   ≥   σ_1x^t-x(y^t+1, z^t),
    x^t+1-x^t   ≥   σ_2x^t+1-x(y^t+1, z^t), 
    y-y'   ≥   σ_3x(y,z)-x(y',z), 
    z^t-z^t+1   ≥   σ_4 x^*(z^t)-x^*(z^t+1),
    z^t-z^t+1   ≥   σ_4x(y^t+1,z^t)-x(y^t+1, z^t+1),

where σ_1=cγ_K=c(p+γ), σ_2=σ_1/(1+σ_1), σ_3=γ_K/σ=(γ+p)/σ and  σ_4=γ_K/p=(γ+p)/p.








Furthermore, suppose Assumption <ref>(b) holds, then there exist positive scalars Δ, σ_5 such that

    (y, Y^*(z))   ≤   σ_5Ax(y, z)-b,


if y∈ y^0+range(A), Ax(y, z)-b≤Δ and (z, X^*)≤Δ, where Y^*(z) denotes the solution set of dual multipliers for (<ref>). 








The proof of inequalities (<ref>) and (<ref>)are standard, we prove them in Appendix <ref>. Also note that (<ref>), (<ref>) and (<ref>) are just Lemma <ref>.
It remains to prove the dual error bound (<ref>). To this end, we first write down, for any z, the optimality conditions for the strongly convex proximal optimization problem (<ref>) (note that γ_K=p+γ>0)
as follows:


    ∇ f(x^*(z))+p(x^*(z)-z)+A^Ty^*(z)-μ(z)+ν(z)    =    0, 
    
    Ax^*(z)    =    b, 
    μ_i(z)(x_i-ℓ_i)   =   0,    i=1,⋯, n,
    ν_i(z)(x_i-u_i)   =   0,    i=1,⋯, n, 
    μ_i(z), ν_i(z)   ≥    0,   i=1, ⋯, n,

where y^*(z), μ(z), ν(z) are the Lagrangian multipliers.
For any z, note that y^*(z) only appears in the term A^Ty^*(z). By replacing y^*(z) with its projection to range(A) if necessary, we can assume, without loss of generality, y^*(z)∈range(A).
If z∈ X^*, by strong convexity, the unique optimal solution of (<ref>) is given by x^*(z)=z.
Moreover, for z∈ X^*, the strict complementarity condition

    μ_i(z)=0     ⟹      x_i^*(z)=z_i>ℓ_i, i=1, ⋯, n, 
    ν_i(z)=0     ⟹      x_i^*(z)=z_i<u_i,  i=1, ⋯, n,
 holds.

Now for any y, z, let I(y,z),  I(z) denote the set of inactive inequality constraints (the inequality constraint holds strictly) in problem (<ref>) at x(y,z) and (<ref>) at x^*(z) respectively. We prove that there exists a Δ>0, such that

    I(y, z)=I(z),

We prove it by contradiction. Suppose the contrary, then there is a sequence Δ_k→ 0 and sequences {y_k}⊆ y^0+range(A), {z_k} with

    Ax(y_k, z_k)-b≤Δ_k,  (z_k,X^*)≤Δ_k,

such that I(y_k, z_k) I(z_k), for all k.
Note that y^*(z_k) is the optimal dual solution for the problem (<ref>) with z=z_k, so we have x(y^*(z_k), z_k)=x^*(z_k) and Ax(y^*(z_k), z_k)-b=0, for all k. By Lemma <ref>, we know that {y_k}, {y^*(z_k)} are bounded.
So we can assume that (passing to a subsequence if necessary)

    y_k→ y^*, y^*∈ℝ^m,


    y^*(z_k)→y̅, y̅∈ℝ^m,


    z_k→ z^*, z^*∈ X^*.

According to Lemma <ref>, we have x^*(z^k)→ x^*(z^*) and x(y_k, z_k)→ x(y^*, z^*).
We also have Ax(y^*, z^*)-b=0, hence x(y^*, z^*)=x(z^*).
Since for any i at most one of μ_i and ν_i can be nonzero, {μ_i(z_k)}_k=1^∞, {ν_i(z^k)}_k=1^∞ are bounded by (<ref>) and the fact that y^*(z_k) is bounded. Hence, passing to a subsequence if necessary, we can assume that there exist μ^*, ν^*, such that

    μ_i(z_k)→μ^*,   ν_i(z_k)→ν^*, i.

Hence, (x^*(z^*), y̅, μ^*, ν^*) is a solution to the KKT system (<ref>). Using the strict complementarity property at z^*∈ X^*, we have μ_i, ν_i 0, i∉ I(z^*). When k→∞, we have

    x_i^*(z_k)∈ (ℓ_i, u_i),       i∈ I(z^*),
    μ_i>0,       i∉ I(z^*) x_i^*(z_k)=ℓ_i,
    ν_i>0,        i∉ I(z^*) x_i^*(z_k)=u_i.



This means that I(z_k)=I(z^*). Similarly, we can consider the KKT conditions for (<ref>) and can show that for large enough k there holds x(y_k,z_k)→ x(y^*,z^*)=x^*(z^*), and I(y_k,z_k)=I(y^*,z^*)=I(z^*). This implies I(y_k,z_k)=I(z^*)=I(z_k) for large k, a contradiction.
So next we assume that I(y, z)=I(z). Also define A_1={i| x_i^*(z)=x_i(y, z)=ℓ_i} and A_2={i| x_i^*(z)=x_i(y, z)=u_i}.
Recall the optimization problem (<ref>) for function K:

    [                                   f(x)+p/2x-z^2+y^T(Ax-b)+Γ/2Ax-b^2;                                          ℓ_i≤ x_i≤ u_i,  i=1, ⋯, n, ]

and its KKT conditions:

    ∇ f(x(y, z))+p(x(y, z)-z)+A^Ty+Γ A^T(Ax(y,z)-b)-μ(y, z)+      ν(y,z) = 0, 
    μ_i(y, z)≥ 0, ℓ_i-x_i(y,z)=0,          i∈ A_1,
    ν_i(y, z)≥ 0,  x_i(y, z)-u_i=0,         i∈ A_2,
    μ_i(y, z)=0, ν_i(y, z)=0,         i∉A_1∪ A_2,
    
    x_i(y,z)∈ [ℓ_i, u_i],         i∈ I(z).


Notice that the optimality conditions (<ref>) for the convex proximal problem (<ref>) can be rewritten (using the information about the active and inactive sets) as:


    ∇ f(x^*(z))+p(x^*(z)-z)+A^Ty^*(z)-μ(z)+ν(z) = 0,        
    
    Ax^*(z)-b=0,        
    μ_i(z)≥ 0, ℓ_i-x_i^*(z)=0,         i∈ A_1,
    ν_i(z)≥ 0,  (x_i^*(z)-u_i)=0,         i∈ A_2,
    μ_i(z), ν_i(z)≥ 0,  x_i^*(z)∈[ℓ_1, u_1],         i∈ I(z).

Let g(x; z)=f(x)+p/2x-z^2+Γ/2Ax-b^2. Consider the following linear system.

    A^Ty^*(z)-μ(z)+ν(z) = -∇ g(x^*(z); z),        
    μ_i(z)≥ 0, ℓ_i-x_i^*(z)=0,         i∈ A_1,
    ν_i(z)≥ 0,  (x_i^*(z)-u_i)=0,         i∈ A_2,
    μ_i(z), ν_i(z)≥ 0,  x_i^*(z)∈[ℓ_1, u_1],         i∈ I(z).

Notice that in the above system, the left-hand-side is linear of y(z), μ(z), ν(z).
By (<ref>), the vector (x(y,z),y,μ(y,z),ν(y,z)) satisfies (<ref>) approximately.
Using Hoffman bound (Proposition <ref>) to (<ref>) at the point (y,μ(y,z),ν(y,z)), we obtain

    (y, Y^*(z))^2   ≤   min_(x^*(z), y(z),μ(z),ν(z))
    y-y(z)^2+(μ(y,z),ν(y,z))-(μ(z),ν(z))^2
    
       ≤    2θ^2∇ g(x(y, z); z)-∇ g(x^*(z); z)^2.
     [3pt]
       ≤    2θ^2((L+p+Γσ^2)^2x(y, z)-x^*(z)^2),

where σ is the largest singular value of A. Using Assumption <ref>(d) and the strong convexity of g(·; z), we obtain

    (y, Y   ^*(   z))^2≤ 2θ^2(L+p+Γσ^2)^2/p+γ⟨∇ f(x(y,z))-∇ f(x^*(z))+p(x(y,z)-x^*(z)), 
          x(y, z)-x^*(z)⟩.
       =   -2θ^2 (L+p+Γσ^2)^2/p+γ⟨ y-y(z),A(x(y, z)-x^*(z))⟩
          +2θ^2 (L+p+Γσ^2)^2/p+γ⟨ (μ(y, z)-μ(z))-(ν(y, z)-ν(z)), x(y, z)-x^*(z)⟩
       ≤    2θ^2((L+p+Γσ^2)^2/p+γAx(y, z)-by-y(z)),

where the equality follows from (<ref>)-(<ref>) and the cross term (<ref>) vanishes because we have

    [ μ_i(y, z)(x_i(y, z)-ℓ_i)=0,     μ_i(z)(x_i^*(z)-ℓ_i)=0,                i=1,2,...,n,;  μ_i(y, z)(x_i^*(z)-ℓ_i)≥0,   μ_i(z)(x_i(y, z)-ℓ_i)≥ 0,                i=1,2,...,n,; ν_i(y, z)(u_i-x_i(y, z))=0,     ν_i(z)(u_i-x_i^*(z))=0,                i=1,2,...,n,;  ν_i(y, z)(u_i-x_i^*(z))≥0,    ν_i(z)(u_i-x_i(y, z))≥0,                i=1,2,...,n. ]

Hence we finish the proof of the dual error bound (<ref>).

Remark. Notice that the error bounds (<ref>), (<ref>) hold for all strongly convex problems over a convex set when a linear term in the objective function is perturbed.
We now give the proof of Lemma <ref>.


First, using Danskin's theorem in convex analysis <cit.>, we have:

    ∇ M(z^t)=p(z^t-x^*(z^t)),

where M(z) is defined by (<ref>).
So it suffices to prove that

    ∇ M(z^t)-∇ M(z^t+1)≤ (σ^-1_4+1)z^t-z^t+1.

But this is a direct corollary of the error bound (<ref>) in Lemma <ref>. The proof is complete.



 §.§ Proof of Theorem <ref>


Using the three descent lemmas in Subsection <ref>, we get

    ϕ^t-ϕ^t+1
    
       ≥    (1/2cx^t+1-x^t^2-αAx^t-b^2+p/2βz^t-z^t+1^2)
          +2 (α(Ax^t-b)^T(Ax(y^t+1, z^t)-b)+p(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1)) )
          +2 (p(z^t+1-z^t)^T(x^*(z^t)-z^t)-pL̃/2z^t+1-z^t^2 )
       =    (1/2cx^t+1-x^t^2-αAx^t-b^2+p/2βz^t-z^t+1^2)
    
    
    +2α(Ax^t-b)^T(Ax(y^t+1, z^t)-b)
          +p(z^t+1-z^t)^T((z^t+1-z^t)-2(x(y^t+1, z^t+1)-x^*(z^t)))
    
    
    -pL̃/2z^t-z^t+1^2
       =   (1/2cx^t+1-x^t^2-αAx^t-b^2+p/2βz^t-z^t+1^2)
    +2α(Ax^t-b)^T(Ax(y^t+1, z^t)-b)
          +p(z^t+1-z^t)^T((z^t+1-z^t)-2(x(y^t+1, z^t+1)-x(y^t+1, z^t))-2(x(y^t+1, z^t)-x^*(z^t)))
          -pL̃/2z^t-z^t+1^2.

Let λ be an arbitrary positive scalar, and by the fact that

    (z^t+1-z^t)/λ-λ(x(y^t+1, z^t)-x^*(z^t))^2≥ 0,

we have

    -2(z^t+1-z^t)^T(x(y^t+1, z^t)-x^*(z^t))≥ -z^t-z^t+1^2/λ-λx(y^t+1, z^t)-x^*(z^t)^2.

Using Cauchy-Schwarz inequality and the error bound (<ref>) in Lemma <ref>, we have

    -2(z^t+1-z^t)^T(x(y^t+1, z^t+1)-x(y^t+1, z^t))   ≥    -z^t-z^t+1x(y^t+1, z^t+1)-x(y^t+1, z^t)
       ≥    -1/σ_4z^t-z^t+1^2.

Substituting these two inequalities into (<ref>), we have

    ϕ^t-ϕ^t+1
       ≥   1/2cx^t+1-x^t^2-(αAx^t-b^2-2α(Ax^t-b)^T(Ax(y^t+1, z^t)-b)+αAx(y^t+1, z^t)-b^2)
          +αAx(y^t+1, z^t)-b^2
    
    +(p/2β+p-p/σ_4-p/λ-pL̃/2)z^t-z^t+1^2
          -pλx(y^t+1, z^t)-x^*(z^t)^2.

By completing the square, we further obtain

    ϕ^t-ϕ^t+1   ≥   1/2cx^t+1-x^t^2-αA(x(y^t+1, z^t)-x^t)^2
    
    
    +αAx(y^t+1, z^t)-b^2
          +(p/2β+p-p/σ_4-p/λ-pL̃/2)z^t-z^t+1^2
    
    
    -pλx(y^t+1, z^t)-x^*(z^t)^2
       ≥    (1/2c-ασ^2/σ_1^2)x^t-x^t+1^2+αAx(y^t+1, z^t)-b^2
    
          
    +(p/2β+p-p/σ_4-p/λ-pL̃/2)z^t-z^t+1^2
    
    
    -pλx(y^t+1, z^t)-x^*(z^t)^2,

where σ is the spectral norm of A, λ is any positive scalar, and the last step is due to the error bound (<ref>) in Lemma <ref>.
Let λ=Dβ for some sufficiently large D (for example, D>6), and set β sufficiently small (β≤β') for some constant β'>0, such that

    p/2β+p-p/σ_4-p/λ-pL̃/2≥p/3β.

Therefore, if we choose α< σ^2_1/4cσ^2, 
then it follows from (<ref>) that


    ϕ^t-ϕ^t+1   ≥   1/4cx^t-x^t+1^2+αAx(y^t+1, z^t)-b^2
          +p/3βz^t-z^t+1^2-pDβx(y^t+1, z^t)-x^*(z^t)^2.

It remains to bound the term pDβx(y^t+1, z^t)-x^*(z^t)^2 in the above expression. The main proof idea is as follows.
In view of the dual error bound (<ref>), when the residuals are sufficiently small, we can use the dual residual Ax(y^t+1, z^t)-b to bound (y^t+1, Y^*(z^t)) and then further use the error bound (<ref>) to bound x(y^t+1, z^t)-x^*(z^t). When some residual is not too small, we will make use of the compactness of the feasible set to bound the term x(y^t+1, z^t)-x^*(z^t)^2.
Define

    M=max_x_1, x_2∈ Px_1-x_2,  ζ=min{Δ, δ(Δ/√(6D))}

and set

    β<min{β', ζ^2/8cpDM^2, ζ^2α/2pDM^2, ασ_3^2/2pDσ_5^2},

where Δ is defined in Lemma <ref> and δ(·) is defined in Corollary <ref> and Corollary <ref>.
We also define the following three conditions:

    x^t-x^t+1^2≤ 8cpDM^2β, 
    Ax(y^t+1, z^t)-b^2≤2pDM^2/αβ,
    x^t+1-z^t^2=(z^t-z^t+1)/β^2≤ 6Dx(y^t+1, z^t)-x^*(z^t)^2.

We now consider two cases.

Case 1. Conditions (<ref>)-(<ref>) hold. In this case, it follows from (<ref>)-(<ref>) that

    x^t-x^t+1≤ζ = min{Δ, δ(Δ/√(6D))}≤Δ,
    Ax(y^t+1, z^t)-b≤ζ = min{Δ, δ(Δ/√(6D))}≤Δ,

which further implies

    x^t+1-z^t=(z^t-z^t+1)/β≤√(6D)x(y^t+1, z^t)-x^*(z^t)≤√(6D)Δ/√(6D)=Δ,

where the last inequality follows from Corollary <ref>.
Therefore, the error bounds (<ref>) and (<ref>) in Lemma <ref> hold and we have


    pDβx(y^t+1, z^t)-x^*(z^t)^2
       ≤    pDβσ_3^-2·(y^t+1, Y^*(z^t))^2 
       ≤    pDβσ_3^-2σ_5^2Ax(y^t+1, z^t)-b^2 
       ≤   α/2Ax(y^t+1, z^t)-b^2,

where the last step follows from (<ref>).
It then follows from (<ref>) that

    ϕ^t-ϕ^t+1   ≥   1/4cx^t-x^t+1^2+α/2Ax(y^t+1, z^t)-b^2+p/3βz^t-z^t+1^2.

Case 2. One of the conditions (<ref>)-(<ref>) is violated. Consider three subcases:


   Case 2.1. x^t-x^t+1^2≥ 8cpDβ M^2. In this case, we have

    1/4cx^t-x^t+1   ^2   -pDβx(y^t+1, z^t)-x^*(z^t)^2 
       ≥   1/8cx^t-x^t+1^2+1/8c8cpDβ M^2-pDβx(y^t+1, z^t)-x^*(z^t)^2 
       ≥   1/8cx^t-x^t+1^2+pDβ M^2-pDβ M^2 
       =    1/8cx^t-x^t+1^2,

where the second step is due to (<ref>).

   Case 2.2. Ax(y^t+1, z^t)-b^2≥2pD β/αM^2. In this case, we use this condition and (<ref>) to obtain

    αAx(y^t+1, z^t)-b^2   -    pDβx(y^t+1, z^t)-x^*(z^t)^2 
       ≥   α/2Ax(y^t+1, z^t)-b^2+α/2·2pD β/αM^2-pDβ M^2 
    
       =   α/2Ax(y^t+1, z^t)-b^2.



  Case 2.3. (z^t-z^t+1)/β^2≥ 6Dx(y^t+1, z^t)-x^*(z^t)^2. In this case, we have

    p/3βz^t-z^t+1^2   -    pDβx(y^t+1, z^t)-x^*(z^t)^2 
       ≥   p/3βz^t-z^t+1^2-pDβx(y^t+1, z^t)-x^*(z^t)^2 
       ≥   p/6βz^t-z^t+1^2+p/6β6β^2Dx(y^t+1, z^t)-x^*(z^t)^2 
          -pDβx(y^t+1, z^t)-x^*(z^t)^2 
       =    p/6βz^t-z^t+1^2.



Considering (<ref>), we have in all three subcases:

    ϕ^t-ϕ^t+1≥1/8cx^t-x^t+1^2+α/2Ax(y^t+1, z^t)-b^2+p/6βz^t-z^t+1^2.


Combining (<ref>), (<ref>), (<ref>) and (<ref>) yields



    pDβx(y^t+1, z^t)-x^*(z^t)^2≤1/8cx^t-x^t+1^2+α/2Ax(y^t+1, z^t)-b^2+p/6βz^t-z^t+1^2.

Then according to (<ref>), we have

    ϕ^t-ϕ^t+1≥1/8cx^t-x^t+1^2+α/2Ax(y^t+1, z^t)-b^2+p/6βz^t-z^t+1^2,  ∀  t≥0.

Since ϕ^t is bounded below, we must have

    max{x^t+1-x^t, Ax(y^t+1, z^t)-b, z^t-x(y^t+1, z^t)}→ 0.

This together with Corollary <ref> shows that the KKT condition for (<ref>) is satisfied in the limit. This completes the proof.


Notice that we can take D=7, then the threshold for the stepsize β is given in (<ref>), where ζ and M are defined in (<ref>). Since ζ is defined by constants Δ, δ(Δ) from the  “local dual error bound”, it is difficult to derive an explicit bound for the stepsize β. In practice, we will need to tune the stepsize β to ensure fast convergence. To give an explicit upper-bound for β, a “global” error bound is needed. This will be the focus of a new paper under planning.

Theorem <ref> establishes the global convergence of Algorithm <ref> to a stationary solution. However, it does not
address the question of convergence rate. The latter is considered in the next section. We will show that the iteration complexity is O(1/ϵ^2) for a fixed problem and
for the special case when the objective function is (nonconvex) quadratic the convergence rate is linear. 





§ ITERATION COMPLEXITY AND LINEAR CONVERGENCE



 §.§ O(1/ϵ^2) iteration complexity

In this subsection, we will see that the iteration complexity of Algorithm <ref> is O(1/ϵ^2) to attain an ϵ-stationary solution.
We first define the ϵ-stationary solution as in <cit.>. Let ι(x) be the indicate function of the set P, i.e. , ι(x)=0 if x∈ P and ι(x)=∞ otherwise.

We say that (x, y) is an ϵ-solution of (<ref>) if Ax-b≤ϵ and  there exists a vector  v∈∇ f(x)+A^Ty+∂ι(x) with v≤ϵ.


There exists a constant B>0 such that  for any t>0, we can find an s∈{0, 1, ⋯, t-1} such that (x^s+1, y^s+1) is a B/√(t)-solution. In other word, we can find an ϵ- solution within B^2/ϵ^2 iterations.

Remark.Note that the constant B depends on some “local” constants and error bounds related to the local geometry of the problem.  Global error bounds are needed for a global complexity analysis, which depends on global constants. This interesting topic will be left as future work in our follow-up paper.

According to (<ref>) , for any t>0, there exists an s∈{0, 1, ⋯, t-1} such thatϕ^s-ϕ^s+1≤ϕ^0-f/t.  Then by inequality (<ref>), letting

    C=(ϕ^0-f)·max{8c, 2/α, 6β/p},

 we have

    x^s-x^s+1^2   <   C/t;
    Ax(y^s+1, z^s)-b^2   <   C/t;
    x^s+1-z^s^2   <   C/t,

According to Algorithm <ref>, we have

    x^s+1= min_x{⟨∇_xK(x^s, z^s; y^s+1), x-x^s ⟩+1/cx-x^s^2+ι(x)}.

Hence, by the optimality condition,

    0∈∇_xK(x^s, z^s;  y^s+1)+2/c(x^s+1-x^s)+∂ι(x^s+1).

Letting

    v=∇_xK(x^s+1, z^s; y^s+1)-∇_xK(x^s, z^s; y^s+1)-2/c(x^s+1-x^s)-Γ(Ax^s+1-b)-p(x^s+1-z^s),

we have

    v∈∇_x(f(x^s+1+A^Ty^s+1))+∂(ι(x^s+1)).

By inequalities (<ref>) and (<ref>), we have

    Ax^s+1-b   ≤   Ax(y^s+1, z^s)-b+A(x^s+1-x(y^s+1, z^s))
       ≤   √(C)/√(t)+σσ_2^-1√(C)/√(t)
       ≤   √(C)(1+σσ_2^-1)/√(t),

where the first inequality is because of triangle inequality. 
Then we have

    v   ≤    (L+p+Γσ^2+2/c)x^s-x^s+1+ΓAx^s+1-b+px^s+1-z^s
       ≤    (L+p+Γσ^2+2/c)√(C)/√(t)+Γ(√(C)/√(t)+√(C)σσ_2^-1/√(t))+p√(C)/√(t)
       ≤   √(B)/√(t),

where the first inequality is due to the triangle inequality and the Lipschitz continuity of ∇_xK(·, z; y), the second inequality is because of inequalities (<ref>), (<ref>) and

    B=((L+p+Γσ^2+2/c)+Γ(1+σσ_2^-1)+p)^2C.

Notice that B>0. Then the result holds for v and B and (x^s+1, y^s+1) is a B/√(t)-solution.




 §.§ Linear convergence for quadratic programming

In this subsection, we consider a nonconvex quadratic program (QP), which is a special case of (<ref>) with f(x) being a quadratic function

    f(x)=1/2x^TQx+r^Tx.

We will strengthen Theorem <ref> in this case by showing that Algorithm <ref> converges linearly to a stationary point of the nonconvex QP problem.
By Theorem <ref>, we have (x^t,X^*)→ 0, (z^t, X^*)→ 0 and x^t+1-x^t→ 0 as t→∞. Since X^* is the union a finite number of polyhedral sets, it follows that the connected components of X^* are properly separated in the sense that there is a positive distance between each pair of distinct connected components of X^*. As a result, the sequences {x^t} and {z^t} will both converge to one unique connected component of X^*. Moreover, it is known <cit.> that for a quadratic programming problem, the objective function value f(x) is constant on each of the connected component of X^*. Let f^* denote the value of f(x) over the connected component of X^* to which x^t (and z^t) converges. Then f(x^t)→ f^* as t→∞.
Note that by the analysis above, when t is sufficiently large, f(z̅^t)=f^* and z̅^t belongs to the connected component that x^t and z^t converge to.
We summarize the above analysis as follows.

Assume the parameters of Algorithm <ref> are chosen to guarantee its convergence (cf. Theorem <ref>).
We have the following.


  * The quadratic cost function f is constant on each of the connected components of X^*.

  * The two sequences {x^t}, {z^t} converge to a same connected component of X^*.

  * Let z̅^t=_z̅∈ X^*=(z^t,X^*) and x^* is any limit point of {x^t}, then f(x^*)=f(z̅^t)=f^*, where f^* is a constant, for all sufficiently large t.


  * There exists a constant ϵ(x^0, y^0, z^0), such that for any z^t, if (z^t, X^*)< ϵ(x^0, y^0, z^0), we always have

    f(z̅^t)=f^*.



Denote

    Φ_p^t(z)   =   K(x^t, z; y^t)-d(y^t, z)
    Φ_d^t(z)   =   M(z)-d(y^t, z) 
    Φ_prx^t(z)   =   M(z)-f^*
    z̅^t   =   _x∈ X^*x-z^t.
    Φ^t(z)   =   Φ^t_p(z)+Φ^t_d(z)+Φ^t_prx(z).

Remark. Notice that Φ^t(z^t) decreases monotonically due to (<ref>) and converges to 0 by Claim <ref>, we have Φ^t(z^t)=ϕ^t-f^* and Φ^t(z^t)≥ 0 for any t≥ 0.
Also notice that

    ϕ^t-ϕ^t+1=Φ^t(z^t)-Φ^t+1(z^t+1).

To prove linear convergence, we make use of some “cost-to-go” estimates from <cit.>.


There exist constants τ_1, τ_2, τ_3>0 such that



    K(x^t+1, z^t; y^t+1)-d(y^t+1, z^t)   ≤   τ_1x^t-x^t+1^2,
    
    
    
    M(z^t)-d(y^t+1, z^t)   ≤   τ_2Ax(y^t+1, z^t)-b^2,
    
    
    
    M(z^t)-f^*   ≤   τ_3z^t-x^*(z^t)^2.



The proof of (<ref>) and (<ref>) is simply to combine Lemma 3.1 of <cit.> with (<ref>), (<ref>) in Lemma <ref>. Specifically, we only need to replace L(x; y), d(y), p^* of <cit.> by K(x, z;y), d(y, z), M(z) respectively with z fixed, since the error estimates of <cit.> are independent of the linear term in K(x,z;y).
To prove the estimate (<ref>), we first notice that f has a Lipschitz continuous gradient and that the classic proximal algorithm belongs to the class of approximate gradient projection algorithm (see <cit.>), namely, we have

    x^*(z^t)=[z^t-∇ M(z^t)+Θ(t)]_+,

where Θ(t) satisfies Θ(t)≤ηx^*(z^t)-z^t for some η>0. Note that the above are similar to the inequalities (3.5) and (3.7) in the proof of the second part of Lemma 3.1 of <cit.>.
Therefore, similar to <cit.>, there exists a constant τ' such that

    f(x^*(z^t))-f(z̅^t)
    ≤τ'(z^t-x^*(z^t)^2+z^t-z̅^t^2),

where z̅^t is the projection of z^t to stationary solution set of problem (<ref>). According to Claim <ref>, z̅^t is in the connected component that {x^t}, {z^t} converge to and f(z̅^t)=f^*. 
From Theorem 2.1 in <cit.>, there exists a constant τ̅>0, such that

    z^t-z̅^t≤τ̅z^t-x^*(z^t).

Combining the above two inequalities and using the definition of M(z), we have

    M(z^t)-f^*=f(x^*(z^t))+p/2x^*(z^t)-z^t^2-f(z̅^t)
    ≤((τ̅^2+1)τ'+p/2)z^t-x^*(z^t)^2.

This completes the proof of (<ref>) by setting τ_3=((τ̅^2+1)τ'+p/2).

The following is the error bound for proximal algorithm, which can be seen in <cit.>.

There exist a constant τ̅ such that

    (z, x^*)≤τ̅z-x^*(z).



For any t≥ 0, we have

    z^t-x^*(z^t)^2≤ B_1(β)z^t-z^t+1^2+B_2(β)x^t-x^t+1^2+B_3(β)Ax(y^t+1, z^t)-b^2,

where

    B_1(β)=3(1/β^2+1/6Dβ^2),   B_2(β)=3(1/σ_2^2+1/8cpDβ)   B_3(β)=3α/2pDβ.



Notice that

    z^t-x^*(z^t)^2
       =   (z^t-x^t+1)+(x^t+1-x(y^t+1, z^t))+(x(y^t+1, z^t)-x^*(z^t))^2
       =   (z^t-z^t+1)/β+(x^t+1-x(y^t+1, z^t))+(x(y^t+1, z^t)-x^*(z^t))^2
       ≤   3(z^t-z^t+1^2/β^2+x^t+1-x(y^t+1, z^t)^2+x(y^t+1, z^t)-x^*(z^t)^2)
       ≤   3(z^t-z^t+1/β^2+x^t-x^t+1^2/σ_2^2+x(y^t+1, z^t)-x^*(z^t)^2),

where (i) is  because z^t+1=z^t+β(x^t+1-z^t), (ii) is due to the convexity of the square norm and (iii) is due to (<ref>).
On the other hand, by (<ref>), we have

    x(y^t+1, z^t)-x^*(z^t)^2≤1/pDβ(1/8cx^t-x^t+1^2+α/2Ax(y^t+1, z^t)-b^2+p/6βz^t-z^t+1^2).

Finally, we substitute (<ref>) to (<ref>) and get the desired result.


For t>0 and β sufficiently small, we have

    ϕ^t-ϕ^t+1≥ C_1(β)x^t-x^t+1^2+C_2(β)Ax(y^t+1, z^t)-b^2+C_3(β)z^t-x^*(z^t)^2+C_4(β)z^t-z^t+1^2,

where C_i(β)>0, i=1, 2, 3, 4 are constants depending on β.


By Lemma <ref>, we have

    z^t-x^*(z^t)^2≤ B_1(β)z^t-z^t+1^2+B_2(β)x^t-x^t+1^2+B_3(β)Ax(y^t+1, z^t)-b^2.

Let B̅=max{B_1(β)/1/8c, B_2(β)/α/2, B_3(β)/p/6β}. We then have

    z^t-x^*(z^t)^2/(2B̅)   ≤   1/16cx^t-x^t+1^2+α/4Ax(y^t+1, z^t)-b^2+p/12βz^t-z^t+1^2
       ≤   1/2(ϕ^t-ϕ^t+1).

Hence, we have

    ϕ^t-ϕ^t+1   ≥   1/8cx^t-x^t+1^2+α/2Ax(y^t+1, z^t)-b^2+p/6βz^t-z^t+1^2
       ≥   1/16cx^t-x^t+1^2+α/4Ax(y^t+1, z^t)-b^2+p/12βz^t-z^t+1^2+z^t-x^*(z^t)^2/(2B̅).

This completes the proof.

Now we can prove the linear convergence of Algorithm <ref> for quadratic programming.


If f(z̅^t)=f^* and max{(z^t, X^*), Ax(y^t+1, z^t)-b}<Δ, we have

    Φ^t(z^t)-Φ^t+1(z^t+1)≥κΦ^t+1(z^t+1),

where Δ is defined in Lemma <ref> and κ>0 is a constant. Hence, Φ^t+1(z^t+1)<1/1+κΦ^t(z^t).




We need to relate the gap Φ^t(z^t)=Φ_p^t(z^t)+Φ_d^t(z^t)+Φ_prx^t(z^t) to the estimates in Lemma <ref>.
Then using Lemmas <ref>-<ref> and Lemma <ref>, we have

    Φ_p^t+1(z^t)   ≤   τ_1x^t-x^t+1^2,
    Φ_d^t+1(z^t)   ≤   τ_2Ax(y^t+1, z^t)-b^2,
    Φ_prx^t+1(z^t)   ≤   τ_3z^t-x^*(z^t)^2.

Recall the inequality (<ref>) in the proof of Lemma <ref>

    K(x^t+1, z^t+1;y^t+1)-K(x^t+1, z^t; y^t+1)≤ -p/2βz^t-z^t+1^2,

and the inequality (<ref>)

    d(y^t+1, z^t)-d(y^t+1, z^t+1)≤ -p/2(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1)).

Then for Φ_p^t+1(z^t), we have

    Φ_p^t+1(z^t+1)-Φ_p^t+1(z^t)   =   ( K(x^t+1, z^t+1;y^t+1)-K(x^t+1, z^t; y^t+1))+( d(y^t+1, z^t)-d(y^t+1, z^t+1))
       ≤    -p/2βz^t-z^t+1^2-p/2(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1)).

For Φ_d^t+1(z^t), we have

    Φ_d^t+1(z^t+1)-Φ_d^t+1(z^t)   =   d(y^t+1, z^t)-d(y^t+1, z^t+1)+M(z^t+1)-M(z^t)
       ≤    -p/2(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1))+M(z^t+1)-M(z^t).


Also, it follows from the definition (<ref>) that

    Φ_prx^t+1(z^t+1)=Φ_prx^t(z^t)+M(z^t+1)-M(z^t).

Combining the above three inequalities, we have

    Φ_p^t+1(z^t+1)   +   Φ_d^t+1(z^t+1)+Φ_prx^t+1(z^t+1)
       ≤   Φ_p^t+1(z^t)+Φ_d^t+1(z^t)+Φ_prx^t+1(z^t)
          -p/2βz^t-z^t+1^2-p/2(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1))
          -p/2(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1))+2(M(z^t+1)-M(z^t))
       =   Φ_p^t+1(z^t)+Φ_d^t+1(z^t)+Φ_prx^t+1(z^t)
          -p/2βz^t-z^t+1^2-p(z^t+1-z^t)^T(z^t+1+z^t-2x(y^t+1, z^t+1))
          +2p(z^t+1-z^t)^T(z^t-x^*(z^t))+p/σ_4z^t-z^t+1^2
       =   Φ_p^t+1(z^t)+Φ_d^t+1(z^t)+Φ_prx^t+1(z^t)
          +(-p/2β+p/σ_4-p)z^t-z^t+1^2+2p(z^t+1-z^t)^T(x(y^t+1, z^t+1)-x^*(z^t)),

where the second step follows from (<ref>).


















Using the inequality 2a^Tb≤a^2+b^2 for any a,b∈ℝ^n, we have

    Φ_p^t+1(z^t+1)   +   Φ_d^t+1(z^t+1)+Φ_prx^t+1(z^t+1)
       ≤   Φ_p^t+1(z^t)+Φ_d^t+1(z^t)+Φ_prx^t+1(z^t)
          +(-p/2β+p/σ_4-p)z^t-z^t+1^2+pz^t-z^t+1^2+p
    x(y^t+1, z^t+1)-x^*(z^t)^2
       =   Φ_p^t+1(z^t)+Φ_d^t+1(z^t)+Φ_prx^t+1(z^t)
          +(-p/2β+p/σ_4)z^t-z^t+1^2+p
    x(y^t+1, z^t+1)-x^*(z^t)^2
       ≤   Φ_p^t+1(z^t)+Φ_d^t+1(z^t)+Φ_prx^t+1(z^t)
          +(-p/2β+p/σ_4)z^t-z^t+1^2+pσ_5^2/σ_3^2Ax(y^t+1, z^t)-b^2,

where the last step is due to the error bounds (<ref>), (<ref>) in Lemma <ref>.
Let

    D_1(β)=τ_1,   D_2(β)=τ_2+pσ_5^2/σ_3^2,   D_3(β)=τ_3,   D_4(β)=-p/2β+p/σ_4,

then we have

    Φ^t+1(z^t+1)    =   Φ_p^t+1(z^t+1)+Φ_d^t+1(z^t+1)+Φ_prx^t+1(z^t+1)
    
    
       ≤   Φ_p^t+1(z^t)+Φ_d^t+1(z^t)+Φ_prx^t+1(z^t)
          +(-p/2β+p/σ_4)z^t-z^t+1^2+pσ_5^2/σ_3^2Ax(y^t+1, z^t)-b^2
    
    
    
       ≤    D_1(β)x^t-x^t+1^2+D_2(β)Ax(y^t+1, z^t)-b^2
          +D_3(β)z^t-x^*(z^t)^2+D_4(β)z^t-z^t+1^2,

where the last step follows from (<ref>)-(<ref>).

Recall the definition of potential function ϕ (cf. (<ref>)). It follows that

    ϕ^t-ϕ^t+1   =   (Φ_p^t+1(z^t)+Φ_d^t+1(z^t)+Φ_prx^t+1(z^t))
    -(Φ_p^t+1(z^t+1)+Φ_d^t+1(z^t+1)+Φ_prx^t+1(z^t+1))
       =   Φ^t(z^t)-Φ^t+1(z^t+1)
       ≥   C_1(β)x^t-x^t+1^2+C_2(β)Ax(y^t+1, z^t)-b^2
          
    +C_3(β)z^t-x^*(z^t)^2+C_4(β)z^t-z^t+1^2 
    
       ≥   κ(D_1(β)x^t-x^t+1^2+D_2(β)Ax(y^t+1, z^t)-b^2.
          .+D_3(β)z^t-x^*(z^t)^2+D_4(β)z^t-z^t+1^2)
       ≥   κΦ^t+1(z^t+1) ,

where κ:=min_i: D_i(β)>0{C_i(β)/D_i(β)} and the first inequality follows from Lemma <ref>. This completes the proof.








































































Combining the above analysis, we can prove the following theorem.

There exists a t_max>0 such that

    Φ^t(z^t)-Φ^t+1(z^t+1)≥κΦ^t+1(z^t+1),  t≥ t_max.



By Claim <ref>,  we have

    f(z̅^t)=f^*, .

Moreover, according to Theorem <ref>, when t is sufficiently large, we have

    (z^t, X^*)   <   Δ,
    Ax(y^t+1, z^t)-b   <   Δ.

Hence there exists a t_max such that the conditions of Lemma <ref> hold after t≥ t_max. The linear convergence follows immediately from Lemma <ref>.


Remark. Note that it follows from Lemma <ref>, (<ref>) and Theorem <ref> that the sequence {z^t} converges R-linearly  after t≥ t_max.
The constant t_max depends on the problem structure and the initial point. It can be shown that

    t_max=max{Φ^0(z^0)τ̅^2B̅/min{Δ, ϵ(x^0, y^0, z^0)}, Φ^0(z^0)/Δ^2α/2}·max{1, 1/κ}.

However, this estimate depends on some unknown constants such as Δ and ϵ(x^0,y^0,z^0) which are implicitly defined by the problem structure and geometry. A global error bound is needed to provide an explicit estimate of t_max. This will be the focus of a forthcoming paper.



§ MULTI-BLOCK CASE AND A LINEARIZED PROXIMAL ADMM

In this section we consider the multi-block case. Consider the following optimization problem:

    [                           f(x_1, x_2, ⋯, x_k);                  [5pt] ∑_j=1^kA_jx_j=b, x∈ P, ]

where x_j∈ℝ^n_j for all 1≤ j≤ k, x=(x_1, x_2, ⋯, x_k)∈ℝ^n and P=∏_i=1^n[ℓ_i, u_i].
We also denote A=[A_1, A_2,⋯, A_k] and σ̅=max_1≤ i≤ kA_i_2, where ·_2 is the spectral norm of a matrix.
We still adopt Assumption <ref> and hence every partial gradient ∇_x_if(·) is L-Lipschitz-continuous.
To solve this problem, we use the following linearized proximal ADMM <ref>, which updates the primal variables blockwise. Let us denote

    x^t(j)=(x_1^t+1, x_2^t+1, ⋯, x_j-1^t+1, x_j^t, ⋯, x_k^t).



Here [·]_+^j means the projection to the set ∏_i∈ N_j[ℓ_i, u_i], where N_j is the set of indexes for the j-th variable block x_j. Note that here we take the stepsize 
    c<1/L+p+Γσ̅^2.
 
To prove the convergence of Algorithm <ref>, we can follow the same line as that for the proof of Theorem <ref>. The only differences are the proof for (<ref>)-(<ref>) in Lemma <ref> and the proof of Lemma <ref>.
For the primal error bounds (<ref>)-(<ref>), we only give the proof of the first one and the others can be proved using the same techniques as that in the proof of Lemma <ref>.

For any t, there exists a constant σ̅_1, such that

    x^t-x^t+1≥σ̅_1x^t-x(y^t+1, z^t).



Since this lemma is not related to the update of y, z, for notation simplicity, we denote K(x)=K(x, z^t; y^t+1).
The proof consists of two parts:


  * When the primal variables are updated via the block coordinate gradient descent scheme, it is a type of approximate gradient projection algorithm, namely,

    x^t+1=[x^t-c∇ K(x^t)+Θ(t)]_+,

where Θ(t) satisfies

    Θ(t)< ηx^t-x^t+1

for some positive constant η>0.

  * Prove the approximate gradient projection algorithm has the primal error bound.

We consider the first part. Notice that

    x^t+1_j   =   [x^t_j-c∇_x_j K(x^t(j))]_+^j,
       =   [x^t_j-c∇_x_j K(x^t)+c(∇_x_j K(x^t(j))-∇_x_j K(x^t))]_+^j
       =   [x^t-c∇_x_j K(x^t)+Θ_j(t)]_+^j,,

where Θ_j(t)=c(∇_x_j K(x^t(j))-∇_x_j K(x^t)).
Due to the Lipschitz continuity of the partial gradient of K, we have

    Θ_j(t)   ≤   c(L+p+Γσ̅^2)x^t(j)-x^t
       ≤   c(L+p+Γσ̅^2)√(∑_i=1^j-1x_i^t+1-x_i^t^2)
       ≤    c(L+p+Γσ̅^2)∑_i=1^j-1x_i^t-x_i^t+1,

where the last inequality is proved by just squaring both sides of the inequality.
Since ∑_i=1^j-1x_i^t-x_i^t+1≤∑_i=1^kx_i^t-x_i^t+1, we have

    Θ(t)   =   (Θ_1(t), ⋯, Θ_k(t))
       ≤   c(L+p+Γσ̅^2)k∑_i=1^kx^t_i-x^t+1_i
       ≤   c(L+p+Γσ̅^2)k^3/2x^t-x^t+1,

where the last inequality is due to Cauchy-Schwartz inequality.
This finishes the proof of (<ref>) with η=c(L+p+Γσ^2)k^3/2.
For the second part, we have

    x^t-   x   ^t+1=x^t-[x^t-c∇ K(x^t)+Θ(t)]_+
       ≥   x^t-[x^t-c∇ K(x^t)]_+-[x^t-c∇ K(x^t)]_+-[x^t-c∇ K(x^t)+Θ(t)]_+
       ≥   σ_1x^t-x(y^t+1, z^t)-[x^t-c∇ K(x^t)]_+-[x^t-c∇ K(x^t)+Θ(t)]_+
       ≥   σ_1x^t-x(y^t+1, z^t)-Θ(t)
       ≥   σ_1x^t-x(y^t+1, z^t)-ηx^t-x^t+1,

where (i) is because of the triangular inequality, (ii) is due to the error bound (<ref>) in Lemma <ref>, (iii) is due to the nonexpansiveness of the projection operator and (iv) is because of (<ref>). Hence setting σ̅_1=σ_1/(1+η) completes the proof.

Next we establish a simple lemma to ensure that the sufficient decrease result Lemma <ref> holds true for the multi-block case.

For any t, we have

    K(x^t, z^t; y^t+1)-K(x^t+1, z^t; y^t+1)≥1/2cx^t-x^t+1^2.



Notice that, by Assumption <ref>(c), the partial gradient of K with respect to any block is c^-1-Lipschitz continuous, so we have

    K(x^t(j), z^t; y^t+1)-K(x^t(j+1), z^t ;y^t+1)≥1/2cx^t_j-x^t+1_j^2,1≤ j≤ k.

Here x^t(0)=x^t.
Summing this from 0 to k-1 and using the fact that

    ∑_j=1^kx^t_j-x^t+1_j^2=x^t-x^t+1^2

yields the desired result.

This shows that the descent condition (<ref>) holds true for the multi-block case, which further implies Lemma <ref> remains valid.
Equipped with these, we conclude that the Algorithm <ref> converges globally.

Suppose Assumption <ref> holds and the parameters c and p satisfy

    1/L+p+Γσ̅^2>c>0,   p>-γ

and that the primal and dual stepsizes β and α are sufficiently small. Then the dual iterates {y^t} are bounded. Moreover, there holds

    lim_t→∞x^t+1-x^t   =   0,
    lim_t→∞(x^t, X^*)   =   0,
    lim_t →∞(z^t, X^*)   =   0,

and every limit point of the sequence {(x^t, y^t)} generated by Algorithm <ref> is a primal-dual stationary point of (<ref>).



§ NUMERICAL RESULTS

In this section, we show some numerical results to show that our algorithm is more efficient than some existing algorithms and the exponential weighting step is helpful for convergence.
We consider the quadratic programming problem:

    (QP)  min_x    1/2x^TQx+r^Tx
    s.t.    Ax=b,
       ℓ_i≤x_i ≤ u_i,

where x∈ℝ^n, Q∈ℝ^n× n, r∈ℝ^n, A∈ℝ^m× n, and b∈ℝ^m.


 §.§ Why do we need the auxiliary sequence {z^t}?

A natural question is whether we can set β=1 in Algorithm <ref> and thus eliminating the sequence {z^t} from the iterations. In this subsection, we give some numerical results showing that β<1 is needed for convergence
and hence the sequence {z^t} is necessary in Algorithm <ref>.
Here we set n=500, m=100. Q=-U^TU, where U is a 500× 500 matrix with each entry following the normal distribution N(0, 1). Every entry of A and r is generated from the normal distribution N(0, 1). b=Ax_0, with x_0∈ℝ^500 uniformly distributed over [0,5]^500.Moreover,  ℓ_i=0 and u_i=1000 for all 1≤ i≤ 500.

In our experiment, we set Γ=1000 and p=5000. Consider the following two cases.

Case 1. We use α=1000, 50, 1 and plot the curves for Ax^t-b in Fig. <ref>.

We see that Algorithm <ref> oscillate after 2× 10^6 iterations for these choices of α.

Case 2. We choose α=50, β=1, 0.02, 0.01 and plot f(x^t), Ax^t-b, x^t+1-z^t and 1/cx^t-x^t+1. This is shown in Fig. <ref>, Fig. <ref> and Fig. <ref>.

We see that Algorithm <ref> converges with β=0.02 and β=0.01. The algorithm with β=0.02 is faster, which suggests that we can try larger β to achieve higher convergence speed.
  
  










 §.§ Comparison with the penalty method


We compare our proximal-ADMM algorithm with the QP-AIPP algorithm of <cit.>  for solving the above QP with n=20,m=5, where Q,A,r,b are randomly generated from uniform distribution 𝒰[0,1], and ℓ_i and u_i are set to be 0 and 1 respectively.

Parameters for our algorithm are set as follows:
Γ=10, α=Γ/4, p=2L_Q+2Γσ^2, c=1/(2(L_Q+Γσ^2)), and β=0.5, where L_Q is the eigenvalue of Q with largest absolute value.

Parameter setting of the QP-AIPP algorithm are as follows (notations in Section 3.2 of <cit.>, AIPP Algorithm): σ=0.1, Γ̅=10^-3, ϵ̅=10^-3, c=10^4, L=cA_F^2+L_Q (L_Q is defined above), and λ=1/(2L_Q).

We use the sum of first order optimality gap and feasibility gap as a measure to evaluate the convergence behavior of the two algorithms, i.e., x^t-[x^t-∇_xL(x^t; λ^t)]_++Ax^t-b.
Here for our algorithm, we set λ^t=y^t while for the penalty method in <cit.>, λ^t=C(Ax^t-b).
We randomly generated two examples and the residuals are plotted in Fig. <ref>. 




 §.§ Comparison with the traditional ADMM

Lastly, we show that the algorithm with linearization steps is faster than the traditional ADMM where any sub-problem is solved to a high accuracy by a first-order algorithm.

We consider the 2-block case:

    [                     min x_1^TQ_1x_1+x_2^TQ_2x_2;                                A_1x_1+A_2x_2=b,;                             x_1∈ P_1, x_2∈ P_2. ]

Here Q_1 and Q_2 are (n/2)× (n/2) symmetric matrices, where  each entry
 is generated from a uniform distribution over [0,1]. A_1, A_2∈ℝ^m × (n/2) and b∈ℝ^m with each entry generated from the uniform distribution  over [0,1]. P_1=P_2=[0, 10]^n/2 are two boxes.
The following shows the convergence of the original ADMM and our smoothed proximal ADMM using the same measure as the above subsections.

This experiment shows that, as a single-loop algorithm, the proposed Proximal ADMM algorithm is more efficient than the original double-loop ADMM algorithm in which the inner-loop can be time consuming.


































































































































































































































































































§ PROOF OF COROLLARY <REF>


We prove by contradiction.
Suppose that there exist sequences {y^i}⊆ỹ^0+range(A),  ⊆ (ỹ^0+Range(A)) and {z^i} such that

    lim_i→∞x^i-(x^i)^+=lim_i→∞z^i-(x^i)^+=lim_i→∞Ax((y^i)^+, z^i)-b=0

while max{limsup_i→∞(x^i,X^*), limsup_i→∞(z^i, X^*)}>0.
Passing to a sub-sequence if necessary, we assume that max{lim_i→∞(x^i ,X^*), lim_i→∞(z^i, X^*), lim_i→∞x((y^i)^+, z^i)-x^*(z^i)}>0.
On one hand, since {y^i}⊆ỹ^0+Range(A) and Ax(y^i, z^i)-b→ 0,  {y^i} is bounded according to  (cf. Lemma <ref>). By further passing to a subsequence if necessary, we can assume that

    x^i→x̅,   z^i→z̅,   (y^i)^+→y̅.

Then we have x̅-x̅^+=0,  z̅-z̅^+=0.

Therefore x̅=x(y̅, z̅).
It follows from Lemma <ref> that
Ax(y̅, z̅)-b=0, which further implies

    (x̅, y̅, z̅)=(x̅^+, y̅^+, z̅^+).

In view of Lemma <ref>, we obtain that x̅=z̅∈ X^*, z̅∈ X^* and (x̅, y̅)∈ W.
Hence, we have (x̅, X^*)=(z̅, X^*)=0
On the other hand, we have max{(x̅, X^*), (z̅, X^*)}>0. A contradiction!








§ PROOF OF COROLLARY <REF>


Again we prove by contradiction.
Suppose the contrary so that there exists sequences {y^i} and {z^i} such that Ax((y^i)^+, z^i)-b→ 0 and limsup_i→∞x((y^i)^+,z^i) -x^*(z^i)>0.
Passing to a sub-sequence if necessary, we can assume  without loss of generality that lim_i→∞x((y^i)^+, z^i)-x^*(z^i)>0.
By Lemma <ref>, because Ax((y^i)^+, z^i)-b→ 0, {(y^i)^+} is bounded. Hence, there exists at least one limit point of {(y^i)^+}. Therefore, further passing to a sub-sequence if necessary, we assume that (y^i)^+→y̅ and z^i→z̅∈ P.
Combining the above analysis and using Lemma <ref>, we have

    Ax(y̅, z̅)-b=0.

It follows from the KKT condition for (<ref>) at z^i that y̅ (the limit of (y^i)^+) is the optimal dual multiplier for the problem:

    [                f(x)+p/2x-z̅^2;          [5pt]    Ax=b, x∈ P, ]

implying x(y̅, z̅)=x^*(z̅). This is a contradiction.



§ PROOF OF INEQUALITIES (<REF>)-(<REF>) IN LEMMA <REF>


We first prove (<ref>).
By the definition of K(x,y;z) (cf. (<ref>)) and Assumption <ref>(b), ∇ K(x,z^t;y^t+1) is Lipschitz continuous in x

    ∇_x K(x^t,z^t;y^t+1)-∇_x K(x(y^t+1,z^t),z^t;y^t+1)≤ (p+L+Γσ^2)x^t-x(y^t+1,z^t).


So the Lipschitz constant for c∇_xK is c(L+p+Γσ^2), where σ is the spectral norm of A.
Since γ_K=p+γ>0, it follows that cK(x,z^t;y^t+1) is strongly convex in x with modulus c(p+γ), so from <cit.> that the following global error bound holds

    x-[x-c∇ K(x,z^t;y^t+1)]_+≥σ_1x-x(y^t+1,z^t),   ∀ x∈ℝ^n,

where σ_1=cγ_K.

Specializing (<ref>) at x=x^t and noticing that

    [x^t-c∇ K(x^t,z^t;y^t+1)]_+=x^t+1,

we obtain (<ref>).

By the triangular inequality, we further obtain

    x^t+1-x(y^t+1,z^t)   ≤   x^t-x^t+1+x^t-x(y^t+1, z^t)
       ≤   x^t-x^t+1+1/σ_1x^t-x^t+1
       ≤   (1+1/σ_1)x^t-x^t+1.

This shows that (<ref>) holds with σ_2=σ_1/(1+σ_1).




plain
