A New Implementable Prediction-Correction Method for Monotone Variational Inequalities with Separable Structure

Abstract

The monotone variational inequalities capture various concrete applications arising in many areas. In this paper, we develop a new prediction-correction method for monotone variational inequalities with separable structure. The new method can be easily implementable, and the main computational effort in each iteration of the method is to evaluate the proximal mappings of the involved operators. At each iteration, the algorithm also allows the involved subvariational inequalities to be solved in parallel. We establish the global convergence of the proposed method. Preliminary numerical results show that the new method can be competitive with Chen's proximal-based decomposition method in Chen and Teboulle (1994).

1. Introduction

The variational inequality (VI (Ω,𝐹)) in the finite-dimensional space is to determine a vector 𝑢∈Ω such that 𝑢^′−𝑢,𝐹(𝑢)≥0,∀𝑢^′∈Ω, (1) where Ω∈ℜ^𝑛 is a nonempty closed convex subset and 𝐹 is a continuous mapping from ℜ^𝑛 into itself. The VI (Ω,𝐹) has found many efficient applications in a broad spectrum of areas such as traffic equilibrium [1] and network economic problems [2]. For solving (1), the proximal point algorithm (PPA), which was proposed by Martinet [3] and further studied by Rockafellar [4, 5], generates the new iterative point 𝑢^𝑘+1 via the following procedure: 𝑢^′−𝑢^𝑘+1,𝐹𝑢^𝑘+1+𝐺𝑢^𝑘+1−𝑢^𝑘≥0,∀𝑢^′∈Ω, (2) where 𝐺∈ℜ^𝑛×𝑛 is a positive definite matrix, playing the role of proximal regularization parameter. Note that the PPA has to solve systems of nonlinear equations in each iteration. In many cases, solving these equations is quite difficult. This difficulty has inspired the burst of approximate versions of the PPAs, in order to approximately solve (2) under certain “relative error.” These new methods include well-known-extragradient type methods (EGM) as special cases. Assume that 𝐹 is Lipschitz continuous; that is, there is 𝑙∈(0,1), such that 𝛽𝐹𝑢^𝑘−𝐹𝑢^𝑘≤𝑙𝑢^𝑘−𝑢^𝑘. (3) Then at each iteration EGM takes the following general form: 𝑢^′−𝑢,𝛽𝐹𝑢^𝑘+𝑢−𝑢^𝑘≥0,∀𝑢^′𝑢∈Ω,_′−𝑢^𝑘+1,𝛽𝐹(𝑢)+𝑢^𝑘+1−𝑢^𝑘≥0,∀𝑢^′∈Ω. (4)

In this paper, we consider the following variational inequalities: find a vector 𝑤∈𝒟 such that 𝑤^′−𝑤,𝐹(𝑤)≥0,∀𝑤^′∈𝒟, (5) with 𝑥𝑦(𝑤)𝑓(𝑥)𝑔𝑦,,𝑤∶=,𝐹∶=𝒟=𝑥,𝑦∣𝑥∈𝒳,𝑦∈𝒴,𝐴𝑥+𝐵𝑦=𝑏 (6) where 𝒟∈ℜ^𝑛+𝑝 is a nonempty closed convex subset and 𝑓∶𝒳→ℜ^𝑛 and 𝑔∶𝒴→ℜ^𝑝 are monotone operators. Problem (5) is referred to as a structured variational inequality (SVI) [6].

By attaching a Lagrange multiplier vector 𝜆∈ℜ^𝑚 to the linear constraints 𝐴𝑥+𝐵𝑦=𝑏, the VI problem (5) is converted into the following form: 𝑥^′−𝑥,𝑓(𝑥)−𝐴^𝑇𝜆𝑦^′−𝑦,𝑔𝑦−𝐵^𝑇𝜆𝜆^′−𝜆,𝐴𝑥+𝐵𝑦−𝑏≥0,∀𝑢^′∈Ω, (7) where Ω=𝒳×𝒴×ℜ^𝑚. (8) The compact form is 𝑢^′−𝑢,𝐹(𝑢)≥0,∀𝑢^′∈Ω, (9) with 𝑥𝑦𝜆(𝑢)𝑓(𝑥)𝑢∶=,𝐹∶=−𝐴^𝑇𝜆𝑔𝑦−𝐵^𝑇𝜆𝐴𝑥+𝐵𝑦−𝑏. (10)

