A METHOD WITH CONVERGENCE RATES FOR OPTIMIZATION

PROBLEMS WITH VARIATIONAL INEQUALITY CONSTRAINTS∗

HARSHAL D. KAUSHIK† AND FARZAD YOUSEFIAN†

Abstract. We consider a class of optimization problems with Cartesian variational inequality

(CVI) constraints, where the objective function is convex and the CVI is associated with a monotone
mapping and a convex Cartesian product set. This mathematical formulation captures a wide range
of optimization problems including those complicated by the presence of equilibrium constraints,
arXiv:2007.15845v2 [[Link]] 14 Feb 2021

complementarity constraints, or an inner-level large scale optimization problem. In particular, an

important motivating application arises from the notion of efficiency estimation of equilibria in multi-
agent networks, e.g., communication networks and power systems. In the literature, the iteration
complexity of the existing solution methods for optimization problems with CVI constraints appears
to be unknown. To address this shortcoming, we develop a first-order method called averaging ran-
domized block iteratively regularized gradient (aRB-IRG). The main contributions include: (i) In
the case where the associated set of the CVI is bounded and the objective function is nondifferen-
tiable and convex, we derive new non-asymptotic suboptimality and infeasibility convergence rate
statements in a mean sense. We also obtain deterministic variants of the convergence rates when we
suppress the randomized block-coordinate scheme. Importantly, this paper appears to be the first
work to provide these rate guarantees for this class of problems. (ii) In the case where the CVI set
is unbounded and the objective function is smooth and strongly convex, utilizing the properties of
the Tikhonov trajectory, we establish the global convergence of aRB-IRG in an almost sure and a
mean sense. We provide the numerical experiments for computing the best Nash equilibrium in a
networked Cournot competition model.

Key words. first-order methods, variational inequalities, complexity analysis, efficiency of Nash
equilibria, iterative regularization, randomized block-coordinate

AMS subject classifications. 65K15, 49J40, 90C33, 91A10, 90C06

1. Introduction. Traditionally, the mathematical models and algorithms for

constrained optimization have been much focused on the cases where the functional
constraints are in the form of inequalities, equalities, or easy-to-project sets. However,
in a breadth of emerging applications in the control theory and economics, the system
constraints are too complex to be characterized in those forms. This may arise in
several network application domains where the optimization model is complicated by
the presence of equilibrium constraints, complementarity constraints, or an inner-level
large scale optimization problem. Accordingly, the goal in this paper lies in addressing
the following constrained optimization problem:

minimize f (x)
(PfVI )
subject to x ∈ SOL(X, F ).

Rn → R is a convex function and X ⊆ Rn is given as a Cartesian

Here, f : Q product,
d Pd
i.e., X , i=1 Xi , where Xi ⊆ Rni is convex for all i = 1, . . . , d and i=1 ni = n. We
consider F : X → Rn to be a monotone mapping. The term SOL(X, F ) denotes the
solution set of the variational inequality VI(X, F ) defined as follows: A vector x ∈ X
is said to be a solution to VI(X, F ) if for any y ∈ X, we have F (x)T (y − x) ≥ 0.
∗ Submitted to the editors DATE. A very preliminary version of this work appeared in the 2019

American Control Conference (ACC), IEEE, Philadelphia, PA, USA, 2019, pp. 3420–3425. https:
//[Link]/10.23919/ACC.2019.8815256.
Funding: Farzad Yousefian gratefully acknowledges the support of the NSF through CAREER
grant ECCS-1944500.
† School of Industrial Engineering & Management, Oklahoma State University, Stillwater, OK

74074, USA ([Link]@[Link], [Link]@[Link]).

1
2 H. D. KAUSHIK AND F. YOUSEFIAN

Variational inequalities, first introduced in late 1950s, are an immensely powerful

mathematical tool that can serve as a unifying framework for capturing a wide range of
applications arising in operations research, finance, and economics (cf. [13, 38, 11, 48]).
Importantly, as it will be discussed shortly, the block structure of the set X allows
for addressing Nash games as well as high-dimensional optimization problems. We
note that the problem (PfVI ) can represent a variety of the standard problems in
optimization and VI regimes. For example, when F (x) := 0n , the problem (PfVI ) is
equivalent to the canonical optimization problem minx∈X f (x). Also, when f (x) := 0,
the problem (PfVI ) is equivalent to solving VI(X, F ). More detailed examples that can
be reformulated as the problem (PfVI ) are presented below.
1.1. Motivating examples.
Example 1.1 (Efficiency estimation of equilibria). The main motivating applica-
tion arises from the notion of efficiency of equilibria in multi-agent networks, including
communication networks and power systems. In the noncooperative regimes, the sys-
tem behavior is governed by a collection of decisions (i.e., equilibrium) made by a set
of independent and self-interested agents. As a result of this noncooperative behavior
(i.e., game) among the agents, the global performance of the system may become
worse than the case where the agents cooperatively seek an optimal decision. A well-
known example is the Prisoners’ Dilemma where the costs of the players incurred by
the Nash equilibrium are superior to their costs when they cooperate [35]. Indeed,
it has been well-received in economics and computer science communities that Nash
equilibria of a game may not attain full efficiency. This perception has led to a surge
of research for understanding the quality of an equilibrium in noncooperative games.
In particular, addressing this question becomes imperative for network design [4] in
the areas of routing [10] and load balancing [39]. In such networks, a protocol de-
signer seeks the best equilibrium with respect to a global performance measure, i.e.,
the function f in (PfVI ). To this end, the notion of price of stability is defined as
the ratio of the best objective function value over the set of equilibria to the best
objective function value under no competition [34]. In regard to the choice of the
objective function f in (PfVI ), different approaches have been considered, including
the utilitarian function and the egalitarian function [34]. In the utilitarian approach,
f is defined as the summation of the individual objective functions of the agents,
while in the egalitarian approach, the maximum of the individual cost functions is
considered. In the context of network resource allocation where a monetary value is
measured, the utilitarian approach is also referred to as Marshallian aggregate surplus
(e.g., see [19]). In the following, we describe the details for the problem of selecting
the best equilibrium in Nash games, where we employ the utilitarian approach. Con-
sider a canonical Nash game among d players where the ith
player is associated with a
strategy x(i) ∈ Xi ⊆ Rni and a cost function gi x(i) ; x(−i) , where x(−i) denotes the
collection of actions of other players. Nash games arise in a wide range of problems
including communication networks [1, 2, 50], cognitive radio networks [47, 43, 27], and
power markets [22, 23, 45]. The game is defined formally as the following collection
of problems for all i = 1, . . . , d:
 
minimizex(i) gi x(i) ; x(−i)


Pi x(−i)
subject to x(i) ∈ Xi .
 
A Nash equilibrium (NE) is a tuple of strategies x∗ , x∗ (1) ; x∗ (2) ; . . . ; x∗ (d) where
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 3

no player can obtain a lower cost by deviating from his own strategy, given that the
strategies of the other players remain unchanged. It is known that (cf. Proposition
1.4.2 [13]) when for all i, Xi is a closed convex set and gi is a differentiable convex
(i)
function with respect
to x , the resulting equilibrium conditions of the Nash game
(−i)
given by (Pi x ), are compactly captured by a Cartesian VI(X, F ) where X ,
Qd (i) (−i)


i=1 Xi and F (x) , (F1 (x); . . . ; Fd (x)) with Fi (x) , ∇x(i) gi x ; x . The set


SOL(X, F ) will then represent the set of Nash equilibria to the game (Pi x(−i) ).
The best NE problem employing the utilitarian approach is formulated as follows:
(1.1)
Xd  
minimize gi x(i) ; x(−i)
i=1
d
!
Y     
(1) (−1) (d) (−d)
subject to x ∈ SOL Xi , ∇x(1) g1 x ;x ; . . . ; ∇x(d) gd x ;x .
i=1

In section 5, we solve the model (1.1) for a class of networked Nash-Cournot games.
Example 1.2 (High-dimensional constrained convex optimization). Another class
of problems that can be captured by the model (PfVI ) is as follows:

minimize f (x)
(1.2) subject to hj (x) ≤ 0 for all j = 1, . . . , J
Yd
Ax = b, x∈X, Xi ,
i=1

where x ∈ Rn , A ∈ Rm×n , b ∈ Rm , hj : Rn → R for all j, and n is possibly very large.

In the following, we show a case where the problem (1.2) can be cast as (PfVI ):
Lemma 1.3. Let the problem (1.2) be feasible and hj (x) be a continuously differ-
entiable convex function for all j. Let the set Xi ∈ Rni be nonempty, closed, and
convex for all i. Then, the problem (1.2) is equivalent to (PfVI ) where F : Rn → Rn
PJ
is defined as F (x) , AT (Ax − b) + j=1 max{0, hj (x)}∇hj (x).
Proof. See Appendix A.1.
Example 1.4 (Optimization problems with complementarity constraints). An-
other class of problems that can be addressed in this work is as follows:
minimize f (x)
(1.3)
subject to xT F (x) = 0, x ≥ 0, F (x) ≥ 0,

where F : Rn → Rn is a mapping. Then, problem (1.3) can be cast as (PfVI ) where

the set X is the nonnegative orthant, i.e., X , Rn+ (see Proposition 1.1.3 in [13]).
Example 1.5 (Optimization problems with nonlinear equality constraints). Con-
sider the following optimization problem:
minimize f (x)
(1.4)
subject to F (x) = 0.

where F : Rn → Rn is a mapping. Defining X , Rn , SOL(X, F ) is equal to the

feasible solution set of the problem (1.4). This implies that the problem (1.4) can be
captured by the formulation (PfVI ).
4 H. D. KAUSHIK AND F. YOUSEFIAN

1.2. Existing methods and the research gap. We first begin by providing
a brief overview of the solution methods for addressing a VI problem. Starting from
the seminal work of Lemke and Howson [29] and Scarf [41], who developed the first
solution methods for computing equilibria, in the past few decades, there has been a
surge of research on the development and analysis of the computational methods for
solving VIs. Perhaps this interest lies in the strong interplay between the VIs and the
formulation of optimization and equilibrium problems arising in many communication
and networking problems [42]. Korpelevich’s celebrated extragradient method [26]
and its extensions [32, 20, 7, 17, 8, 52, 15, 9, 54] were developed which require weaker
assumptions than their gradient counterparts. In the past decade, there has been a
trending interest in addressing VIs in the stochastic regimes. Among these, Jiang
and Xu [18] developed the stochastic approximation methods for solving VIs with
strongly monotone and smooth mappings. This work was later extended to the case
with merely monotone mappings [27, 21, 16] and nonsmooth mappings [53].

Algorithm 1.1 The existing SR scheme for solving problem (PfVI ) when f := 21 k · k2
1: Input: Set X, mapping F , and an initial regularization parameter η0 > 0.
2: for t = 0, 1, . . . do
3: Compute x∗ηt defined as x∗ηt ∈ SOL (X, F + ηt In ).
4: Update ηt to ηt+1 such that ηt+1 < ηt .
5: end for