For the purpose of parallel computing, the proximal alternating directions method (PADM) generates 𝑢^𝑘=(𝑥^𝑘,𝑦^𝑘,𝜆^𝑘)∈Ω as follows [7, 8]: first find an 𝑥^𝑘∈𝒳 such that 𝑥^′−𝑥^𝑘,𝑓𝑥^𝑘−𝐴^𝑇𝜆^𝑘−𝛽𝐴𝑥^𝑘+𝑦^𝑘𝑥−𝑏+𝑟^𝑘−𝑥^𝑘≥0,∀𝑥∈𝒳. (11) Then find an 𝑦^𝑘∈𝒴 such that 𝑦^′−𝑦^𝑘,𝑔𝑦^𝑘−𝐵^𝑇𝜆^𝑘−𝛽𝐴𝑥^𝑘+𝑦^𝑘𝑦−𝑏+𝑠^𝑘−𝑦^𝑘≥0,∀𝑦∈𝒴. (12) Finally, update 𝜆^𝑘 via 𝜆^𝑘=𝜆^𝑘−𝛽𝐴𝑥^𝑘+𝐵𝑦^𝑘−𝑏. (13) Here 𝑟≥0 and 𝑠≥0 are given proximal parameters; 𝛽≥0 is a given penalty parameter for the linearly constraints. Note that when 𝑟=𝑠=0 in (11)-(12), the classical alternating directions method (ADM) is recovered. To make the PADM (11)–(13) more efficient and flexible, some strategies have been developed. For example, allow 𝑟, 𝑠, and 𝛽 to vary from iteration to iteration according to certain strategies [8–10]; produce the new iterate based on the minor correction to the predictor. A simple and effective correction scheme is (see, e.g., [11, 12]) 𝑢^𝑘+1=𝑢^𝑘−𝛼_𝑘𝑢^𝑘−𝑢^𝑘, (14) where 𝛼_𝑘>0 is a chosen step size.

The PADM (11)–(13) is often easy to implement under the assumption that the decomposed subproblems have closed-form solutions or can be efficiently solved up to a high precision. However, in some cases, matrixes 𝐴 and 𝐵 are not identity matrices, and the two subproblems in PADM (11)-(12) are difficult to solve because the evaluation of (𝐴^𝑇𝐴+(1/𝛽)𝑓)⁻¹(𝐴𝜐) and (𝐵^𝑇𝐵+(1/𝛽)𝑔)⁻¹(𝐵𝜐) could be costly. To overcome this difficulty, we propose a new implementable prediction-correction method for the SVI. At each iteration, we first decompose the problem to two small problems with respect to 𝑥 and 𝑦, respectively. The two subproblems are all easy to solve under the assumption that the resolvent operators of 𝑓 and 𝑔 are easy to evaluate, where the resolvent operator of mapping 𝑇 is defined as (𝐼+𝜆𝑇)⁻¹(𝜐). Then, we update the Lagrange multipliers and make a correction step to ensure the algorithm's convergence.

The SVI has been studied extensively both in the theoretical frameworks and applications. Recently, Han [13] proposed a hybrid entropic proximal decomposition method for the SVI. Han's method is based on logarithmic-quadratic functions and combined with self-adaptive strategy. He [14] presented a parallel splitting augmented Lagrangian method which can be extended to solve the system of equilibrium problems with three separable operators. Xu et al. [15] proposed two classes of correction methods for the SVI in which the mapping 𝐹 does not have an explicit form. Besides, Xu and Wu [16] also studied a class of linearized proximal alternating direction methods and showed that the relaxation factor can have the same restriction region as for the general ADM. Yuan and Li [17] developed a logarithmic-quadratic-proximal- (LQP-) based decomposition method by applying the LQP terms to regularize the ADM subproblems; then Bnouhachem et al. [18] studied a new inexact LQP alternating direction method by solving a series of related systems of nonlinear equations.

The rest of this paper is organized as follows. In Section 2, we review some preliminaries which are useful for further analysis. In Section 3, we present the new implementable prediction-correction method for SVI, and the global convergence result is established. Numerical experiments and some conclusions are addressed in Sections 4 and 5, respectively.

2. Preliminaries

In this section, we make some standard assumptions and summarize some basic properties of VI which will be used in the subsequent discussions.

Assumption(A1)𝒳,𝒴 are simple closed convex sets. A set which is said to be simple means that the projection onto the set is easy to compute, where the projection of a point 𝜐 onto the closed convex set Ω, denoted by 𝑃_Ω(𝜐), is defined as the nearest point 𝑢∈Ω to 𝜐; that is, 𝑃_Ω(𝜐)=argmin{‖𝑢−𝜐‖∣𝑢∈Ω}. (15)(A2)The mapping 𝐹 is point-to-point, monotone, and continuous. A mapping 𝐹∶ℜ^𝑛→ℜ^𝑛 is said to be monotone on Ω if ⟨𝑢−𝜐,𝐹(𝑢)−𝐹(𝜐)⟩≥0,∀𝑢,𝜐∈Ω. (16)(A3)The solution set of SVI (Ω,𝐹), denoted by Ω^∗, is nonempty.

Properties. Let 𝐺 be a symmetric positive definite matrix; the 𝐺-norm of the vector 𝑢 is denoted by ‖𝑢‖_𝐺∶=⟨𝑢,𝐺𝑢⟩. In particular, when 𝐺=𝐼, ‖𝑢‖∶=⟨𝑢,𝑢⟩ is the Euclidean norm of 𝑢. For a matrix 𝐴, ‖𝐴‖ denotes its norm ‖𝐴‖∶=max{‖𝐴𝑥‖∶‖𝑥‖≤1}.

The following well-known properties of the projection operator will be used in the coming analysis.

Lemma 1. Let Ω∈ℜ^𝑛 be a nonempty closed convex set; let 𝑃_Ω(⋅) be the projection operator onto Ω under the 𝐺-norm. Then 𝑢^′−𝑃_Ω𝑢^′,𝐺𝑢−𝑃_Ω𝑢^′≤0,∀𝑢^′∈ℜ^𝑛𝑃,∀𝑢∈Ω,_Ω(𝑢)−𝑃_Ω𝑢^′_𝐺≤𝑢−𝑢^′_𝐺,∀𝑢,𝑢^′∈ℜ^𝑛,𝑢−𝑃_Ω𝑢^′2_𝐺≤𝑢−𝑢^′2_𝐺−𝑢^′−𝑃_Ω𝑢^′2_𝐺,∀𝑢^′∈ℜ^𝑛,∀𝑢∈Ω. (17)

For any arbitrary positive scalar 𝛽 and 𝑢∈Ω, let 𝑒(𝑢,𝛽) denote the residual function associated with the mapping 𝐹; that is, 𝑒𝑢,𝛽=𝑢−𝑃_Ω𝑢−𝛽𝐹(𝑢). (18)

Lemma 2. 𝑢^∗ is a solution of the SVI (Ω,𝐹) if and only if 𝑒(𝑢^∗,𝛽)=0 for any given positive constant 𝛽 (see [2, page 267]).

Lemma 3. Solving SVI (Ω,𝐹) (7) is equivalent to find a zero point of the mapping 𝑒𝑒𝑢,𝛽∶=𝑒¹𝑢,𝛽₂𝑒𝑢,𝛽₃=𝑢,𝛽𝑥−𝑃_𝒳𝑥−𝛽𝑓(𝑥)−𝐴^𝑇𝜆𝑦−𝑃_𝒴𝑦−𝛽𝑔𝑦−𝐵^𝑇𝜆𝛽.𝐴𝑥+𝐵𝑦−𝑏 (19)

3. The New Algorithm

In this section, we present a new prediction-correction method for SVI (Ω,𝐹) and show its global convergence. But, at the beginning, to make the algorithm more succinct, we first define some matrices: 1𝐻=𝑟𝐼000𝑠𝐼000𝛽𝐼1,𝑀=𝐼0𝑟𝐴^𝑇10𝐼𝑠𝐵^𝑇,00𝐼𝑄=𝑟𝐼0𝐴^𝑇0𝑠𝐼𝐵^𝑇100𝛽𝐼. (20)

Obviously, 𝐻 is a symmetric positive definite matrix whenever 𝑟>0,𝑠>0, and 𝛽>0, and we also have 𝑄=𝐻𝑀.

3.1. Description of the Algorithm

Algorithm 4. It is a prediction-correction-based algorithm for the SVI (Ω,𝐹).

Phase 1 (initialization step). Given a small number 𝜖>0, let 𝛾∈(0,2); matrixes 𝑄,𝑀 are defined in (20). Take 𝑢⁰∈ℜ^{𝑛+𝑝+𝑚}; set 𝑘=0. Choose the parameters 𝑟>0,𝑠>0, and 𝛽>0 such that 𝑟>2𝛽𝐴^𝑇𝐴,𝑠>2𝛽𝐵^𝑇𝐵. (21)

Phase 2 (prediction step). Generate the predictor 𝑥^𝑘 via solving the following projection equation: 𝑥^𝑘=𝑃_𝒳𝑥^𝑘−1𝑟𝑓𝑥^𝑘−𝐴^𝑇𝜆^𝑘. (22) Then find an 𝑦^𝑘∈𝒴 such that 𝑦^𝑘=𝑃_𝒴𝑦^𝑘−1𝑠𝑔𝑦^𝑘−𝐵^𝑇𝜆^𝑘. (23) Finally, update 𝜆^𝑘 via 𝜆^𝑘=𝜆^𝑘−𝛽𝐴𝑥^𝑘+𝐵𝑦^𝑘−𝑏. (24)

Phase 3 (correction step). Correct the predictor, and generate the new iterate 𝑢^𝑘+1 via 𝑢^𝑘+1=𝑢^𝑘−𝛼_𝑘𝑀𝑢^𝑘−𝑢^𝑘, (25) where 𝛼_𝑘=𝛾𝛼^∗_𝑘,𝛼^∗_𝑘=𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘𝑀𝑢^𝑘−𝑢^𝑘2_𝐻. (26)