Despite much advances in the theory and algorithms for VIs, solving the problem
(PfVI ) has remained challenging. To the best of our knowledge, the computational
complexity of the existing solution methods for addressing (PfVI ) is unknown. In
addressing the standard constrained optimization problems, Lagrangian duality and
relaxation rules have often proven to be very successful [6]. However, when it comes
to solving (PfVI ), the duality theory cannot be practically employed. This is primarily
because unlike in the standard constrained optimization problems where the objective
function provides a metric for distinguishing solutions, there is no immediate analog
in the VI problems. Inspired by the contributions of Andrey Tikhonov in 1980s on ad-
dressing illposed optimization problems, the existing methods for solving (PfVI ) share
in common a sequential regularization (SR) scheme presented by Algorithm 1.1. The
SR scheme is a two-loop framework where at each iteration, given a fixed parameter
ηt , a regularized VI denoted by VI (X, F + ηt In ) is required to be solved. In the spe-
cial case where f (x) := 12 kxk2 , it can be shown when ηt → 0, under the monotonicity
of the mapping F and closedness and convexity of the set X, any limit point of the
Tikhonov trajectory denoted by {x∗ηt }, where x∗ηt ∈ SOL (X, F + ηt In ), converges to
the least `2 –norm vector in SOL(X, F ) (cf. Chapter 12 in [13]). The SR approach is
associated with two main drawbacks: (i) It is a computationally inefficient scheme,
as it requires solving a series of increasingly more difficult VI problems. (ii) The
iteration complexity of the SR scheme in addressing the problem (PfVI ) is unknown.
Accordingly, the main goal in this work lies in the development of an efficient scheme
equipped with computational complexity analysis for solving the problem (PfVI ).
1.3. Summary of contributions. Our main contributions are as follows:
(i) Development of a single timescale method equipped with convergence rate guaran-
tees: In addressing (PfVI ), we develop an efficient first-order method called averaging
randomized block iteratively regularized gradient (aRB-IRG). The proposed method
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 5

is single timescale in the sense that, unlike the SR approach, it does not require solv-
ing a VI at each iteration. Instead, it only uses evaluations of the mapping F and
the subgradient of the objective function f at each iteration. In the first part of the
paper, we consider the case where the set X is bounded. We let f be a subdiffer-
entiable merely convex function and F be a monotone mapping. In Theorem 3.3,
we derive a suboptimality convergence rate in terms of the expected value of the
objective function. We also derive a convergence rate for the infeasibility that is char-
acterized by the expected value of a dual gap function. We also derive deterministic
variants of the aforementioned convergence rates when we suppress the randomized
block-coordinate scheme. In the second part of the paper, we consider the case where
the set X is unbounded and f is smooth and strongly convex. Utilizing the properties
of the Tikhonov trajectory, we establish the global convergence of the scheme in an
almost sure and a mean sense. To the best of our knowledge, this work appears to
be the first paper that provides the two rate statements for problems of the form
(PfVI ). In particular, the complexity analysis in this work contributes to the exist-
ing convergence theory in several previous papers including [49, 27, 21, 53, 24, 55].
Moreover, in the special case where the VI constraints represent the optimal solu-
tion set of an optimization problem, (PfVI ) captures a class of bilevel optimization
problems. This class of problems has been studied in a number of recent papers in
deterministic [46, 5, 40, 14], stochastic [3], and distributed regimes [51]. However, the
complexity analysis in the aforementioned papers lacks a suboptimality rate, or lacks
an infeasibility rate, or requires much stronger assumptions such as strong convexity
and smoothness of f .
(ii) Advancing the convergence rate properties of the randomized block-coordinate
schemes: Block-coordinate schemes, and specifically their randomized variants, have
been widely studied in addressing the standard optimization problems (e.g., see
[33, 37, 44, 12, 54]). However, in addressing VI problems, there are only a hand-
ful of recent papers, including [28, 54], that employ this technique and are equipped
with rate guarantees. The aforementioned papers address standard VI problems that
can be viewed a special case of the model (PfVI ) where f (x) := 0. In this work, we
extend the convergence and rate analysis of the randomized block-coordinate schemes
to the much broader regime of optimization problems with CVI constraints.
Outline of the paper: The paper is organized as follows. The proposed algo-
rithm is presented in section 2. The complexity analysis is provided in section 3. In
section 4, we provide the convergence analysis when the set X is unbounded. The
experimental results are in section 5, and the conclusions follow in section 6.
Notation and preliminary definitions: Throughout, a vector x ∈ Rn is as-
sumed to be a column vector and xT denotes the transpose of x. We
use x(i) ni
P∈d R to
th (1) (d)
denote the i block-coordinate of vector x where x = x ; . . . ; x and i=1 ni =
√
n. The Euclidean norm of a vector x is denoted by kxk, i.e., kxk , xT x. For a
mapping F : Rn → Rn , we denote the ith block-coordinate of F by Fi : Rn → Rni ,
i.e., F (x) = (F1 (x); . . . ; Fd (x)). A mapping F : Rn → Rn is said to be monotone
on a convex set X ⊆ Rn if for any x, y ∈ X, we have (F (x) − F (y))T (x − y) ≥ 0.
The mapping F is said to be µ–strongly monotone on a convex set X ⊆ Rn if µ > 0
and for any x, y ∈ X, we have (F (x) − F (y))T (x − y) ≥ µkx − yk2 . Also, F is
said to be Lipschitz with parameter L > 0 on the set X if for any x, y ∈ X, we have
kF (x)−F (y)k ≤ Lkx−yk. A continuously differentiable function f : Rn → R is called
µ–strongly convex on a convex set X if f (x) ≥ f (y) + ∇f (y)T (x − y) + µ2 kx − yk2 .
Function f is µ–strongly convex if and only if ∇f is µ–strongly monotone on X. For
6 H. D. KAUSHIK AND F. YOUSEFIAN

a convex function f with the domain dom(f ), the subgradient of f at x ∈ dom(f ) is

denoted by ∇f˜ (x) and it satisfies f (x) + ∇f
˜ (x)T (y − x) ≤ f (y) for all y ∈ dom(f ).
The subdifferential set of f at x is the set of all subgradients of f at x and is denoted
by ∂f (x). The Euclidean projection of vector s onto a convex set S is denoted by
PS (s), where PS (s) , argminy∈S ks − yk. We use In to denote the identity matrix
of size n × n. The probability of an event Z is denoted by Prob(Z) and the expec-
tation of a random variable z is denoted by E[z]. We use Rn+ and Rn++ to denote
{x ∈ Rn | x ≥ 0} and {x ∈ Rn | x > 0}, respectively.
2. Outline of the algorithm. In this section, we state the main assumptions
and present the proposed scheme for solving the optimization problem (PfVI ).
Assumption 2.1. Consider the problem (PfVI ) under the following conditions:
(a) The set Xi is nonempty, closed, and convex for all i = 1, . . . , d.
(b) The function f is convex and has bounded subgradients over the set X.
(c) The mapping F : Rn → Rn is continuous, monotone, and bounded over the set X.
(d) The optimal solution set of problem (PfVI ) in nonempty.
Assumption 2.1(b) implies that f is Lipschitz continuous over the set X. Under this
assumption, we address a broad class of problems of the form (PfVI ) where the objective
function is possibly nondifferentiable and nonstrongly convex. In the following, we
discuss the conditions under which Assumption 2.1(d) is satisfied.
Remark 2.2 (Existence of an optimal solution). Suppose Assumption 2.1(a), (b),
and (c) hold. The existence of an optimal solution to the problem (PfVI ) can be es-
tablished under different conditions. We provide two instances as follows: (i) Suppose
there exists a vector x̄ ∈ X such that the set X̄ , {x ∈ X : F (x)T (x − x̄) ≤ 0} is
bounded. Then, from Proposition 2.2.3 in [13], SOL(X, F ) is nonempty and compact.
Consequently, the Weierstrass’ Theorem implies the existence of an optimal solution
to the problem (PfVI ). (ii) Suppose the set X is compact. Then, from Corollary 2.2.5
in [13], the set SOL(X, F ) is nonempty and compact. Again, Assumption 2.1(d) is
guaranteed by the Weierstrass’ Theorem.
Throughout, we let CF > 0 denote the bound on the Euclidean norm of the mapping
F , i.e., kF (x)k ≤ CF for all x ∈ X. Also, we let Cf > 0 denote the bound on the
˜ (x)k ≤ Cf for all ∇f
norm of the subgradients of f , i.e., k∇f ˜ (x) ∈ ∂f (x) and x ∈ X.
The outline of the proposed method is presented by Algorithm 2.1. At iteration k, a
block-coordinate index ik is selected randomly as follows:
Assumption 2.3 (Block-coordinate selection rule). At each iteration k ≥ 0,
the random variable ik is generated from an independent and identically distributed
discrete probability distribution such that Prob (ik = i) = pi where pi > 0 for i ∈
Pd
{1, . . . , d} and i=1 pi = 1.
Then, the ithk block-coordinate of the iterate xk is updated using (2.1). Here, γk
denotes the stepsize at iteration k and ηk denotes the regularization parameter at it-
eration k. We note that these sequences are updated iteratively. Here, we incorporate
the information of the mapping F and the subgradient mapping ∇f ˜ by employing an
iterative regularization scheme. At each iteration, a projection operation is performed
onto a randomly selected set Xik . For this reason, the proposed algorithm finds some
relevance with the prior study on randomized projection methods such as [31]. We
will show that the convergence and rate analysis of the proposed method mainly rely
on the choices of {γk } and {ηk }. Accordingly, one key research objective in this sec-
tion will center around the development of suitable update rules for the two sequences
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 7

Algorithm 2.1 aRB-IRG

1: Input: A random initial point x0 ∈ X, x̄0 := x0 , initial stepsize γ0 > 0, initial
regularization parameter η0 > 0, a scalar 0 ≤ r < 1, and S0 := γ0r .
2: for k = 0, 1, . . . do
3: Generate a realization of random variable ik according to Assumption 2.3.
4: Evaluate Fik (xk ) and ∇ ˜ i f (xk ) where ∇f
˜ (xk ) ∈ ∂f (xk ).
k
5: Update xk as follows:
(   
(i ) ˜ i f (xk )
(i) PXik xk k − γk Fik (xk ) + ηk ∇ k
if i = ik ,
(2.1) xk+1 := (i)
xk if i 6= ik .

6: Obtain γk+1 and ηk+1 (cf. Theorem 3.3 and Theorem 4.11 for the update rules).
7: Update the averaged iterate x̄k as follows:
r
r Sk x̄k + γk+1 xk+1
(2.2) Sk+1 := Sk + γk+1 , x̄k+1 := .
Sk+1

8: end for

so that we can establish the convergence and derive rate statements. To obtain the
rate results, we employ an averaging step using the equations given by (2.2), where
the sequence {x̄k } is obtained as a weighted average of {xk }. The averaging weights
are characterized by the stepsize γk and a scalar r ∈ R. It will be shown that the rate
results can be provided when 0 ≤ r < 1 (cf. Theorem 3.3).
Remark 2.4. Importantly, unlike Algorithm 1.1, Algorithm 2.1 is a single
timescale scheme that does not require solving any inner-level VI problem. In par-
ticular, the update rule given by (2.1) mainly requires evaluations of random blocks
of the mappings F and ∇f ˜ . For this reason, (2.1) is computationally more efficient
than the step 3 in Algorithm 1.1.
2.1. Preliminaries. In the following, we provide some definitions and prelimi-
nary results that will be used to analyze the convergence of Algorithm 2.1.
Definition 2.5 (Distance function). For any x, y ∈ Rn , function D(x, y) is
Pd 2
defined as D(x, y) , i=1 p−1
i x(i) − y (i) , where pi is given by Assumption 2.3.
Remark 2.6. Under Assumption 2.3, we can relate the distance function D with
the `2 –norm as follows: pmin D(x, y) ≤ kx − yk2 ≤ pmax D(x, y) for all x, y ∈ Rn ,
where pmin , min1≤i≤d {pi } and pmax , max1≤i≤d {pi }.
One of the main challenges in the convergence analysis of computational methods
for solving VI problems lies in the lack of access to a standard metric to quantify the
quality of the solution iterates. This is in contrast with solving the standard optimiza-
tion problems where the objective function can serve as an immediate performance
metric for the underlying algorithm. Addressing this challenge in the literature of VI
problems has led to the study of so-called gap functions (cf. [13, 53]). Of these, in the
analysis of this section, we use the dual gap function defined as follows:
Definition 2.7 (The dual gap function [30]). Let a nonempty closed set X ⊆ Rn
and a mapping F : X → Rn be given. Then, for any x ∈ X, the dual gap function
GAP : X → R ∪ {+∞} is defined as GAP(x) , supy∈X F (y)T (x − y).
8 H. D. KAUSHIK AND F. YOUSEFIAN