Phase 4 (convergence verification). If ‖𝑢^𝑘−𝑢^𝑘+1‖≤𝜖, stop; otherwise set 𝑘∶=𝑘+1; go to Phase 2.

Remark 5. Note that (22) does not involve 𝑦^𝑘 and that (23) is independent on the 𝑥^𝑘 generated by (22). Hence the two projections (22) and (23) are eligible for parallel computation.

Remark 6. It is easy to check that 𝑢^𝑘=(𝑥^𝑘,𝑦^𝑘,𝜆^𝑘) is a solution of SVI (Ω,𝐹) if and only if 𝐴𝑥^𝑘=𝐴𝑥^𝑘,𝐵𝑦^𝑘=𝐵𝑦^𝑘, and 𝜆^𝑘=𝜆^𝑘. Thus, it is reasonable to take the magnitude of ‖𝑢^𝑘−𝑢^𝑘+1‖≤𝜖 as the stopping criterion.

Remark 7. The strategy of choosing the step size 𝛼_𝑘 in the correction step which coincides with the strategy in He's papers, see, for example, [19], will be explained in detail in the following section.

Remark 8. Our method and the methods proposed in [6, 15, 20] are all in the prediction-correction algorithmic framework, where at each iteration they make a prediction step to produce a predictor and a correction step to generate the new iterate via correcting this predictor.

3.2. Contractive Properties

Now, we start to prove some properties of the sequence {𝑢^𝑘}. The first lemma quantifies the discrepancy between the point 𝑢^𝑘 and a solution point of SVI (Ω,𝐹).

Lemma 9. Let {𝑢} be generated by (22)–(24), and let the matrix 𝑀 be given in (20). Then one has 𝑢^′−𝑢^𝑘,𝐹(𝑢)−𝑄𝑢^𝑘−𝑢^𝑘≥0,∀𝑢^′∈Ω. (27)

Proof. Note that 𝑢^𝑘 generated by (22)–(24) are actually solutions of the following VIs: 𝑥^′−𝑥^𝑘,𝑓𝑥^𝑘−𝐴^𝑇𝜆^𝑘−𝑟𝑥^𝑘−𝑥^𝑘𝑦≥0,∀𝑥∈𝒳,_′−𝑦^𝑘,𝑔𝑦^𝑘−𝐵^𝑇𝜆^𝑘−𝑠𝑦^𝑘−𝑦^𝑘𝜆≥0,∀𝑦∈𝒴,_′−𝜆^𝑘,𝐴𝑥^𝑘+𝐵𝑦^𝑘1−𝑏−𝛽𝜆^𝑘−𝜆^𝑘≥0,∀𝜆∈ℜ^𝑚. (28)(29)(30) Combining (28)–(30) together, we have 𝑥^′−𝑥^𝑘,𝑓𝑥^𝑘−𝐴^𝑇𝜆^𝑘−𝐴^𝑇𝜆^𝑘−𝜆^𝑘−𝑟𝑥^𝑘−𝑥^𝑘𝑦^′−𝑦^𝑘,𝑔𝑦^𝑘−𝐵^𝑇𝜆^𝑘−𝐵^𝑇𝜆^𝑘−𝜆^𝑘−𝑠𝑦^𝑘−𝑦^𝑘𝜆^′−𝜆,𝐴𝑥^𝑘+𝐵𝑦^𝑘1−𝑏−𝛽𝜆^𝑘−𝜆^𝑘≥0,∀𝑢^′∈Ω. (31) Using the notations of 𝐹 (see (10)) and 𝑄 (see (20)), the earlier inequality can be rewritten into 𝑢^′−𝑢^𝑘,𝐹(𝑢)−𝑄𝑢^𝑘−𝑢^𝑘≥0,∀𝑢^′∈Ω. (32)
The assertion (27) is thus proved.

The following lemma plays a key role in proving the convergence of the algorithm.

Lemma 10. Let matrixes 𝑄, 𝐻 be defined in (20), if the parameters 𝑟>0,𝑠>0, and 𝛽>0 in (22)–(24) satisfy 𝑟>2𝛽𝐴^𝑇𝐴,𝑠>2𝛽𝐵^𝑇𝐵. (33) Then for the matrix 𝑄 in (27), one has ⟨(𝑢≥𝜇𝑢−𝑢,𝑄𝑢−)⟩1−2‖𝑢−𝑢‖²_𝐻∀𝑢≠𝑢∈ℜ^{𝑛+𝑝+𝑚}, (34) with 𝜇=𝐴max2𝛽^𝑇𝐴𝑟,2𝛽𝐵^𝑇𝐵𝑠∈(0,1). (35)