Remark 2.8. When X 6= ∅, Definition 2.7 implies that the dual gap function is
nonnegative over X. It is also known that when F is continuous and monotone and
the set X is closed and convex, GAP(x∗ ) = 0 if and only if x∗ ∈ SOL(X, F ) (cf. [20]).
Thus, we conclude that under Assumption 2.1, the dual gap function is well-defined.
Definition 2.9 (Regularized mapping). Given a vector x ∈ X, a subgradient
˜ (x) ∈ ∂f (x), and an integer k ≥ 0, the regularized mapping Gk : X → R is defined
∇f
˜ (x). The ith block-coordinate of Gk is denoted by Gk,i .
as Gk (x) , F (x) + ηk ∇f
Definition 2.10 (History of the method). Throughout, we let the history of the
algorithm to be denoted by Fk , {x0 , i0 , i1 , . . . , ik−1 } for k ≥ 1, with F0 , {x0 }.
Next, we show that x̄k generated by Algorithm 2.1 is a well-defined weighted average.
Lemma 2.11 (Weighted averaging). Let {x̄k } be generated by Algorithm 2.1. Let
γr
us define the weights λk,N , PN k γ r for k ∈ {0, . . . , N } and N ≥ 0. Then, for any
j=0 j
PN
N ≥ 0, we have x̄N = k=0 λk,N xk . Also, when X is a convex set, we have x̄N ∈ X.
Proof. See Appendix A.2.
In the following, we define two terms that characterize the error between the true
maps with their randomized block variants.
Definition 2.12 (Randomized block error terms). Let Ui ∈ Rn×ni for i =
1, . . . , d be the collection of matrices such that In = [U1 , . . . , Ud ] ∈ Rn×n . Consider
the following definitions for k ≥ 0:

(2.3) ∆k , F (xk ) − p−1

ik Uik Fik (xk ),
˜ (xk ) − p−1 Ui ∇
δk , ∇f ik k
˜ i f (xk ).
k

Lemma 2.13 (Properties of ∆k and δk ). Consider Definition 2.12. We have:

 k | Fk ] = E[δk | Fk ] = 0.