Proof. For any 𝑢≠𝑢, we have ⟨𝑢−𝑢,𝑄(𝑢−𝑢)⟩=‖𝑢−𝑢‖²_𝐻+𝜆−𝜆,𝐴(𝑥−𝑥)+𝜆−𝜆,𝑦−𝑦. (36) According to the Cauchy-Schwarz inequality, we get (𝑥)=1⟨𝜆−𝜆,𝐴𝑥−⟩+⟨𝜆−𝜆,𝑦−𝑦⟩22(𝑥𝑦1𝜆−𝜆,𝐴𝑥−)+2𝜆−𝜆,𝐵𝑦−≥−22𝛽𝜇‖𝐴(𝑥−𝑥)‖²+𝜇𝜆2𝛽𝜆−²−122𝛽𝜇𝐵𝑦−𝑦²+𝜇𝜆2𝛽𝜆−²1=−22𝛽𝜇‖𝐴(𝑥−𝑥)‖²+2𝛽𝜇×𝐵𝑦−𝑦²+𝜇𝛽𝜆−𝜆². (37) With the 𝜇 defined in (35), we have 2𝛽𝜇‖𝐴(𝑥−𝑥)‖²≤𝜇𝑟‖𝑥−𝑥‖²,2𝛽𝜇𝐵𝑦−𝑦²≤𝜇𝑠𝑦−𝑦². (38) Substituting (38) into (37), combining (36), the assertion (34) is proved.

Lemma 11. Suppose that 𝑢^∗=(𝑥^∗,𝑦^∗,𝜆^∗)∈Ω is a solution point of (9) and the sequences {𝑢^𝑘+1} are corrected by an undeterminate step size denoted by 𝛼 instead of (26); that is, 𝑢^𝑘+1=𝑢^𝑘−𝛼𝑀𝑢^𝑘−𝑢^𝑘. (39) Then one has 𝜗^𝑘(𝛼)≥𝜑^𝑘(𝛼), (40) where 𝜗^𝑘(𝛼)=𝑢^𝑘−𝑢^∗2_𝐻−𝑢^𝑘+1−𝑢^∗2_𝐻,𝜑^𝑘(𝛼)=2𝛼𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘−𝛼²𝑀𝑢^𝑘−𝑢^𝑘2_𝐻. (41)

Proof. One can see that 𝜗^𝑘(𝛼)=𝑢^𝑘−𝑢^∗2_𝐻−𝑢^𝑘+1−𝑢^∗2_𝐻=𝑢^𝑘−𝑢^∗2_𝐻−𝑢^𝑘−𝛼𝑀𝑢^𝑘−𝑢^𝑘−𝑢^∗2_𝐻=2𝛼𝑢^𝑘−𝑢^∗,𝐻𝑀𝑢^𝑘−𝑢^𝑘−𝛼²𝑀𝑢^𝑘−𝑢^𝑘2_𝐻. (42) On the other hand, since 𝑄=𝐻𝑀, using the monotonicity of 𝐹 and Lemma 9, we have 𝑢^𝑘−𝑢^∗,𝐻𝑀𝑢^𝑘−𝑢^𝑘=𝑢^𝑘−𝑢^∗,𝑄𝑢^𝑘−𝑢^𝑘=𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘+𝑢^𝑘−𝑢^∗,𝑄𝑢^𝑘−𝑢^𝑘≥𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘+𝑢^𝑘−𝑢^∗,𝐹𝑢^𝑘≥𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘+𝑢^𝑘−𝑢^∗,𝐹𝑢^∗≥𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘. (43) Combining (42)-(43) together, we have 𝜗^𝑘(𝛼)=2𝛼𝑢^𝑘−𝑢^∗,𝑄𝑢^𝑘−𝑢^𝑘−𝛼²𝑀𝑢^𝑘−𝑢^𝑘2_𝐻≥2𝛼𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘−𝛼²𝑀𝑢^𝑘−𝑢^𝑘2_𝐻=𝜑^𝑘(𝛼). (44) Thus, 𝜑^𝑘(𝛼) is a lower bound of 𝜗^𝑘(𝛼) for any 𝛼>0.

Remark 12. Note that 𝜑^𝑘(𝛼) is a quadratic function of 𝛼 and it reaches its maximum at 𝛼^∗_𝑘=𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘𝑀𝑢^𝑘−𝑢^𝑘2_𝐻. (45) Hence, it is reasonable to use the step size strategy (26). The parameter 𝛾 in (26) plays the role of a relaxation or scaling parameter. We can easily see that 𝛾∈(0,2) can ensure convergence.

Now, we prove the Fejér monotonicity of the iterative sequence {𝑢^𝑘} generated by the algorithm.

Theorem 13. Suppose that 𝑢^∗=(𝑥^∗,𝑦^∗,𝜆^∗)∈Ω is a solution point of (9) and the sequences {𝑢^𝑘} are generated by the algorithm. Then 𝑢^𝑘+1−𝑢^∗2_𝐻≤𝑢^𝑘−𝑢^∗2_𝐻−12𝑟()𝜇2−𝑟1−2𝑢^𝑘−𝑢^𝑘2_𝐻. (46)

Proof. According to Lemma 11, 𝑢^𝑘+1−𝑢^∗2_𝐻≤𝑢^𝑘−𝑢^∗2_𝐻−𝜑^𝑘𝛼_𝑘=𝑢^𝑘−𝑢^∗2_𝐻−2𝛼_𝑘𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘−𝛼²_𝑘𝑀𝑢^𝑘−𝑢^𝑘2_𝐻=𝑢^𝑘−𝑢^∗2_𝐻−𝛾2−𝛾𝛼^∗_𝑘×𝑢^𝑘−𝑢^𝑘,𝑄𝑢^𝑘−𝑢^𝑘≤𝑢^𝑘−𝑢^∗2_𝐻−𝛾2−𝛾𝛼^∗_𝑘𝜇1−2‖𝑢−𝑢‖²_𝐻. (47) Moreover, it follows from (26) that the step size 𝛼^∗_𝑘=𝑢^𝑘−𝑢^𝑘2_{𝑄+𝑄𝑇}2𝑀𝑢^𝑘−𝑢^𝑘2_𝐻=𝑢^𝑘−𝑢^𝑘2_{𝑄+𝑄𝑇}2𝑢^𝑘−𝑢^𝑘2_{𝑀𝑇𝐻𝑀}=12𝑢^𝑘−𝑢^𝑘2_{𝑀𝑇𝐻𝑀}+𝑟𝑥^𝑘−𝑥^𝑘2+𝑠𝑦^𝑘−𝑦^𝑘2+1𝛽𝜆^𝑘−𝜆^𝑘2−1𝑟𝐴^𝑇𝜆^𝑘−𝜆^𝑘2−1𝑠𝐵^𝑇𝜆^𝑘−𝜆^𝑘2×𝑢^𝑘−𝑢^𝑘2_{𝑀𝑇𝐻𝑀}⁻¹≥12𝑢^𝑘−𝑢^𝑘2_{𝑀𝑇𝐻𝑀}+𝑟𝑥^𝑘−𝑥^𝑘2+𝑠𝑦^𝑘−𝑦^𝑘2+1𝛽−1𝑟𝐴^𝑇𝐴−1𝑠𝐵^𝑇𝐵𝜆^𝑘−𝜆^𝑘2×𝑢^𝑘−𝑢^𝑘2_{𝑀𝑇𝐻𝑀}⁻¹. (48) Based on the conditions (33), we have 𝑟𝑥^𝑘−𝑥^𝑘2+𝑠𝑦^𝑘−𝑦^𝑘2+1𝛽−1𝑟𝐴^𝑇𝐴−1𝑠𝐵^𝑇𝐵𝜆^𝑘−𝜆^𝑘2≥0. (49) Hence, 𝛼^∗_𝑘≥12𝑢^𝑘−𝑢^𝑘2_{𝑀𝑇𝐻𝑀}𝑢^𝑘−𝑢^𝑘2_{𝑀𝑇𝐻𝑀}=12. (50) Substituting (50) into (47), we have 𝑢^𝑘+1−𝑢^∗2_𝐻≤𝑢^𝑘−𝑢^∗2_𝐻−12𝑟()𝜇2−𝑟1−2𝑢^𝑘−𝑢^𝑘2_𝐻. (51) Thus, we obtain the assertion of this theorem.

Based on the earlier results, we are now ready to prove the global convergence of the algorithm.

Theorem 14. The sequence {𝑢^𝑘} generated by the proposed algorithm converges to a solution of SVI (Ω,𝐹).