2
(a) E[∆
(b) E k∆k k2 | Fk ≤ p−1 −1
 2

2
min − 1 CF and E kδk k | Fk ≤ pmin − 1 Cf .

Proof. See Appendix A.3

We will use the next result in deriving the suboptimality and infeasibility rate results.
Lemma 2.14 (Bounds on the harmonic series). Let 0 ≤ α < 1 be a given scalar.
1 +1)1−α PN (N +1)1−α
Then, for any integer N ≥ 2 1−α − 1, we have (N2(1−α) 1
≤ k=0 (k+1) α ≤ 1−α .
Proof. See Appendix A.4
3. Convergence rate analysis. In the following result, we derive an inequality
that will be later used to construct bounds on the objective function value and the
dual gap function at the averaged sequence generated by Algorithm 2.1.
Lemma 3.1. Consider the sequence {xk } in Algorithm 2.1. Suppose {γk } and
{ηk } are strictly positive sequences. Let Assumption 2.1 and Assumption 2.3 hold.
Let the auxiliary sequence {uk } ⊂ X be defined as uk+1 , PX (uk − γk (∆k + ηk δk )),
where u0 ∈ X is an arbitrary vector. Then, for all y ∈ X and all k ≥ 0, we have:
r−1
˜ (xk )T (xk − y) ≤ γk
γkr F (y)T (xk − y) + γkr ηk ∇f D (xk , y) + kuk − yk2


2
γkr−1
D (xk+1 , y) + kuk+1 − yk +γk (xk − uk )T (∆k + ηk δk )
2 r


−
2
2
+ γkr+1 k∆k k2 + ηk2 kδk k2 + 0.5pi−1 γk1+r k Gk,ik (xk )k .


(3.1) k
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 9

Proof. Let k ≥ 1 be fixed. From Definition 2.5 and (2.1), for any y ∈ X, we have:

2 d 2
(i ) (i)
X
(3.2) D (xk+1 , y) = p−1
ik
k
xk+1 − y (ik ) + p−1
i xk − y (i) .
i=1, i6=ik

2
(i )
Next, we find a bound on the term k
xk+1 − y (ik ) . From the block structure of
X, we have y (ik ) ∈ Xik . Invoking the nonexpansiveness property of the projection
mapping, the update rule (2.1), Definition 2.9, and the preceding relation, we obtain:
2 2
(i ) (i )
k
xk+1 − y (ik ) ≤ xk k − γk Gk,ik (xk ) − y (ik ) .

Combining the preceding relation with (3.2), we obtain:

d 2 2
(i) (i )
X
D (xk+1 , y) ≤ p−1
i xk − y (i) + p−1
ik xk k − y (ik )
i=1, i6=ik
 T
(ik ) 2
− 2 p−1
ik γ k x k − y (ik )
Gk,ik (xk ) + p−1 2
ik γk k Gk,ik (xk )k .

Invoking Definition 2.5 again, we obtain:

(3.3)
 T
(ik ) 2
D (xk+1 , y) ≤ D (xk , y) − 2 p−1
ik γ k x k − y (ik )
Gk,ik (xk ) + p−1 2
ik γk k Gk,ik (xk )k .

From Definition 2.12 and Definition 2.9, we can write:

 T
(i )
p−1
ik xk k − y (ik ) Gk,ik (xk ) = p−1 T
ik (xk − y) (Uik Gk,ik (xk ))
 
= p−1
ik (x k − y) T
U i k
F i k
(x k ) + η k U i k
˜ i f (xk )
∇ k
 
= (xk − y)T F (xk ) − ∆k + ηk ∇f ˜ (xk ) − ηk δk = (xk − y)T (Gk (xk ) − ∆k − ηk δk ) .

Combining the preceding inequality and relation (3.3), we obtain:

(3.4)
2
D (xk+1 , y) ≤ D (xk , y) −2γk (xk − y)T (Gk (xk ) − ∆k − ηk δk ) + p−1 2
ik γk k Gk,ik (xk )k .

Consider the definition of the auxiliary sequence {uk } in Lemma 3.1. Invoking the
nonexpansiveness property of the projection again, we can obtain:

kuk+1 − yk2 ≤ kuk − γk (∆k + ηk δk ) − yk2

= kuk − yk2 −2γk (uk − y)T (∆k + ηk δk ) + γk2 k∆k + ηk δk k2 .

Thus, we have:

kuk+1 − yk2 ≤ kuk − yk2 −2γk (uk − y)T (∆k + ηk δk ) + 2γk2 k∆k k2 + 2γk2 ηk2 kδk k2 .

Adding the preceding inequality and the inequality (3.4) together, we obtain:

2γk (xk − y)T Gk (xk ) ≤ D (xk , y) + kuk − yk2 − D (xk+1 , y) + kuk+1 − yk2



10 H. D. KAUSHIK AND F. YOUSEFIAN

(3.5)
2
+2γk (xk − uk )T (∆k + ηk δk ) + 2γk2 k∆k k2 + ηk2 kδk k2 + p−1 2


ik γk k Gk,ik (xk )k .

From the monotonicity property of the mapping F and Definition 2.9, we have:
˜ (xk )T (xk − y).
(xk − y)T Gk (xk ) ≥ (xk − y)T F (y) + ηk ∇f
This provides a lower bound on the left-hand side of (3.5). The inequality (3.1) is
γkr−1
obtained by substituting this bound in (3.5) and multiplying both sides by 2 .
In the following, we develop upper bounds for suboptimality and infeasibility of
the weighted average iterate generated by Algorithm 2.1. Both of these error bounds
are characterized in terms of the stepsize and the regularization parameter.
Proposition 3.2 (Error bounds for Algorithm 2.1). Let the sequence {x̄k } be
generated by Algorithm 2.1, where 0 ≤ r < 1. Suppose {γk } and {ηk } are strictly
positive and nonincreasing sequences. Let Assumption 2.1 and Assumption 2.3 hold
and assume that the set X is bounded, i.e., kxk ≤ M for all x ∈ X and some M > 0.
(a) Let x∗ be an optimal solution to problem (PfVI ). Then, for all N ≥ 1:
r−1
4M 2 γN PN  
−1 r+1 2 2 2
ηN + η
k=0 k γ k C F + η C
k f
(3.6) E[f (x̄N )] − f (x∗ ) ≤ PN .
pmin k=0 γk r

(b) Consider the dual gap function in Definition 2.7. Then, for all N ≥ 1:
PN  
r−1
4M 2 γN + k=0 γkr 2pmin ηk Cf M + γk CF2 + γk ηk2 Cf2
(3.7) E[GAP (x̄N )] ≤ PN .
pmin k=0 γkr
Proof. We define the following terms for all k ≥ 0, that appear in (3.1):
Θk,1 , γkr (xk − uk )T (∆k + ηk δk ), Θk,2 , γkr+1 k∆k k2 + ηk2 kδk k2 ,

2
(3.8) Θk,3 , 0.5p−1 1+r
ik γ k k Gk,ik (xk )k .
Next, we estimate the expected values of these terms. Consider the notation of Fk
given by Definition 2.10. Note that xk is Fk –measurable. Also, from the definition of
uk in Lemma 3.1, uk is Fk –measurable. Note, however, that Θk,j is Fk+1 –measurable
for all j ∈ {1, 2, 3}. Taking these into account and using the total probability law, for
any k ≥ 0 and j ∈ {1, 2, 3}, we have E[Θk,j ] = EFk [Eik [Θk,j | Fk ]]. From this relation
and Lemma 2.13, we have for any k ≥ 0:
E[Θk,2 ] = p−1

r+1 2
CF + ηk2 Cf2 .


(3.9) E[Θk,1 ] = 0, min − 1 γk

Also, using Definition 2.9 and the triangle inequality, we can write:
d
X  
2
Eik [Θk,3 | Fk ] = pi 0.5p−1 1+r
i γk k Gk,i (xk )k
i=1
d 
X 
≤ γk1+r ˜ i f (xk )k2 = γ 1+r kF (xk )k2 + ηk2 k∇f
kFi (xk )k2 + ηk2 k∇ ˜ (xk )k2 .
k
i=1

From the preceding inequality, we obtain:

E[Θk,3 ] ≤ γk1+r CF2 + ηk2 Cf2 .


(3.10)
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 11

We are now ready to show the inequalities (3.6) and (3.7) as follows:
(a) Consider (3.1). From the definition of subgradients of the convex function f , we
˜ (xk )T (xk − y). Thus, from (3.8) we obtain for any y ∈ X:
have that f (xk ) − f (y) ≤ ∇f

γkr−1
γkr F (y)T (xk − y) + γkr ηk (f (xk ) − f (y)) ≤ D (xk , y) + kuk − yk2


2
γkr−1
D (xk+1 , y) + kuk+1 − yk2 + Θk,1 + Θk,2 + Θk,3 .


−
2

Let us substitute y := x∗ , where x∗ denotes an optimal solution to (PfVI ). Note that

x∗ must be a feasible solution to (PfVI ), i.e., F (x∗ )T (xk − x∗ ) ≥ 0. Thus, we obtain:

γkr−1
γkr ηk (f (xk ) − f (x∗ )) ≤ D (xk , x∗ ) + kuk − x∗ k2


2
γkr−1
D (xk+1 , x∗ ) + kuk+1 − x∗ k2 + Θk,1 + Θk,2 + Θk,3 .


(3.11) −
2
r−1  
γk−1 2
Dividing both sides by ηk and adding and subtracting 2ηk−1 D (xk , x∗ ) + kuk − x∗ k
in the right-hand side of (3.11), we obtain for k ≥ 1:
r−1 
γk−1 2

γkr (f (xk ) − f (x∗ )) ≤ D (xk , x∗ ) + kuk − x∗ k
2ηk−1
γkr−1 
2

− D (xk+1 , x∗ ) + kuk+1 − x∗ k
2ηk
(3.12) !
1 γkr−1 γ r−1 
2

+ − k−1 D (xk , x∗ ) + kuk − x∗ k + ηk−1 (Θk,1 + Θk,2 + Θk,3 ) .
2 ηk ηk−1

γ r−1 γ r−1
Since r − 1 < 0 and that {γk } and {ηk } are nonincreasing, we have kηk − ηk−1
k−1
≥ 0.
∗
Also, from the boundedness of the set X, since xk , x , and uk belong to X, using
Remark 2.6 and the triangle inequality, we have:

(3.13)
2 ∗ 2 ∗ 2

8M 2
D (xk , x∗ ) + kuk − x∗ k ≤ p−1
min kx k − x k + kuk − x k ≤ 4M 2
p−1
min + 1 ≤ .
pmin

Summing over (3.12) from k = 1 to N and using (3.13), we obtain:

N
X γ0r−1  2

γkr (f (xk ) − f (x∗ )) ≤ D (x1 , x∗ ) + ku1 − x∗ k
2η0
k=1
r−1 N
γ r−1
  X
γN
+ 4M 2 p−1
min − 0 + ηk−1 (Θk,1 + Θk,2 + Θk,3 ) ,
ηN η0
k=1

where we drop the nonpositive term. From relation (3.11) when k = 0, we have:

γ0r−1
γ0r (f (x0 ) − f (x∗ )) ≤ D (x0 , x∗ ) + ku0 − x∗ k2


2η0
12 H. D. KAUSHIK AND F. YOUSEFIAN

γ0r−1
D (x1 , x∗ ) + ku1 − x∗ k2 + η0−1 (Θ0,1 + Θ0,2 + Θ0,3 ) .


−
2η0

Adding
PN the last two inequalities, multiplying and dividing the left-hand side by
r
k=0 k , and then, invoking Lemma 2.11 and convexity of f , we obtain:
γ

N
!
X γ0r−1  2

γkr (f (x̄N ) − f (x∗ )) ≤ D (x0 , x∗ ) + ku0 − x∗ k
2η0
k=0
r−1 N
γ r−1
 
γN X
+ 4M 2 p−1
min − 0 + ηk−1 (Θk,1 + Θk,2 + Θk,3 ) .
ηN η0
k=0

Taking the expectation on both sides and invoking (3.13), we obtain:

4M 2 p−1 r−1
min γN
PN
∗ ηN + k=0 ηk−1 E[Θk,1 + Θk,2 + Θk,3 ]
E[f (x̄N )] − f (x ) ≤ PN .
k=0 γkr

From the relations (3.9) and (3.10), we obtain:

4M 2 p−1 r−1   
min γN
PN
ηN + k=0 ηk−1 p−1 r+1
min γk CF2 + ηk2 Cf2
E[f (x̄N )] − f (x∗ ) ≤ PN ,
r
k=0 γk

which implies the inequality (3.6).

(b) From the Cauchy-Schwarz inequality, the definitions of Cf and M , and the triangle
˜ (xk )T (y − xk ) ≤ ∇f
inequality, we have ∇f ˜ (xk ) kxk − yk ≤ 2Cf M. Adding the
preceding inequality with the relation (3.1), from (3.8) we obtain:

γkr−1
γkr F (y)T (xk − y) ≤ D (xk , y) + kuk − yk2


2
γkr−1
D (xk+1 , y) + kuk+1 − yk2 + 2γkr ηk Cf M + Θk,1 + Θk,2 + Θk,3 .


(3.14) −
2
r−1  
γk−1 2
Adding and subtracting the term 2 D (xk , y) + kuk − yk , we obtain:

r−1 
γk−1 2

γkr F (y)T (xk − y) ≤ D (xk , y) + kuk − yk
2
r−1 
γ 2
 1
 2

− k D (xk+1 , y) + kuk+1 − yk + γkr−1 − γk−1
r−1
D (xk , y) + kuk − yk
2 2
+ 2γkr ηk Cf M + Θk,1 + Θk,2 + Θk,3 .

Substituting the bound given by (3.13) in the preceding relation, we obtain:

r−1 
γk−1 2

γkr F (y)T (xk − y) ≤ D (xk , y) + kuk − yk
2
γkr−1 
2

D (xk+1 , y) + kuk+1 − yk + 4M 2 p−1 r−1 r−1


− min γk − γk−1
2
+ 2γkr ηk Cf M + Θk,1 + Θk,2 + Θk,3 .
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 13

Summing both sides from k = 1 to N , we obtain:

N
X γ0r−1  2

γkr F (y)T (xk − y) ≤ D (x1 , y) + ku1 − yk + 4M 2 p−1 r−1
− γ0r−1


min γN
2
k=1
N
X
(3.15) + (2γkr ηk Cf M + Θk,1 + Θk,2 + Θk,3 ) .
k=1

Writing the inequality (3.14) for k = 0, we have:

γ0r−1
γ r−1
γ0r F (y)T (x0 − y) ≤ D (x0 , y) + ku0 − yk2 − 0 D (x1 , y) + ku1 − yk2


2 2
(3.16) + 2γ0r η0 Cf M + Θ0,1 + Θ0,2 + Θ0,3 .

Adding (3.15) and (3.16) together, we obtain:

N
X γ0r−1  2

γkr F (y)T (xk − y) ≤ D (x0 , y) + ku0 − yk + 4M 2 p−1 r−1
− γ0r−1


min γN
2
k=0
N
X
(3.17) + (2γkr ηk Cf M + Θk,1 + Θk,2 + Θk,3 ) .
k=0

PN
Recalling x̄N = k=0 λk,N xk in Lemma 2.11, applying the bound given by (3.13),
and using the triangle inequality, we obtain:
N
! N
X X
γkr F (y)T (x̄N − y) ≤4M 2 p−1 r−1
min γN + (2γkr ηk Cf M + Θk,1 + Θk,2 + Θk,3 ) .
k=0 k=0

Taking the supremum with respect to y over the set

PNX from the left-hand side, invoking
Definition 2.7, and then dividing both sides by k=0 γkr , we obtain:
PN
4M 2 p−1 r−1
min γN + k=0 (2γkr ηk Cf M + Θk,1 + Θk,2 + Θk,3 )
GAP (x̄N ) ≤ PN .
r
k=0 γk

Taking the expectation on both sides, using the relations (3.9) and (3.10), and rear-
ranging the terms, we obtain the inequality (3.7).
We are now ready to present the convergence rate results of the proposed method.
Theorem 3.3 (Convergence rate statements for Algorithm 2.1). Consider Al-
gorithm 2.1. Let Assumption 2.1 and Assumption 2.3 hold and assume that the set
X is bounded such that kxk ≤ M for all x ∈ X and some M > 0. Suppose for all
k ≥ 0, γk := √γk+1
0 η0
and ηk := (k+1) b , where γ0 > 0, η0 > 0, and 0 < b < 0.5. Then,

for any 0 ≤ r < 1, the following results hold:

2
(i) Let x∗ be an optimal solution to the problem (PfVI ). Then, for all N ≥ 2 1−r − 1:
  
2 2 2
2 − r  4M γ 0 2 CF + η C
0 f 1
(3.18) E[f (x̄N )] − f (x∗ ) ≤ +  .
pmin η0 γ0 0.5 − 0.5r + b (N + 1)0.5−b
14 H. D. KAUSHIK AND F. YOUSEFIAN

2
(ii) Consider the dual gap function in Definition 2.7. Then, for all N ≥ 2 1−r − 1:

(3.19)  
 
2 − r  4M 2 γ0 CF2 + η02 Cf2 2pmin Cf M η0  1
E[GAP (x̄N )] ≤ + + .
pmin γ0 0.5 − 0.5r 1 − 0.5r − b (N + 1)b

Proof. Let us define the following terms:

N r−1 N
X 4M 2 γN
X
ΛN,1 , pmin γkr , ΛN,2 , , ΛN,3 , CF2 + η02 Cf2 ηk−1 γkr+1 ,
ηN
k=0 k=0
N N
r−1

X X
ΛN,4 , 4M 2 γN , ΛN,5 , CF2 + η02 Cf 2
γkr+1 , ΛN,6 , 2pmin Cf M ηk γkr .
k=0 k=0

Note that from (3.6) and (3.7), we have:

ΛN,2 + ΛN,3 ΛN,4 + ΛN,5 + ΛN,6
(3.20) E[f (x̄N )] − f (x∗ ) ≤ , E[GAP (x̄N )] ≤ .
ΛN,1 ΛN,1
Next, we apply Lemma 2.14 to estimate the terms ΛN,i . Substituting γk and ηk by
their update rules, we obtain:
N
X γ0r pmin γ0r (N + 1)1−0.5r
ΛN,1 = pmin ≥ ,
(k + 1)0.5r 2(1 − 0.5r)
k=0
4M 2 (N + 1)0.5(1−r)+b 4M 2 (N + 1)0.5(1−r)
ΛN,2 = , Λ N,4 = ,
η0 γ01−r γ01−r
   
XN CF2 + η02 Cf2 γ01+r γ01+r CF2 + η02 Cf2 (N + 1)1−0.5(1+r)+b
ΛN,3 = ≤ ,
η (k + 1)0.5(1+r)−b
k=0 0
η0 (1 − 0.5(1 + r) + b)
 
N
X γ0r+1 C 2
F + η 2 2
0 Cf γ0r+1 (N + 1)1−0.5(1+r)
= CF2 + η02 Cf2


ΛN,5 ≤ ,
(k + 1)0.5(1+r)
k=0
1 − 0.5(1 + r)
N
X 1 2pmin Cf M η0 γ0r (N + 1)1−0.5r−b
ΛN,6 = 2pmin Cf M η0 γ0r 0.5r+b
≤ .
(k + 1) 1 − 0.5r − b
k=0

For these inequalities to hold, we need to ensure that the conditions of Lemma 2.14
are met. Accordingly, we must have 0 ≤ 0.5r < 1, 0 ≤ 0.5(1 + r) − b < 1, 0 ≤
0.5r + b < 1, and 0 ≤ 0.5(1 + r) < 1. These relations hold because 0 ≤ r < 1 and
0 < b < 0.5. Another set of conditions when applying Lemma 2.14 includes N ≥
max 21/(1−0.5r) , 21/(1−0.5(1+r)+b) , 21/(1−0.5r−b) , 21/(1−0.5(1+r)) − 1. This relation is
2
indeed satisfied as a consequence of N ≥ 2 1−r − 1, 0 < b < 0.5, and 0 ≤ r < 1. We
conclude that all the necessary conditions for applying Lemma 2.14 and obtaining the
aforementioned bounds for the terms ΛN,i are satisfied. To show that the inequalities
(3.18) and (3.19) hold, it suffices to substitute the preceding bounds on the terms
ΛN,i into the two inequalities given by (3.20). The details are as follows:

4M 2 (N + 1)0.5−0.5r+b

ΛN,2 + ΛN,3 2−r
E[f (x̄N )] − f (x∗ ) ≤ = r
ΛN,1 pmin γ0 (N + 1) 1−0.5r
η0 γ01−r
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 15
  
 C 2 + η 2 C 2 (N + 1)0.5−0.5r+b
γ01+r

F 0 f
+ .
η0 0.5 − 0.5r + b

The inequality (3.18) is obtained by rearranging the terms in the preceding relation.

4M 2 (N + 1)0.5−0.5r

ΛN,4 + ΛN,5 + ΛN,6 2−r
E[GAP (x̄N )] ≤ ≤ r
ΛN,1 pmin γ0 (N + 1) 1−0.5r
γ01−r
  
CF2 + η02 Cf2 γ0r+1 (N + 1)0.5−0.5r r
2pmin Cf M η0 γ0 (N + 1) 1−0.5r−b
+ + .
0.5 − 0.5r 1 − 0.5r − b

Then, (3.19) can be obtained by rearranging the terms in the preceding inequality.
Remark 3.4 (Iteration complexity of Algorithm 2.1). As an immediate result
from Theorem 3.3, choosing γk := √γk+1
0
and ηk := √
4
η0
k+1
, we obtain:
 
1
E[f (x̄N ) − f (x∗ )] = E[GAP (x̄N )] = O √
4
.
N

This implies that Algorithm 2.1 achieves an iteration complexity of O −4 in solving

(PfVI ), where  > 0 denotes the expected tolerance in both of the suboptimality and
infeasibility metrics.
The rate statements derived in Theorem 3.3 are in a mean sense. In the following, we
consider a deterministic variant of Algorithm 2.1 where we suppress the randomized
block-coordinate scheme. The outline of this deterministic method is presented by
Algorithm 3.1. In Corollary 3.5, we show that non-asymptotic deterministic rate
statements can be derived for Algorithm 3.1.

Algorithm 3.1 a-IRG

1: Input: An arbitrary initial point x0 ∈ X, x̄0 := x0 , initial stepsize γ0 > 0, initial
regularization parameter η0 > 0, a scalar 0 ≤ r < 1, and S0 := γ0r .
2: for k = 0, 1, . . . do
3: Evaluate F (xk ) and ∇f ˜ (xk ) where ∇f
˜ (xk ) ∈ ∂f (xk ).
4: For all i ∈ {1, . . . , d}, do the following updates:
  
(i) (i) ˜ i f (xk ) .
(3.21) xk+1 := PXi xk − γk Fi (xk ) + ηk ∇

5: Obtain γk+1 and ηk+1 (cf. Corollary 3.5 for the update rules).
6: Update the averaged iterate x̄k as follows:
r
r Sk x̄k + γk+1 xk+1
(3.22) Sk+1 := Sk + γk+1 , x̄k+1 := .
Sk+1

7: end for

Corollary 3.5 (Convergence rate statements for Algorithm 3.1). Consider

Algorithm 3.1. Let Assumption 2.1 hold and assume that the set X is bounded such
that kxk ≤ M for all x ∈ X and some M > 0. Suppose for k ≥ 0, γk := √γk+10
and
16 H. D. KAUSHIK AND F. YOUSEFIAN

η0
ηk := (k+1) b , where γ0 > 0, η0 > 0, and 0 < b < 0.5. Then, for any 0 ≤ r < 1, the

following results hold:

2
(i) Let x∗ be an optimal solution to the problem (PfVI ). Then, for all N ≥ 2 1−r − 1:
  
2 − r 4M 2 γ0 CF2 + η02 Cf2 1
(3.23) f (x̄N ) − f (x∗ ) ≤  +  .
η0 γ0 0.5 − 0.5r + b (N + 1)0.5−b

2
(ii) Consider the dual gap function in Definition 2.7. Then, for all N ≥ 2 1−r − 1:
   
4M 2 γ0 CF2 + η02 Cf2 2C M η 1
f 0 
(3.24) GAP (x̄N ) ≤ (2 − r)  + + .
γ0 0.5 − 0.5r 1 − 0.5r − b (N + 1)b

Proof. See Appendix A.6.

4. Addressing the case where X is unbounded. The convergence and rate
statements provided by Theorem 3.3 require the set X to be bounded. We, however,
note that in some applications, e.g., in the models presented in Example 1.4 and
Example 1.5, this assumption may not hold. Accordingly, in this section, our aim is
to analyze the convergence of Algorithm 2.1 when X is unbounded. To this end, we
consider the following main assumption:
Assumption 4.1. Consider problem (PfVI ) under the following conditions:
(a) The set Xi is nonempty, closed, and convex for all i = 1, . . . , d.
(b) The function f is continuously differentiable and µf –strongly convex over X.
(c) The mapping F : Rn → Rn is continuous and monotone over X.
(d) The solution set SOL(X, F ) is nonempty.
Remark 4.2 (Existence and uniqueness of the optimal solution). Under Assump-
tion 4.1, the constraint set of (PfVI ), i.e., SOL(X, F ), is nonempty, closed, and convex.
The convexity of this set is implied by Theorem 2.3.5 in [13] and its closedness prop-
erty is obtained by the continuity of the mapping F and closedness of the set X.
Because in the problem (PfVI ), the objective function f is strongly convex and that
the constraint set is nonempty, closed, and convex, we conclude from Proposition 1.1.2
in [6] that the problem (PfVI ) has a unique optimal solution. Throughout this section,
we let x∗ denote this unique optimal solution.
4.1. Preliminaries. In this part, we provide some preliminary results that will
be used in the convergence analysis. We begin by defining a generalized variant of
the Tikhonov trajectory that is associated with the problem of interest in this paper.
Definition 4.3 (Tikhonov trajectory). Consider the problem (PfVI ) under As-
sumption 4.1. Let {ηk } be a sequence of strictly positive scalars for all k ≥ 0, and
x∗ηk ∈ X denote the unique solution to the regularized variational inequality problem
given by VI (X, F + ηk ∇f ). Then, the sequence x∗ηk is defined as the Tikhonov


trajectory associated with the problem (PfVI ).

Remark 4.4. The uniqueness of the solution of VI (X, F + ηk ∇f ) in Definition 4.3
is due to the strong monotonicity of the mapping F + ηk ∇f and closedness and
convexity of the set X (see Theorem 2.3.3 in [13]). Definition 4.3 generalizes the
notion of Tikhonov trajectory provided in [13] in the following way: in [13], x∗ηk is
defined as the solution to the regularized problem VI (X, F + ηk In ). This is indeed
the special case where we choose f (x) := 21 kxk2 in Definition 4.3.
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 17

To analyze the convergence, we utilize the properties of the Tikhonov trajectory. The
following result ascertains the asymptotic convergence of this trajectory to the optimal
solution of the problem (PfVI ). It also provides an upper bound on the error between
any two successive vectors of the trajectory.
Lemma 4.5. Consider Definition 4.3 and let Assumption 4.1 hold. Let {ηk } be a
sequence such that limk→∞ ηk = 0 and ηk > 0 for all k ≥ 0. Then:
(a) The Tikhonov trajectory {x∗ηk } converges to a unique limit point, that is x∗ .
C̄f ηk−1
(b) There exists C̄f > 0 such that x∗ηk − x∗ηk−1 ≤ µf 1− ηk for all k ≥ 1.
Proof. See Appendix A.5.
The following lemmas will be employed to establish the asymptotic convergence result.
Lemma 4.6 (Theorem 6, page 75 in [25]). Let {ut } ⊂ Rn denote a sequence of
vectors where limt→∞ ut = û. Also, let {αk } denote a sequence of P
strictly positive
P∞ k
n αt ut
scalars such that k=0 αk = ∞. Suppose vk ∈ R is defined by vk , Pt=0 k for all
t=0 αt
k ≥ 0. Then, lim vk = û.
k→∞

Lemma 4.7 (Lemma 10, page 49 in [36]). Let {vk } be a sequence of nonnegative
random variables, where E[v0 ] < ∞, and let {αk } and {βk } be deterministic scalar
P|v
sequencesPsuch that E[vk+1
∞ ∞
0 , . . . , vk ] ≤ (1 − αk )vk + βk for all k ≥ 0, 0 ≤ αk ≤ 1,
βk ≥ 0, k=0 αk = ∞, k=0 βk < ∞, and limk→∞ αβkk = 0. Then, vk → 0 almost
surely and lim E[vk ] = 0.
k→∞

4.2. Convergence analysis. As a key step toward performing the convergence

analysis for Algorithm 2.1 when the set X is unbounded, next we derive a recursive
inequality for the distance between the generated sequence {xk } by the algorithm and
the Tikhonov trajectory {x∗ηk }. To this end, we first make the following assumption:
Assumption 4.8. Consider the problem (PfVI ) under the following assumptions:
(a) There exist nonnegative scalars LF and BF such that for all x, y ∈ X:

kF (x) − F (y)k2 ≤ L2F kx − yk2 + BF .

(b) The gradient mapping ∇f is Lipschitz with parameter Lf > 0.

Remark 4.9. By allowing LF or BF to be zero, Assumption 4.8 provides a unifying
structure for considering both smooth and nonsmooth cases. In particular, when
LF = 0, part (a) refers to a bounded, but possibly non-Lipschitzian mapping F . Also,
when BF = 0, part (a) refers to a Lipschitzian, but possibly unbounded mapping F .
The following recursive relation will play a key role in establishing the convergence.
Lemma 4.10 (A recursive error bound for Algorithm 2.1). Consider the se-
quence {xk } in Algorithm 2.1. Let Assumption 4.1, Assumption 2.3, and Assump-
tion 4.8 hold. Suppose {γk } and {ηk } are nonincreasing and strictly positive where
µf pmin
limk→∞ ηk = 0 and γηkk ≤ 2p for all k ≥ 0. Then, for all k ≥ 1:
max (LF +η0 Lf )
2 2 2

 pmax  pmin µf γk ηk   
E D xk+1 , x∗ηk |Fk ≤ D xk , x∗ηk−1


1−
pmin 2
2 2
C̄f (µf γ0 η0 + 2/pmin ) ηk−1

(4.1) + − 1 + 2γk2 BF .
µ3f pmin γk ηk ηk
18 H. D. KAUSHIK AND F. YOUSEFIAN

Proof. From Definition 2.5, we have:

2 d 2
(ik ) (ik ) X (i) (i)
D xk+1 , x∗ηk = p−1 − x∗ηk p−1 xk − x∗ηk


(4.2) ik xk+1 + i .
i=1, i6=ik

(i ) 2
(i )
k
Next, we find a bound on the term xk+1 − x∗ηk k . From the properties of the natu-
ral map (cf. Proposition 1.5.8 in [13]), Definition 2.9, and that x∗ηk ∈ X, we have x∗ηk =
PX x∗ηk − γk Gk x∗ηk . From Assumption 4.1(a) and that x∗ηk ∈ SOL (X, Gk ) ⊆ X,



(i )
we have x∗ηk k ∈ Xik . Invoking the nonexpansiveness property of the projection map-
ping, (2.1), and the preceding relation, we obtain:
(ik ) 2 (ik ) 2
(i ) (i )
− xη∗k ≤ xk k − γk Gk,ik (xk ) − x∗ηk + γk Gk,ik x∗ηk
k


xk+1 .

Combining the preceding relation with (4.2), we obtain:

d 2 2
X (i) (i) (i ) (ik )
D xk+1 , x∗ηk ≤ p−1 xk − x∗ηk + p−1 xk k − x∗ηk


i ik
i=1, i6=ik
 T
(ik ) ∗(ik )
− 2 p−1 Gk,ik (xk ) − Gk,ik x∗ηk



ik γ k x k − x ηk

2
+ p−1 2
ik γk Gk,ik (xk ) − Gk,ik xηk
∗
.

Invoking Definition 2.5, from the preceding relation we obtain:

(ik ) T
 
(ik )
D xk+1 , x∗ηk ≤ D xk , x∗ηk − 2 p−1 − x∗ηk Gk,ik (xk ) − Gk,ik x∗ηk





ik γ k x k

2
+ p−1 2
ik γk Gk,ik (xk ) − Gk,ik xηk
∗
.

Taking the conditional
expectation from the both sides of preceding relation and
noting that D xk , x∗ηk is Fk –measurable, we obtain the following inequality:
h
2i
E D xk+1 , x∗ηk |Fk ≤ D xk , x∗ηk + γk2 E p−1 Gk,ik (xk ) − Gk,ik x∗ηk



ik
  T 
−1 (ik ) ∗(ik ) ∗



(4.3) − 2γk E pik xk − xηk Gk,ik (xk ) − Gk,ik xηk .

Next, we estimate the second and third expectations in the preceding relation:
  T 
−1 (ik ) ∗(ik ) ∗



E p ik x k − x ηk Gk,ik (xk ) − Gk,ik xηk
d
(i) T
 
(i)
X
pi p−1 xk − x∗ηk Gk,i (xk ) − Gk,i x∗ηk



= i
i=1

T
= xk − x∗ηk Gk (xk ) − Gk x∗ηk



(4.4) .

We can also write:

h i
2
E p−1 Gk,ik (xk ) − Gk,ik x∗ηk


ik
d
X 2 2
pi p−1 Gk,i (xk ) − Gk,i x∗ηk = Gk (xk ) − Gk x∗ηk



(4.5) = i .
i=1
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 19

From Assumption 4.8, taking into account that Gk is (ηk µf )–strongly monotone, and
combining (4.3), (4.4), and (4.5) we obtain:
2
E D xk+1 , x∗ηk |Fk ≤ D xk , x∗ηk − 2µf γk ηk xk − x∗ηk



 
2
+ 2γk2 L2F + ηk2 L2f xk − x∗ηk + BF .

From Remark 2.6 and that {ηk } is a nonincreasing sequence, we obtain:

E D xk+1 , x∗ηk |Fk ≤ 1 − 2µf γk ηk pmin + 2γk2 pmax L2F + η0 2 L2f D xk , x∗ηk



+ 2γk2 BF .
µf ηk pmin
From the assumption γk ≤ 2pmax (L2F +η02 L2f )
and the preceding inequality, we obtain:

E D xk+1 , x∗ηk |Fk ≤ (1 − µf γk ηk pmin ) D xk , x∗ηk + 2γk2 BF .




(4.6)

The preceding relation is not yet fully recursive as the term x∗ηk on the right-hand
side must change to x∗ηk−1 . Next, we find an upper bound for D xk , x∗ηk in terms


 
of D xk , x∗ηk−1 . Note that we have ku + vk2 ≤ (1 + θ)kuk2 + 1 + θ1 kvk2 for any

vectors u, v ∈ Rn and θ > 0. Utilizing this inequality, by setting u := xk − x∗ηk−1 ,

p µ γ η
v := x∗ηk−1 − x∗ηk , and θ := min 2f k k we obtain:

2
 pmin µf γk ηk  2
xk − x∗ηk ≤ 1+ xk − x∗ηk−1
 2 
2 2
+ 1+ x∗ηk−1 − x∗ηk .
pmin µf γk ηk

Together with Lemma 4.5(b) and Remark 2.6, we have:

 pmin µf γk ηk   
pmin D xk , x∗ηk ≤ 1 + pmax D xk , x∗ηk−1
2
 2 2
C̄f

2 ηk−1
+ 1+ 1− .
pmin µf γk ηk µ2f ηk

Dividing both sides by pmin and substituting this in (4.6), we obtain:

 pmax  pmin µf γk ηk   
E D xk+1 , x∗ηk |Fk ≤ D xk , x∗ηk−1


(1 − γk ηk µf pmin ) 1 +
pmin 2
2 2
C̄f
 
2 ηk−1
+ 2 1+ 1− + 2γk2 BF .
µf pmin pmin µf γk ηk ηk
pmin µf γk ηk
pmin µf γk ηk
(4.1) is obtained by noting that (1 − γk ηk µf pmin ) 1 + 2 ≤ 1− 2 .
In the following result, we provide a class of update rules for the stepsize and
the regularization sequences such that Algorithm 2.1 attains both an almost sure
convergence and a convergence in the mean sense.
Theorem 4.11 (Convergence of Algorithm 2.1 when X is unbounded). Consider
(PfVI ).
Let the sequence {x̄k } be generated by Algorithm 2.1. Let Assumption 4.1,
Assumption 2.3, and Assumption 4.8 hold. Suppose the random block-coordinate ik
in Assumption 2.3 is drawn from a uniform distribution for all k ≥ 0. Let the stepsize
20 H. D. KAUSHIK AND F. YOUSEFIAN

{γk } and the regularization parameter {ηk } be given by γk := γ0 (k + 1)−a and ηk :=

η0 (k + 1)−b , respectively, where γ0 > 0, η0 > 0, 0 < b < 0.5 < a, and a + b < 1. Then,
the following results hold for all 0 ≤ r < 1:
(i) The sequence {x̄k } converges almost surely to the unique optimal solution of (PfVI ).
(ii) We have that limk→∞ E[kx̄k − x∗ k] = 0.
Proof. The proof is done in two main steps. In the first step, we show that
the non-averaged sequence {xk } converges to x∗ in an almost sure sense and that
limk→∞ E[kxk − x∗ k] = 0. In the second step, we show that these results hold for the
weighted average sequence {x̄k } as well.
Step 1: The proof of this step is done by applying Lemma 4.7 to the recursive inequal-
ity (4.1) with pi := d1 for all i ∈ {1, . . . , d}. The details are as follows. First, we note
that from the update rules of γk and ηk and that a > b, we have limk→∞ γηkk = 0. Thus,
µf pmin
there exists an integer k0 ≥ 1 such that for all k ≥ k0 , we have γηkk ≤ 2p 2 2 .
max (LF +η0 Lf )
2

This implies that the conditions of Lemma 4.10 are satisfied and the inequality (4.1)
holds for all k ≥ k0 . To apply Lemma 4.7, we define the following terms for all k ≥ 1:
  µf γk ηk
vk , D xk , x∗ηk−1 , αk , ,
!  2d 2
dC̄f2 (µf η0 γ0 + 2d) ηk−1
βk , − 1 + 2γk2 BF .
µ3f γk ηk ηk

Since γk ηk → 0, there exists an integer k1 ≥ k0 such that forPany k ≥ k1 we have

∞
0 ≤ αk ≤ 1. From the assumption that a + b < 1, we have that k=k1 αk = ∞. Next,
P∞
we show that k=k1 βk < ∞. From the update rules of γk and ηk and invoking the
Taylor series expansion, for k ≥ 2 we can write:
 b  
ηk−1 1 b b(b − 1) 1 b(b − 1)(b − 2) 1
−1= 1+ −1= 1+ + + + ... −1
ηk k k 2! k2 3! k3
  ∞
b (1 − b) (1 − b)(2 − b) (1 − b)(2 − b)(3 − b) bX 1
= 1− + − + ... ≤ ,
k 2!k 3!k 2 4!k 3 k i=0 k 2i

where the inequality is obtained using b < 1 and neglecting the negative terms. This
 2
4b 2 2
implies that ηk−1 b ηk−1
≤ 2b


ηk − 1 ≤ −2
k(1−k ) and thus ηk − 1 ≤ 3k k2 for all k ≥ 2.
Using the preceding relation, invoking the definition of βk ,
and the update formulas
βk = O k −(2−a−b) + O k −2a . From the assumptions


of γk and ηk , we have that P
∞
on a and b, we obtain that k=k1 βk < ∞. Also, from the assumption a > b, we get
limk→∞ βk /αk = 0. Therefore,
 all conditions of Lemma 4.7 are satisfied.
h  As such,
i we
have that D xk , x∗ηk−1 → 0 almost surely and also limk→∞ E D xk , x∗ηk−1 = 0.
From Remark 2.6 and that ik is drawn uniformly, we obtain:
2 2
2
kxk − x∗ k ≤ 2 xk − x∗ηk−1 + 2 x∗ηk−1 − x∗
2   2
(4.7) = D xk , x∗ηk−1 + 2 x∗ηk−1 − x∗ ,
d
where the first inequality is obtained from the triangle inequality. Taking the limit
 k → ∞and invoking Lemma 4.5(a),
from both sides of the preceding relation when
2
we obtain limk→∞ kxk − x∗ k ≤ d2 limk→∞ D xk , x∗ηk−1 . From the almost sure con-
 
vergence of D xk , x∗ηk−1 to zero, we conclude that {xk } converges to x∗ almost
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 21

surely. To show the convergence in mean, let us take the expectation from both
sides of (4.7). Noting that the Tikhonov trajectory is deterministic, we obtain that
h i h  i 2
2
E kxk − x∗ k ≤ d2 E D xk , x∗ηk−1 + 2 x∗ηk−1 − x∗ . Now, taking the limit from
both sides of thehpreceding
 relation
i when k → ∞, invoking Lemma h 4.5(a), and
i re-
∗ ∗ 2
calling limk→∞ E D xk , xηk−1 = 0, we conclude that limk→∞ E kxk − x k = 0.
Invoking Jensen’s inequality, we can conclude that limk→∞ E[kxk − x∗ k] = 0.
Step 2: Invoking Lemma 2.11 and using the triangle inequality, we have:
k
X k
X k
X
(4.8) kx̄k − x∗ k = λt,k xt − x∗ = λt,k (xt − x∗ ) ≤ λt,k kxt − x∗ k ,
t=0 t=0 t=0
Pk
where λt,k , γtr / j=0 γjr . In view of Lemma 4.6, let us define ut , kxt − x∗ k,
Pk P∞
vk , t=0 λt,k kxt − x∗ k, and αt , γtr . Note that since ar ≤ 1, we have t=0 αt =
∞
γ0r t=0 (t + 1)−ar = ∞. Also, from Step 1 we have that û , limt→∞ ut = 0 in
P
an almost sure sense. Thus, from Lemma 4.6, we conclude that {vk } converges to
zero almost surely. Thus, (4.8) implies that {x̄k } converges to x∗ almost surely.
Next, we apply Lemma 4.6 again, but in a slightly different fashion to show that
limk→∞ E[kx̄k − x∗ k] = 0. From (4.8), we have:
k
X
(4.9) E[kx̄k − x∗ k] ≤ λt,k E[kxt − x∗ k] .
t=0
Pk
In view of Lemma 4.6, let us define ut , E[kxt − x∗ k], vk , t=0 λt,k E[kxt − x∗ k],
and αt , γtr . First, note that from Step 1, we have û , limt→∞ ut = 0.
In view of Lemma 4.6, limk→∞ vk = 0. Thus, from (4.9), we conclude that
limk→∞ E[kx̄k − x∗ k] = 0. Hence, the proof is completed.
5. Experimental results. In this section, we revisit the problem of finding the
best Nash equilibrium formulated as in (1.1). We consider a case where the Nash
game is characterized as a Cournot competition over a network. Cournot game is
one of the most extensively studied economic models for competition among multiple
firms, including imperfectly competitive power markets as well as rate control over
communication networks [19, 21, 13]. Consider a collection of d firms who compete
to sell a commodity over a network with J nodes. The decision of each firm i ∈
{1, . . . , d} includes variables yij and sij , denoting the generation and sales of the firm
i at the node j, respectively. Considering the definitions yi , (yi1 ; . . . ; yiJ ) and
si , (si1 ; . . . ; siJ ), we can compactly denote the decision variable of the ith firm as
x(i) , (yi ; si )
∈ R2J . The goal of the ith firm lies in minimizing the net cost function
gi x(i) , x(−i) over the network defined as follows:

  XJ J
X
gi x(i) ;x(−i) , cij (yij ) − sij pj (s̄j ) ,
j=1 j=1

where cij : R → R denotes the production cost function of the firm i at the node
Pd
j, s̄j , i=1 sij denotes the aggregate sales from all the firms at the node j, and
pj : R → R denotes the price function with respect to the aggregate sales s̄j at the node
j. We assume that the cost functions are linear and the price functions are given as
σ
pj (s̄j ) , αj − βj (s̄j ) where σ ≥ 1 and αj and βj are positive scalars. Throughout,
22 H. D. KAUSHIK AND F. YOUSEFIAN

(γ0 , η0 ) =
ln (sample ave. gap) (0.1, 0.1) (0.1, 1) (1, 0.1)

10 10 10
ln(infeasibility)

ln(infeasibility)

ln(infeasibility)
8 8 8
6 6 6
4 Alg. 1.1 4 4
Alg. 2.1 (r=0)
2 Alg. 2.1 (r=0.5)
Alg. 2.1 (r=1)
2 2
0 50 100 150 200 250 300 0 50 100 150 200 250 300 0 50 100 150 200 250 300
time (seconds) time (seconds) time (seconds)
ln (sample ave. obj.)

10 10 10
9 9 9
8 8 8
ln(Objective)

ln(Objective)

ln(Objective)
7 f 7 7
6 g1 6 6
5 g2
g3
5 5
4 g4 4 4
3 50 100 150 200 250 300 3 50 100 150 200 250 300 3 50 100 150 200 250 300
time (seconds) time (seconds) time (seconds)

Fig. 1: Algorithm 2.1 in terms of infeasibility and the objective function value

we assume that the transportation costs are negligible. We let the generation be
capacitated as yij ≤ Bij , where Bij is a positive scalar for i ∈ {1, . . . , d} and j ∈
{1, . . . , J}. Lastly, for any firm i, the total sales must match with the total generation.
Consequently, the strategy set of the firm i is given as follows:
 
 XJ XJ 
Xi , (yi ; si ) | yij = sij , yij , sij ≥ 0, yij ≤ Bij , for all j = 1, . . . , J .
 
j=1 j=1

Following the model (1.1), we employ the Marshallian aggregate surplus function
Pd

defined as f (x) , i=1 gi x(i) ; x(−i) . We note that the convexity of the function f
is implied by σ ≥ 1 and the monotonicity of mapping F is guaranteed when either
σ = 1, or when 1 < σ ≤ 3 and d ≤ 3σ−1 σ−1 (cf. section 4 in [21]).
The set-up: In the experiment, we consider a Cournot game among 4 firms over 3
nodes. We let the slopes of the linear cost functions take values between 10 and 50. We
assume that αj := 50 and βj := 0.05 for all j, Bij := 120 for all i and j, and σ := 1.01.
To report the performance of Algorithm 2.1 in terms of the suboptimality, we plot a
sample average approximation of E[f (x̄N )] using the sample size of 25. With regard
to the infeasibility, we compute a sample average approximation of E[GAP (x̄N )] using
the same sample size. Following Remark 3.4, we use γk := √γk+1 0
and ηk := √ 4
η0
k+1
. To
select the block-coordinates in Algorithm 2.1, we use a discrete uniform distribution.
Results and insights: Figure 1 shows the experimental results. Here, in the top
three figures, we compare the performance of Algorithm 2.1 with that of Algorithm 1.1
in terms of infeasibility measured by the sample averaged gap function. Importantly,
the proposed algorithm performs significantly better than the SR scheme. This claim
is supported by considering the different values of the parameter r and the initial
conditions of the proposed scheme in terms of the initial stepsize γ0 and the initial
regularization parameter η0 . The three figures in the bottom row of Figure 1 demon-
strate the performance of Algorithm 2.1 in terms of reaching a stability in the objective
values. This includes the Marshallian objective function f as well as the individual
objective functions gi . Note that all the objective values in Figure 1 appear to reach
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 23

to a desired level of stability after around 60 seconds. This interesting observation

could be linked to the impact of the averaging scheme (2.2). Generally, it is expected
that the trajectories of the objective function values in Figure 1 be noisy due to the
randomness in the block-coordinate selection rule. However, the weighted averaging
scheme employed in Algorithm 2.1 appears to induce much robustness with respect to
this uncertainty, resulting in an accelerated convergence for the proposed algorithm.
6. Conclusions. Motivated by the applications arising from noncooperative
multi-agent networks, we consider a class of optimization problems with Cartesian
variational inequality (CVI) constraints. The computational complexity of the solu-
tion methods for addressing this class of problems appears to be unknown. We develop
a single timescale algorithm equipped with non-asymptotic suboptimality and infeasi-
bility convergence rates. Moreover, in the case where the set associated with the CVI
is unbounded, we establish the global convergence of the sequence generated by the
proposed algorithm. We apply the method in finding the best Nash equilibrium in
a networked Cournot competition. Our experimental results show that the proposed
method outperforms the classical sequential regularized schemes.
Appendix A. Additional proofs.
A.1. Proof of Lemma 1.3. Let us define the function φ : Rn → R as φ(x) ,
1 2 1
PJ 2
2 kAx−bk + 2 j=1 (max{0, hj (x)}) . We first note that φ is a differentiable function
such that ∇φ(x) = F (x) where F is given by Lemma 1.3 (e.g., see page 380 in [6]).
Next, we also note that φ is convex. To see this, note that from the convexity of

2
hj (x), the function h+ j (x) , max{0, hj (x)} is convex. Then, the function hj (x)
+

can be viewed as a composition of s(u) , u2 for u ∈ R and the convex function

h+ +
j . Since hj is nonnegative on its domain and s(u) is nondecreasing on [0, +∞),

2
we have that h+ j (x) is a convex function. As such, φ is a convex function as well.
Consequently, from the first-order optimality conditions for convex programs, we have
SOL(X, F ) = argminx∈X φ(x). To show the desired equivalence between problems
(PfVI ) and (1.2), it suffices to show that X = argminx∈X φ(x) where X denotes the
feasible set of problem (1.2). To show this statement, first we let x̄ ∈ X . Then, from
the definition of φ(x), we have φ(x̄) = 0. This implies that x̄ ∈ argminx∈X φ(x).
Thus, we have X ⊆ argminx∈X φ(x). Second, let x̃ ∈ argminx∈X φ(x). The feasibility
assumption of the set X implies that there exists an x0 ∈ X such that Ax0 = b and
hj (x0 ) ≤ 0 for all j. This implies that φ(x0 ) = 0. From the nonnegativity of φ and
that x̃ ∈ argminx∈X φ(x), we must have φ(x̃) = 0 and x̃ ∈ X. Therefore, we obtain
Ax̃ = b, hj (x̃) ≤ 0 for all j, and x̃ ∈ X. Thus, we have argminx∈X φ(x) ⊆ X . Hence,
we conclude that X = argminx∈X φ(x) = SOL(X, F ) and the proof is completed.
PN
A.2. Proof of Lemma 2.11. We use induction to show x̄N = k=0 λk,N xk
for any N ≥ 0. For N = 0, the relation holds due to the initialization x̄0 := x0 in
Algorithm 2.1 and that λ0,0 = 1. Next, let the relation hold for some N ≥ 0. From
PN
the hypothesis, equation (2.2), and that SN = k=0 γkr for all N ≥ 0, we can write:
r PN +1 r N +1
SN x̄N + γN +1 xN +1 k=0 γk xk
X
x̄N +1 = = PN +1 r
= λk,N +1 xk ,
SN +1 k=0 γk k=0

implying that the induction hypothesis holds for N + 1. Thus, we conclude that the
desired
PN averaging formula holds for all N ≥ 0. To complete the proof, note that since
k=0 λk,N = 1, under the convexity of the set X, we have x̄N ∈ X.
24 H. D. KAUSHIK AND F. YOUSEFIAN

A.3. Proof of Lemma 2.13. (a) From Definition 2.12, we can write:
d
X d
X
E[∆k | Fk ] = F (xk ) − pi p−1
i Ui Fi (xk ) = F (xk ) − Ui Fi (xk ) = 0.
i=1 i=1

The relation E[δk | Fk ] = 0 can be shown in a similar fashion.

(b) We can write:
d
 X 2
E k∆k k2 | Fk = pi F (xk ) − p−1

i Ui Fi (xk )
i=1
d
X  
2
= pi kF (xk )k2 + p−2 −1 T
i kUi Fi (xk )k − 2pi F (xk ) Ui Fi (xk )
i=1
d
X d
X
2
= kF (xk )k2 + p−1 kFi (xk )k2 ≤ p−1

2
i kUi Fi (xk )k − 2 min − 1 CF .
i=1 i=1

The relation E kδk k2 | Fk ≤ p−1

 
2
min − 1 Cf can be shown using a similar approach.

A.4. Proof of Lemma 2.14. Given 0 ≤ α < 1, let us define the function
φ : R++ → R as φ(x) , x−α for all x > 0. Since α > 0, the function φ is nonincreasing.
We can write:
N N +1 Z N +1
X 1 X 1 dx (N + 1)1−α − 1 (N + 1)1−α
α
= 1 + α
≤ 1 + α
=1+ ≤ ,
(k + 1) k 1 x 1−α 1−α
k=0 k=2

implying the desired upper bound. To show that the lower bound holds, we can write:
N N +1 Z N +2 Z N +1
X 1 X 1 dx dx (N + 1)1−α − 0.5(N + 1)1−α
= ≥ ≥ ≥ ,
(k + 1)α kα 1 xα 1 xα 1−α
k=0 k=1
1
where the last inequality is obtained using the assumption that N ≥ 2 1−α − 1. There-
fore, the desired lower bound holds as well. This completes the proof.
A.5. Proof of Lemma 4.5. (a) From the definition of x∗ and x∗ηk (cf. Defini-
tion 4.3), we have that:
(A.1) F (x∗ )T (x − x∗ ) ≥ 0 for all x ∈ X,
T
F x∗ηk + ηk ∇f x∗ηk y − x∗ηk ≥ 0





(A.2) for all y ∈ X.
For x := x∗ηk and y := x∗ , adding the resulting two relations together, we obtain:

T

T ∗
ηk ∇f x∗ηk x∗ − x∗ηk ≥ F (x∗ ) − F x∗ηk x − x∗ηk .

From the monotonicity of the mapping F and the preceding relation, we obtain that

T ∗
∇f x∗ηk x − x∗ηk ≥ 0. Also, from the strong convexity of f , we have:

T ∗
µf 2
f (x∗ ) ≥ f x∗ηk + ∇f x∗ηk x − x∗ηk + x∗ − x∗ηk


.
2
From the preceding relations, we obtain:

µf 2
(A.3) f (x∗ ) ≥ f x∗ηk + x∗ − x∗ηk for all k ≥ 0.
2
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 25

Thus, f (x∗ ) ≥ f x∗ηk for all k ≥ 0. Recall that from Remark 4.2, under Assump-

tion 4.1, x∗ ∈ X and x∗ηk ∈ X both exist and are unique. Therefore, f x∗ηk is

bounded above for all k ≥ 0. From this statement and invoking the coercive property
of f (implied by the strong convexity of f ), we can conclude that {x∗ηk } is a bounded
sequence. Therefore, it must have at least one limit point. Let {x∗ηk }k∈K be an ar-
bitrary subsequence such that limk→∞, k∈K x∗ηk = x̂. We show that x̂ ∈ SOL(X, F ).
Taking the limit from both sides of (A.2) with respect to the aforementioned sub-
sequence and using the continuity of F and ∇f , we obtain that for all y ∈ X,
T
(F (x̂) + limk→∞, k∈K ηk ∇f (x̂)) (y − x̂) ≥ 0. Note that the mapping ∇f (x̂) is
bounded. This is because x̂ ∈ X (due to the closedness of X) and that ∇f is con-
tinuous on the set X. Therefore, from the preceding inequality and limk→∞ ηk = 0,
T
we obtain F (x̂) (y − x̂) ≥ 0 for all y ∈ X, implying that x̂ ∈ SOL(X, F ) and so,
x̂ is a feasible solution to (PfVI ). Next, we show that x̂ is the optimal solution to
µ 2
(PfVI ). From (A.3), continuity of f , and neglecting the term 2f x∗ − x∗ηk , we ob-
tain f (x∗ ) ≥ f limk→∞, k∈K x∗ηk = f (x̂). Hence, from the uniqueness of x∗ , all the

limit points of {x∗ηk } fall in the singleton {x∗ } and the proof is completed.
(b) If x∗ηk = x∗ηk−1 , the desired relation holds. Suppose for k ≥ 1, we have x∗ηk 6= x∗ηk−1 .
From x∗ηk−1 ∈ SOL (X, F + ηk−1 ∇f ) and x∗ηk ∈ SOL (X, F + ηk ∇f ), we have that:
    T  
F x∗ηk−1 + ηk−1 ∇f x∗ηk−1 x − x∗ηk−1 ≥ 0 for all x ∈ X,
T
F x∗ηk + ηk ∇f x∗ηk y − x∗ηk ≥ 0





for all y ∈ X.

Adding the resulting two relations together, for x := x∗ηk and y := x∗ηk−1 we have:
  T  
−F x∗ηk − ηk ∇f x∗ηk + F (x∗ηk−1 ) + ηk−1 ∇f x∗ηk−1 x∗ηk − x∗ηk−1 ≥ 0.

  T  
The monotonicity of F implies that F x∗ηk − F x∗ηk−1 x∗ηk − x∗ηk−1 ≥ 0.

Adding this relation to the preceding inequality, we have:

  T  
ηk ∇f x∗ηk − ηk−1 ∇f x∗ηk−1 x∗ηk−1 − x∗ηk ≥ 0.

 T  
Adding and subtracting the term ηk ∇f x∗ηk−1 x∗ηk−1 − x∗ηk , we obtain:

 T  
(ηk − ηk−1 ) ∇f x∗ηk−1 x∗ηk−1 − x∗ηk ≥
  
 T  ∗ 
(A.4) ηk ∇f x∗ηk−1 − ∇f x∗ηk xηk−1 − x∗ηk .

From the strong convexity of function f , we have:

  
 T  ∗  2
(A.5) ∇f x∗ηk−1 − ∇f x∗ηk xηk−1 − x∗ηk ≥ µf x∗ηk − x∗ηk−1 .

From (A.4) and (A.5), we can write:

 T   2
(ηk − ηk−1 ) ∇f x∗ηk−1 x∗ηk−1 − x∗ηk ≥ ηk µf x∗ηk − x∗ηk−1 .
26 H. D. KAUSHIK AND F. YOUSEFIAN

Using the Cauchy-Schwarz inequality, we obtain:

  2
|ηk − ηk−1 | ∇f x∗ηk−1 x∗ηk−1 − x∗ηk ≥ ηk µf x∗ηk − x∗ηk−1 ,

Since x∗ηk 6= x∗ηk−1 , dividing the both sides by ηk x∗ηk − x∗ηk−1 , we obtain:

ηk−1  
(A.6) 1− ∇f x∗ηk−1 ≥ µf x∗ηk − x∗ηk−1 .
ηk

From part (a), the trajectory {x∗ηk } is bounded. Also, for any k ≥ 0, x∗ηk ∈ X by the
definition. Since X is closed, there exists a compact set S⊂X such that {x∗ηk } ⊂ S.
 and the continuity of ∇f imply that there exists C̄f > 0 such that
This statement
∇f x∗ηk−1 ≤ C̄f for all k ≥ 1. Thus, from (A.6), we obtain the desired inequality.

A.6. Proof of Corollary 3.5. Let us rewrite (PfVI ) as the equivalent problem:

minimize f (x)
(A.7)
subject to x ∈ SOL(Y, F ),
Qd0
where Y , i=1 Yi and d0 , 1 and Y1 , X. Note that this setting immediately
implies that Y = X. Now, let us consider Algorithm 2.1 for solving (A.7) where
we assume that x0 ∈ X is an arbitrary fixed vector. Since d0 = 1, Assumption 2.3
holds with Prob (ik = 1) = 1 for all k ≥ 0. This setting implies that Algorithm 2.1
reduces to a deterministic scheme where the step 5 in Algorithm 2.1 is equivalent to
the following update rule:
  
(A.8) ˜ (xk ) ,
xk+1 := PX xk − γk F (xk ) + ηk ∇f

where we used Y = Y1 = X. Next, we note that from Qdthe properties of the Euclidean
projection mapping, for any z ∈ X where X , i=1 Xi , we have that PX (z) =
Qd (i)


i=1 PXi z . In view of this property, the equation (A.8) compactly represents the
d updates given by (3.21). Therefore, Algorithm 3.1 is equivalent to Algorithm 2.1
and thus, all the results in Theorem 3.3 will hold with pmin = 1. Note that in both
(3.18) and (3.19), the expectation is eliminated. This completes the proof.

REFERENCES

[1] T. Alpcan and T. Başar, A game-theoretic framework for congestion control in general
topology networks, in Proceedings of the 41st IEEE Conference on Decision and Control,
December 2002, pp. 1218–1224.
[2] T. Alpcan and T. Başar, Distributed algorithms for Nash equilibria of flow control games,
in Advances in Dynamic Games, vol. 7 of Annals of the International Society of Dynamic
Games, Birkhäuser Boston, 2003, pp. 473–498.
[3] M. Amini and F. Yousefian, An iterative regularized mirror descent method for ill-posed
nondifferentiable stochastic optimization, (2019), [Link]
[4] E. Anshelevich, A. Dasgupta, J. Kleinberg, E. Tardos, T. Wexler, and T. Roughgar-
den, The price of stability for network design with fair cost allocation, SIAM Journal on
Computing, 38 (2008), pp. 1602—-1623.
[5] A. Beck and S. Sabach, A first order method for finding minimal norm-like solutions of
convex optimization problems, Mathematical Programming, 147 (2014), pp. 25–46.
[6] D. P. Bertsekas, Nonlinear Programming, Athena Scientific, Bellmont, MA, 3th ed., 2016.
 A METHOD WITH COMPLEXITY FOR OPTIMIZATION WITH VI CONSTRAINTS 27

[7] Y. Censor, A. Gibali, and S. Reich, The subgradient extragradient method for solving varia-
tional inequalities in Hilbert space, Journal of Optimization Theory and Applications, 148
(2011), pp. 318–335.
[8] Y. Censor, A. Gibali, and S. Reich, Extensions of Korpelevich’s extragradient method for the
variational inequality problem in Euclidean space, Optimization, 61 (2012), pp. 1119–1132.
[9] Y. Chen, G. Lan, and Y. Ouyang, Accelerated schemes for a class of variational inequalities,
Mathematical Programming, 165 (2017), pp. 113–149.
[10] J. R. Correa, A. S. Schulz, and N. E. Stier-Moses, Selfish routing in capacitated networks,
Mathematics of Operations Research, 29 (2004), pp. 961–976.
[11] C. D. Dang and G. Lan, On the convergence properties of non-Euclidean extragradient meth-
ods for variational inequalities with generalized monotone operators, Computational Op-
timization and Applications, 60 (2015), pp. 277–310.
[12] C. D. Dang and G. Lan, Stochastic block mirror descent methods for nonsmooth and stochastic
optimization, SIAM Journal on Optimization, 25 (2015), pp. 856–881.
[13] F. Facchinei and J.-S. Pang, Finite-dimensional Variational Inequalities and Complemen-
tarity Problems. Vols. I,II, Springer Series in Operations Research, Springer-Verlag, New
York, 2003.
[14] G. Garrigos, L. Rosasco, and S. Villa, Iterative regularization via dual diagonal descent,
Journal of Mathematical Imaging and Vision, 60 (2018), pp. 189–215.
[15] A. N. Iusem, A. Jofré, R. I. Oliveira, and P. Thompson, Extragradient method with variance
reduction for stochastic variational inequalities, SIAM Journal on Optimization, 27 (2016),
pp. 686–724.
[16] A. N. Iusem, A. Jofré, and P. Thompson, Incremental constraint projection methods for
monotone stochastic variational inequalities, Mathematics of Operations Research, 44
(2019), pp. 236–263.
[17] A. N. Iusem and M. Nasri, Korpelevich’s method for variational inequality problems in Banach
spaces, Journal of Global Optimization, 50 (2011), pp. 59–76.
[18] H. Jiang and H. Xu, Stochastic approximation approaches to the stochastic variational in-
equality problem, IEEE Transactions on Automatic Control, 53 (2008), pp. 1462–1475.
[19] R. Johari, Efficiency Loss in Market Mechanisms for Resource Allocation, PhD thesis, MIT,
2004.
[20] A. Juditsky, A. Nemirovski, and C. Tauvel, Solving variational inequalities with stochastic
mirror-prox algorithm, Stochastic Systems, 1 (2011), pp. 17–58.
[21] A. Kannan and U. V. Shanbhag, Distributed computation of equilibria in monotone Nash
games via iterative regularization techniques, SIAM Journal on Optimization, 22 (2012),
pp. 1177–1205.
[22] A. Kannan, U. V. Shanbhag, and H. M. Kim, Strategic behavior in power markets under
uncertainty, Energy Systems, 2 (2011), pp. 115–141.
[23] A. Kannan, U. V. Shanbhag, and H. M. Kim, Addressing supply-side risk in uncertain power
markets: stochastic Nash models, scalable algorithms and error analysis, Optimization
Methods and Software, 28 (2013), pp. 1095–1138.
[24] H. Kaushik and F. Yousefian, A randomized block coordinate iterative regularized subgradient
method for high-dimensional ill-posed convex optimization, in Proceedings of the American
Control Conference, IEEE, July 2019, pp. 3420–3425, [Link]
8815256, [Link]
[25] K. Knopp, Theory and applications of infinite series, Blackie & Son Ltd., Bishopbriggs, Glas-
gow G64 2NZ, Scotland, 1951.
[26] G. M. Korpelevich, An extragradient method for finding saddle points and for other problems,
Eknomika i Matematicheskie Metody, 12 (1976), pp. 747—-756.
[27] J. Koshal, A. Nedić, and U. V. Shanbhag, Regularized iterative stochastic approximation
methods for stochastic variational inequality problems, IEEE Transactions on Automatic
Control, 58 (2013), pp. 594–609.
[28] J. Lei, U. V. Shanbhag, J.-S. Pang, and S. Sen, On synchronous, asynchronous, and ran-
domized best-response schemes for stochastic Nash games, Mathematics of Operations
Research, 45 (2020), pp. 157–190.
[29] C. E. Lemke and J. T. Howson Jr., Equilibrium points of bimatrix games, Journal of the
Society for Industrial and Applied Mathematics, 12 (1964), pp. 413–423.
[30] P. Marcotte and D. Zhu, Weak sharp solutions of variational inequalities, SIAM Journal on
Optimization, 9 (1998), pp. 179–189.
[31] A. Nedić, Random algorithms for convex minimization problems, Mathematical Programming,
129 (2011), pp. 225–253.
[32] A. Nemirovski, Prox-method with rate of convergence O(1/t) for variational inequalities with
28 H. D. KAUSHIK AND F. YOUSEFIAN

Lipschitz continuous monotone operators and smooth convex-concave saddle point prob-
lems, SIAM Journal on Optimization, 15 (2004), pp. 229–251.
[33] YU. Nesterov, Efficiency of coordinate descent methods on huge-scale optimization problems,
SIAM Journal on Optimization, 22 (2012), pp. 341–362.
[34] N. Nisan, T. Roughgarden, E. Tardos, and V. V. Vazirani, Algorithmic Game Theory,
Cambridge University Press, New York, NY, USA, 2007.
[35] M. J. Osborne and A. Rubinstein, A Course in Game Theory, MIT Press, Cambridge,
Massachusetts, 1994.
[36] B. T. Polyak, Introduction to Optimization, Optimization Software, Inc., New York, 1987.
[37] P. Ricktárik and M. Takáč, Iteration complexity of randomized block-coordinate descent
methods for minimizing a composite function, Mathematical Programming, 144 (2014),
pp. 1–38.
[38] R. T. Rockafellar and R. J.-B. Wets, Variational Analysis, Springer-Verlag Berlin Heidel-
berg, 1998.
[39] T. Roughgarden, Stackelberg scheduling strategies, SIAM Journal on Computing, 33 (2004),
pp. 332–350.
[40] S. Sabach and S. Shtern, A first order method for solving convex bilevel optimization prob-
lems, SIAM Journal on Optimization, 27 (2017), pp. 640–660.
[41] H. Scarf, The approximation of fixed points of a continuous mapping, SIAM Journal on
Applied Mathematics, 15 (1967), pp. 1328–1343.
[42] G. Scutari, D. P. Palomar, F. Facchinei, and J.-S. Pang, Convex optimization, game
theory, and variational inequality theory, IEEE Signal Processing Magazine, 27 (2010),
pp. 35–49.
[43] G. Scutari, D. P. Palomar, F. Facchinei, and J.-S. Pang, Monotone games for cognitive
radio systems, (2012), pp. 83–112.
[44] S. Shalev-Shwartz and T. Zhang, Stochastic dual coordinate ascent methods for regularized
loss minimization, Journal of Machine Learning Research, 14 (2013), pp. 567–599.
[45] U. V. Shanbhag, G. Infanger, and P. W. Glynn, A complementarity framework for forward
contracting under uncertainty, Operations Research, 59 (2011), pp. 810–834.
[46] M. V. Solodov, An explicit descent method for bilevel convex optimization, Journal of Convex
Analysis, 14 (2007), pp. 227–237.
[47] J. Wang, G. Scutari, and D. P. Palomar, Robust MIMO cognitive radio via game theory,
IEEE Transactions on Signal Processing, 59 (2011), pp. 1183–1201.
[48] M. Wang and D. P. Bertsekas, Incremental constraint projection methods for variational
inequalities, Mathematical Programming, 150 (2015), pp. 321–363.
[49] H.-K. Xu, Viscosity approximation methods for nonexpansive mappings, Journal of Mathe-
matical Analysis and Applications, 298 (2004), pp. 279–291.
[50] H. Yin, U. V. Shanbhag, and P. G. Mehta, Nash equilibrium problems with scaled conges-
tion costs and shared constraints, IEEE Transactions on Automatic Control, 56 (2011),
pp. 1702–1708.
[51] F. Yousefian, Bilevel distributed optimization in directed networks, in Proceedings of the
American Control Conference (accepted), 2021, [Link]
[52] F. Yousefian, A. Nedić, and U. V. Shanbhag, Optimal robust smoothing extragradient al-
gorithms for stochastic variational inequality problems, in Proceedings of the 53rd IEEE
Conference on Decision and Control, IEEE, Dec. 2014, pp. 5831–5836.
[53] F. Yousefian, A. Nedić, and U. V. Shanbhag, On smoothing, regularization, and av-
eraging in stochastic approximation methods for stochastic variational inequality prob-
lems, Mathematical Programming, 165 (2017), pp. 391–431, [Link]
s10107-017-1175-y.
[54] F. Yousefian, A. Nedić, and U. V. Shanbhag, On stochastic mirror-prox algorithms for
stochastic Cartesian variational inequalities: Randomized block coordinate and optimal
averaging schemes, Set-Valued and Variational Analysis, 26 (2018), pp. 789–819, https:
//[Link]/10.1007/s11228-018-0472-9.
[55] F. Yousefian, A. Nedić, and U. V. Shanbhag, On stochastic and deterministic quasi-Newton
methods for nonstrongly convex optimization: Asymptotic convergence and rate analy-
sis, SIAM Journal on Optimization, 30 (2020), pp. 1144–1172, [Link]
17M1152474.