Proof. We prove the convergence of the proposed algorithm by following the standard analytic framework of contraction-type methods. It follows from (46) that {𝑢^𝑘} is bounded, and we have that lim_𝑘→∞𝑢^𝑘−𝑢^𝑘_𝐻=0. (52) Consequently, lim_𝑘→∞𝑥^𝑘−𝑥^𝑘=0,lim_𝑘→∞𝑦^𝑘−𝑦^𝑘=0,lim_𝑘→∞𝐴𝑥^𝑘+𝐵𝑦^𝑘−𝑏=lim_𝑘→∞1𝛽𝜆^𝑘−𝜆^𝑘=0, (53) since (see (22) and (23)) 𝑥^𝑘=𝑃_𝒳𝑥^𝑘−1𝑟𝑓𝑥^𝑘−𝐴^𝑇𝜆^𝑘+𝑥^𝑘−𝑥^𝑘+1𝑟𝐴^𝑇𝜆^𝑘−𝜆^𝑘,𝑦^𝑘=𝑃_𝒴𝑦^𝑘−1𝑠𝑔𝑦^𝑘−𝐵^𝑇𝜆^𝑘+𝑦^𝑘−𝑦^𝑘+1𝑠𝐵^𝑇𝜆^𝑘−𝜆^𝑘,𝜆^𝑘=𝜆^𝑘−𝛽𝐴𝑥^𝑘+𝐵𝑦^𝑘−𝑏. (54) It follows from (53) that lim_𝑘→∞𝑥^𝑘−𝑃_𝒳𝑥^𝑘−1𝑟𝑓𝑥^𝑘−𝐴^𝑇𝜆^𝑘=0,lim_𝑘→∞𝑦^𝑘−𝑃_𝒴𝑦^𝑘−1𝑠𝑔𝑦^𝑘−𝐵^𝑇𝜆^𝑘=0,lim_𝑘→∞𝐴𝑥^𝑘+𝐵𝑦^𝑘−𝑏=0. (55) Because 𝑢^𝑘 is also bounded, it has at least one cluster point. Let 𝑢^∞ be a cluster point of 𝑢^𝑘, and let 𝑢^𝑘𝑗 be the subsequence converging to 𝑢^∞. It follows from (55) that lim_𝑗→∞𝑥^𝑘𝑗−𝑃_𝒳𝑥^𝑘𝑗−1𝑟𝑓𝑥^𝑘𝑗−𝐴^𝑇𝜆^𝑘𝑗=0,lim_𝑗→∞𝑦^𝑘𝑗−𝑃_𝒴𝑦^𝑘𝑗−1𝑠𝑔𝑦^𝑘𝑗−𝐵^𝑇𝜆^𝑘𝑗=0,lim_𝑗→∞𝐴𝑥^𝑘𝑗+𝐵𝑦^𝑘𝑗−𝑏=0. (56) Consequently, 𝑥^∞−𝑃_𝒳𝑥^∞−1𝑟𝑓𝑥^∞−𝐴^𝑇𝜆^∞𝑦=0,_∞−𝑃_𝒴𝑦^∞−1𝑠𝑔𝑦^∞−𝐵^𝑇𝜆^∞=0,𝐴𝑥^∞+𝐵𝑦^∞−𝑏=0. (57) Using the continuity of 𝐹 and the projection operator 𝑃_Ω(⋅), we have that 𝑢^∞ is a solution of SVI (Ω,𝐹).
On the other hand, by taking limits over the subsequences in (52) and using lim_𝑗→∞𝑢^𝑘𝑗=𝑢^∞, we have that, for any 𝑘>𝑘_𝑗, it follows from (46) that 𝑢^𝑘−𝑢^∞_𝐻≤𝑢^𝑘𝑗−𝑢^∞_𝐻. (58) Thus, the sequence {𝑢^𝑘} converges to 𝑢^∞, which is a solution of SVI (Ω,𝐹).

4. Numerical Experiments

In this section, we present some numerical experiments results to show the effectiveness of the proposed algorithm. The codes are run on a notebook computer with Inter(R) Core(TM) 2 CPU 2.0 GHZ and RAM 2.00 GM under MATLAB Version 2009b.

We consider the following optimization problem: min_{𝑥∈ℜ𝑛,𝑦∈ℜ𝑝}12𝑥^𝑇1𝑃𝑥+2𝑦^𝑇𝑄𝑦s.t.𝐴𝑥+𝐵𝑦=𝑏, (59) where  𝑃∈ℜ^𝑛×𝑛,𝑄∈ℜ^𝑝×𝑝 are symmetric positive semidefinite matrixes, 𝐴∈ℜ^𝑚×𝑛, 𝐵∈ℜ^𝑚×𝑝,  and 𝑏∈ℜ^𝑚.

Using the KKT condition, the problem (59) can be converted into the following variational inequality: find 𝑤^∗=(𝑥^∗,𝑦^∗,𝜆^∗)∈ℜ^{𝑛+𝑝+𝑚} such that 𝑥^′−𝑥^∗,𝑃𝑥^∗−𝐴^𝑇𝜆^∗𝑦^′−𝑦^∗,𝑄𝑦^∗−𝐵^𝑇𝜆^∗𝜆^′−𝜆^∗,𝐴𝑥^∗+𝐵𝑦^∗−𝑏≥0,∀𝑤^′∈ℜ^{𝑛+𝑝+𝑚}. (60) In this example, we randomly created the input data of the tested collection in the following manner.(i)𝑃 and 𝑄 were generated randomly with eigenvalues in [5,10] according to the following MATLAB scripts: 𝑃=rand(𝑛);[𝑄1,𝑅1]=qr(𝑃), 𝑆=5+5∗rand(𝑛,1); 𝑃=𝑄1∗diag(𝑆)∗𝑄1^′, 𝑄=rand(𝑝); [𝑄2,𝑅2]=qr(𝑄), 𝑆=5+5∗rand(𝑚,1); 𝑃=𝑄2∗diag(𝑆)∗𝑄2^′.(ii)𝐴 and 𝐵 were generated randomly with singular values in [0,3], and the maximum singular value is 3 according to the following MATLAB scripts: 𝐴=rand(𝑚,𝑛);[𝑈,𝑆,𝑉]=svd(𝐴), 𝑆=𝑆/𝑆(1,1)∗3; 𝐴=𝑈∗𝑆∗𝑉^′, 𝐵=rand(𝑚,𝑝); [𝑈,𝑆,𝑉]=svd(𝐵), 𝑆=𝑆/𝑆(1,1)∗3; 𝐵=𝑈∗𝑆∗𝑉^′.(iii)𝑏 is generated randomly with 𝑏=rand(𝑚,1)∗10. According to the data generation, we have ‖𝐴^𝑇𝐴‖=9 and ‖𝐵^𝑇𝐵‖=9.

To apply (22)–(25) to solve (59), instead of choosing the step length 𝛼_𝑘 judiciously as (24), we can simply choose 𝛼_𝑘≡1 by takeing 𝑟=1/𝛼^∗_𝑘 (since 𝛼^∗_𝑘>1/2 when 𝑢≠𝑢, we have 𝑟∈(0,2) which satisfies the requirement). Then, we obtain the following subproblems which are all easy enough to have closed-form solutions: 𝑥^𝑘=(𝑟𝐼+𝑃)⁻¹𝑟𝑥^𝑘+𝐴^𝑇𝜆^𝑘,𝑦^𝑘=(𝑠𝐼+𝑄)⁻¹𝑠𝑦^𝑘+𝐵^𝑇𝜆^𝑘,𝜆^𝑘=𝜆^𝑘−𝛽𝐴𝑥^𝑘+𝐵𝑦^𝑘,𝑥−𝑏_𝑘+1=𝑥^𝑘−1𝑟𝐴^𝑇𝜆^𝑘−𝜆^𝑘,𝑦^𝑘+1=𝑦^𝑘−1𝑠𝐵^𝑇𝜆^𝑘−𝜆^𝑘,𝜆^𝑘+1=𝜆^𝑘. (61)

For comparison, we also solve it by the parallel decomposition method (denoted by PDM) that has been studied extensively in the literature (e.g., [21, 22]). For PDM, the restrictions on the proximal parameters are the same as our algorithm. By applying PDM to (59), we obtain the following subproblems which are also easy enough to have closed-form solutions: 𝑥^𝑘+1=(𝑟𝐼+𝑃)⁻¹𝑟𝑥^𝑘−𝛽𝐴^𝑇𝐴𝑥^𝑘+𝐵𝑦^𝑘1−𝑏−𝛽𝜆^𝑘,𝑦^𝑘+1=(𝑠𝐼+𝑄)⁻¹𝑠𝑦^𝑘−𝛽𝐵^𝑇𝐴𝑥^𝑘+𝐵𝑦^𝑘1−𝑏−𝛽𝜆^𝑘,𝜆^𝑘+1=𝜆^𝑘−𝛽𝐴𝑥^𝑘+1+𝐵𝑦^𝑘+1−𝑏. (62)

We report the numerical experiments by building their performance profiles in terms of the number of iterations and the total of computational time. Here, we take 𝛽=3+(𝑛/10),𝑟=𝑠=20𝛽 for the two algorithms. We set the initial vector (𝑥⁰,𝑦⁰,𝜆⁰)=(0,0,0), and the stopping criterion is Tol=max𝑥^𝑘+1−𝑥^𝑘_∞,𝑦^𝑘+1−𝑦^𝑘_∞,𝜆^𝑘+1−𝜆^𝑘_∞≤10⁻⁴. (63)

The computational results are given in Table 1 for different choices of 𝑚, 𝑛, and 𝑝. We reported the number of iterations (Iter.) and the computing time in seconds (CPU(s)) when the mentioned stopping criterion is achieved.

Table 1 

Numerical results for the example.

The data in Table 1 indicates clearly that the proposed method is efficient compared with the classical PDM in [21, 22]. We can observe that the iteration numbers and the CPU time of the two algorithms are almost the same.

5. Conclusions

In this paper, we proposed a new implementable algorithm for solving the monotone variational inequalities with separable structure. At each iteration, the algorithm performs two easily implementable projections parallelly to produce a predictor and then makes a simple correction to generate the new iterate. Under some mild conditions, we proved the global convergence of the new method. We also give some numerical experiments to show that the proposed method is applicable and valid.