Some Recent Trends in Variational Inequalities and Optimization Problems with Applications

and Applied Analysis Some Recent Trends in Variational Inequalities and Optimization Problems with Applications Guest Editors: Abdellah Bnouhachem, AbdelouahedHamdi, and Xu Minghua

Variational inequalities theory, which was introduced in the sixties, has emerged as an interesting and fascinating branch of applicable mathematics with a wide range of applications in industry, finance, economics, social, and pure and applied sciences. This field is dynamic and is experiencing an explosive growth in both theory and applications; as a consequence, research techniques and problems are drawn from various fields. The ideas and techniques of variational inequalities are being applied in a variety of diverse areas of sciences and prove to be productive and innovative. It has been shown that this theory provides the most natural, direct, simple, unified, and efficient framework for a general treatment of a wide class of unrelated linear and nonlinear problems. Variational inequalities have been extended and generalized in several directions using novel and new techniques. In parallel, optimization methods based on proximal point and proximal-like type methods have attracted a large number of researchers in the last three decades. In the same spirit, we can cite, for instance, the alternating direction multipliers method, which is based on the augmented lagrangian algorithm, which itself can be seen as a direct application of the proximal point algorithm to the dual problem of a constrained optimization problem.
The aim of this special issue is to present new approaches and theories for variational inequalities arising in mathematics and applied sciences. This special issue includes 14 highquality peer-reviewed papers that deal with different aspects of variational inequalities. These papers contain some new, novel, and innovative techniques and ideas. We hope that all the papers published in this special issue can motivate and foster further scientific works and development of the research in the area of theory, algorithms, and applications of variational inequalities.
The summaries of the 14 papers in this issue are listed as follows.
The paper of C. Chen et al. considers a class of linearly constrained separable convex programming problems without coupled variables. They weaken some conditions to obtain convergence of the alternating direction method of multipliers and they propose also a relaxed ADMM involving an additional computation of optimal step size and establish its global convergence under mild conditions. The paper of H. Sun and Y. Wang revisits the global error bound for the generalized nonlinear complementarity problem over a polyhedral cone (GNCP) and sharpens the global error bound for the GNCP under weaker conditions, which improves the existing error bound estimation for the problem.
The paper of M. Ma concerns the design and the convergence analysis of algorithms to split variational inequality and equilibrium problems.
The paper of Y. Wang and C. Wang gives a new modified Ishikawa type iteration algorithm for common fixed points of total asymptotically strict pseudocontractive semigroups.
Strong and weak convergence are proved under mild conditions. Furthermore, the main results presented in this work extend and improve some recent results.
The paper of X. Fu presents an implementable proximal step by a slight relaxation to the subproblem of proximal point algorithm (PPA) to solve linearly constrained convex programming. Self-adaptive strategies are proposed to improve the convergence in practice. The paper also discusses some applications and performs some numerical experiments to confirm the efficiency of the proposed method.
The paper of H. Xu establishes the strong convergence of prediction-correction and relaxed hybrid steepest-descent method (PRH method) for variational inequalities under some suitable conditions that simplify the proof. Further, the author shows the efficiency of the proposed algorithm through a well-designed set of practical numerical experiments.
The paper of M. Xu et al. considers the study of some matrix optimization problems using the proximal alternating direction method. The authors show that the restriction on the proximal parameters can be relaxed for solving these kinds of problems and give some numerical experiments to conclude that their modified method presents better performance than the classical proximal alternating direction method.
The paper of J.-L. Jiang et al. considers the locations of multiple facilities in the space R , with the aim of minimizing the sum of weighted distances between facilities and regional customers, where the proximity between a facility and a regional customer is evaluated by the closest distance. And the authors propose a new location-allocation heuristic scheme to solve their problem. Convergence is proved under mild assumptions; and furthermore some preliminary numerical results are reported to show the effectiveness of the new algorithm.
The paper of A. Roldan et al. studies the existence and uniqueness of coincidence point for nonlinear mappings of any number of arguments under a weak ( , )-contractivity condition in partial metric spaces. The obtained results generalize, extend, and unify several classical and very recent related results in the literature in metric spaces and in partial metric spaces.
The paper of Z. Jia et al. extends the convergence analysis given by Han and Yuan for alternating direction method of multipliers (ADMM) from the strongly convex to a more general case. Further, the authors prove under the assumption that the individual functions are composites of strongly convex functions and linear functions that the classical ADMM for separable convex programming with two blocks can be extended to the case with more than three blocks.
The paper of M. Li and Y. You presents a simple proof for the same convergence rate of the relaxed proximal point algorithm (PPA) in both ergodic and nonergodic senses.
The paper of W.-S. Du et al. extends, generalizes, and improves several fundamental results on the existence (and uniqueness) of coincidence points and fixed points for well-known maps in the literature. Furthermore, some fixed coincidence point theorems for multivalued nonself maps in the context of complete metric spaces are given.
The paper of F. Ma et al. develops, studies, and implements a new prediction-correction method for monotone variational inequalities with separable structure. At each iteration, the proposed algorithm also allows the involved subvariational inequalities to be solved in parallel.
The paper of A. Barbagallo and P. Mauro concerns a dynamic oligopolistic market equilibrium problem in the realistic case in which the presence of capacity constraints and production excesses are allowed and, moreover, the production function depends not only on the time but also on the equilibrium distribution. The authors prove the equivalence between this equilibrium definition and a suitable evolutionary quasi-variational inequality, and they study the analysis of existence, regularity, and sensitivity of solutions.

Introduction
This paper concerns the following problem: where ∈ × is a given symmetric matrix, matrices ∈ × and ∈ × are symmetric and scalars, and are the problem data, ⪰ 0 denotes that is a positive semidefinite matrix, Tr denotes the trace of a matrix, and ‖ ⋅ ‖ denotes the Frobenius norm; that is, 2 ) and + ∩ is nonempty. Throughout this paper, we assume that the Slater's constraint qualification condition holds so that there is no duality gap if we use Lagrangian techniques to find the optimal solution to problem (1). Problem (1) is a type of matrix nearness problem, that is, the problem of finding a matrix that satisfies some properties and is nearest to a given one. Problem (1) can be called the least squares covariance adjustment problem or the least squares semidefinite programming problem and solved by many methods [1][2][3][4]. In a least squares covariance adjustment problem, we make adjustments to a symmetric matrix so that it is consistent with prior knowledge or assumptions and a valid covariance matrix [2,5,6]. The matrix nearness problem has many applications especially in several areas of numerical linear algebra, finance industry, and statistics in [6]. A recent survey of matrix nearness problems can be found in [7]. It is clear that the matrix nearness problem considered here is a convex optimization problem. It thus follows from the strict feasibility and coercivity of the objective function that the minimum of (1) is attainable and unique.

Abstract and Applied Analysis
In the literature of interior point algorithms, + is called the semidefinite cone and the related problem (1) belongs to the class of semidefinite programming (SDP) and secondorder cone programming (SOCP) [8]. In fact, it is possible to reformulate problem (1) into a mixed SDP and SOCP as in [3,9]: where ⟨ , ⟩ = Tr( ). Thus, problem (1) can be efficiently solved by standard interior-point methods such as SeDuMi [10] and SDPT3 [11] when the number of variables (i.e., entries in the matrix ) is modest, say under 1000 (corresponds to around 32) and the number of equality and inequality constraints is not too large (say 5,000) [2,3,12].
Specially, let where Diag( ) is the vector of diagonal elements of and is the vector of 1s. Then problem (1) can be viewed as the nearest correlation matrix problem. For the nearest correlation matrix problem, a quadratically convergent Newton algorithm was presented recently by Qi and Sun [13], and improved by Borsdorf and Higham [1]. For problem (1) with equality and inequality constraints, one difficulty in finding an efficient method for solving this problem is the presence of the inequality constraints. In [3], Gao and Sun overcome this difficulty by reformulating the problem as a system of semismooth equations with two level metric projection operators and then design an inexact smoothing Newton method to solve the resulting semismooth system. For the problem (1) with large number of equality and inequality constraints, the numerical experiments in [14] show that the alternating direction method (hereafter alternating direction method is abbreviated as ADM) is more efficient in computing time than the inexact smoothing Newton method which additionally requires solving a large system of linear equations at each iteration. The ADM has many applications in solving optimization problems [15,16]. Papers written by Zhang, Han, Li, Yuan, and Bauschke and Borwein show that the ADM can be applied to solve convex feasibility problems [17][18][19].
The proximal ADM is a class of ADM type methods which can also be easily applied to solve the matrix optimization problems. Generally, the proximal parameters (i.e., the parameters and in (14) and (15)) of the proximal ADM are greater than zero. In this paper, we will show that the restriction on the proximal parameters can be relaxed while the proximal ADM is used to solve problem (1). Numerical experiments also show that the proximal ADM with the relaxed proximal parameters generally has a better performance than the classical proximal ADM.
The paper is organized as follows. In Section 2, we give some preliminaries about the proximal alternating direction method. In Section 3, we convert the problem (1) to a structured variational inequality and apply the proximal ADM to solve it. The basic analysis and convergent results of the proximal ADM with relaxed proximal parameters are built in Section 4. Preliminary numerical results are reported in Section 5. Finally, we give some conclusions in Section 6.
First, +1 is found by solving the following problem: Abstract and Applied Analysis 3 where ∈ X. Then, +1 is obtained by solving where ∈ Y. Finally, the multiplier is updated by where > 0 is a given penalty parameter for the linearly constraint + − = 0. Most of the existing ADM methods require that the subvariational inequality problems (11)- (12) should be solved exactly at each iteration. Note that the involved subvariational inequality problem (11)- (12) may not be well-conditioned without strongly monotone assumptions on and . Hence, it is difficult to solve these subvariational inequality problems exactly in many cases. In order to improve the condition of solving the subproblem by the ADM, some proximal ADMs were proposed (see, e.g., [26,27,[30][31][32][33][34]). The classical proximal ADM is one of the attractive ADMs. From a given triple = ( , , ) ∈ W, the classical proximal ADM produces the new iterate +1 = ( +1 , +1 , +1 ) ∈ W by the following procedure.
First, +1 is obtained by solving the following variational inequality problem: where > 0 is the given proximal parameter and ∈ X. Then, +1 is found by solving where > 0 is the given proximal parameter and ∈ Y. Finally, the multiplier is updated by In this paper, we will conclude that problem (1) can be solved by the proximal ADM and the restriction on the proximal parameters > 0, > 0 can be relaxed as > −1/2, > −1/2 when the proximal ADM is applied to solve problem (1). Our numerical experiments later also show that the numerical performance of the proximal ADM with smaller value of proximal parameters is generally better than the proximal ADM with comparatively larger value of proximal parameters.
Problem (29) is often a medium-scale quadratic programming (QP) problem. A variety of methods for solving the QP are commonly used, including interior-point methods and active set algorithm (see [36,37]).
Particularly, if is the following special case: where ≥ 0 expresses that each element of is nonnegative, and are given × symmetric matrices, and ≤ means that − ≥ 0; then ( ) is easy to be carried out and is given by where max( , ) and min( , ) compute the element-wise maximum and minimum of matrix and , respectively.

Main Results
Let { } be the sequence generated by applying the procedure (14)- (16) to problem (18)- (19); then for any = ( , , Λ ) ∈ W, we have that Further, letting where ∈ × is the unit matrix, and then we can get the following lemmas.
Proof. From (41), we have Abstract and Applied Analysis Rearranging the inequality above, we find that Using the Cauchy-Schwarz Inequality on the last term of the right-hand side of (49), we obtain Substituting (50) into (49), we get Thus, the proof is completed.
Based on the Theorem 3, we get the following lemma.
Proof. Since it is easy to check that if > −1/2, > −1/2, then and are symmetric positive-definite matrices. Let > 0 be the smallest eigenvalue of matrix . Then, from (46), we have Following (53), we immediately have that ‖ − * ‖ 2 is nonincreasing and thus the sequence { } is bounded. Moreover, we have Following point (3) of Lemma 4, we have This means that ∞ is a solution point of (18)- (19). Since { } converges to ∞ , we have that, for any given > 0, there exists an integer > 0 such that Furthermore, using the inequality (53), we have Combining (59) and (60), we get that This implies that the sequence { } converges to ∞ . So the proof is completed.

Numerical Experiments
In this section, we implement the proximal ADM to solve the problem (1) and show the numerical performances of proximal ADM with different proximal parameters. Additionally, we compare the classical ADM (i.e., the proximal ADM with proximal parameters = 0 and = 0) with the alternating projections method proposed by Higham [6] numerically and show that the alternating projections method is not equivalent to proximal ADM with zero proximal parameters. All the codes were written in Matlab 7.1 and run on IBM notebook PC R400.
Example 6. In the first numerical experiment, we set the 1 as an × matrix whose entries are generated randomly in [−1, 1]. Let = ( 1 + 1 )/2 and further let the diagonal elements of be 1 that is, = 1, = 1, 2, . . . , . In this test example, we simply let be in the form of (31) Moreover, let 0 = eye( ), 0 = eye( ), Λ 0 = zeroes( ), = 4, and = 10 −6 , where eye( ) and zeroes( ) are both the Matlab functions. For different problem size and different proximal parameters and , Table 1 shows the computational results. There, we report the number of iterations (It.) and the computing time in seconds (CPU.) it takes to reach convergence. The stopping criterion of the proximal ADM is where ‖ ‖ max = max(max(abs( ))) is the maximum absolute value of the elements of the matrix .
Remark 7. Note that if the proximal parameters are equal to zero, that is, = 0 and = 0, then the proximal ADM is the classical ADM.
The computational results are reported in Table 2.
Example 9. Let be in the form of (31) and = 0, = +∞, , = 1, 2, . . . , . Assume that , 0 , 0 , Λ 0 , , , and the stopping criterion are the same as those in Example 6, but the diagonal elements of matrix are replaced by = + (1 − ) × rand, = 1, 2, . . . , , where ∈ (0, 1) is a given number, rand is the Matlab function generating a number randomly in [0, 1]. In the following numerical experiments, we let = 0.2. For different problem size and different proximal parameters and , Table 3 shows the number of iterations and the computing time in seconds it takes to reach convergence. 8 Abstract and Applied Analysis    Example 10. All the data are the same as in Example 9 except that = 0. The computational results are reported in Table 4.
Example 11. Let 1 be an × matrix whose entries are generated randomly in [−0.5, 0.5], = ( 1 + 1 )/2, and let the diagonal elements of be 1. And let where B , B , B are subsets of {( , ) | 1 ≤ , ≤ } denoting the indexes of such entries of that are constrained by equality, lower bounds, and upper bounds, respectively. In this test example, we let the index sets B , B , and B be the same as in Example 5.4 of [3]; that is, B = {( , ) | = 1, 2, . . . , } and B , B ⊂ {( , ) | 1 ≤ < ≤ } consist of the indices of min(̂, − ) randomly generated elements at the th row of , = 1, 2, . . . , witĥ= 5 and̂= 10, respectively. We take = 1 for ( , ) ∈ B , = −0.1 for ( , ) ∈ B , and = 0.1 for ( , ) ∈ B . Moreover, let 0 , 0 , Λ 0 , , , and the stopping criterion be the same as those in Example 6. For different problem size , different proximal parameters and , and different values of̂, Tables 5(a) and 5(b) show the number of iterations and the computing time in seconds it takes to reach convergence, respectively. Numerical experiments show that the proximal ADM with relaxed parameters is convergent. Moreover, we draw the conclusion that the proximal ADM with smaller value of proximal parameters generally converges more quickly than the proximal ADM with comparatively larger value of proximal parameters to solve the problem (1).

Example 12.
In this test example, we apply the proximal ADM with = 0, = 0 (i.e., the classical ADM) to solve the nearest correlation matrix problem, that is, problem (1) with in the form of (5), and compare the classical ADM numerically with the alternating projections method (APM) [6]. The APM computes the nearest correlation matrix to a symmetric ∈ × by the following process: end.
In this numerical experiment, the stopping criterion of the APM is Let the matrix and the initial parameters of classical ADM be the same as those in Example 6. Table 6(a) reports the numerical performance of proximal ADM and the APM for computing the nearest correlation matrix to . Further, let 1 be an × matrix whose entries are generated randomly in [0, 1] and = ( 1 + 1 )/2. The other data are the same as above. Table 6(b) reports the numerical performance of the classical ADM and the APM for computing the nearest correlation matrix to the matrix . Numerical experiments show that the classical ADM generally exhibits a better numerical performance than the APM for the test problems above.

Conclusions
In this paper, we apply the proximal ADM to a class of matrix optimization problems and find that the restriction of proximal parameters can be relaxed. Moreover, numerical experiments show that the proximal ADM with relaxed parameters generally has a better numerical performance in solving the matrix optimization problem than the classical proximal alternating direction method.

Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.

Introduction
The aim of the paper is to improve the results obtained in [1] concerning the dynamic oligopolistic market equilibrium problem in presence of production excesses by introducing the dependence on the equilibrium commodity shipment in the production function (see K * ( * ) in (4)) and, as a consequence, studying the so-called elastic model. This is a more realistic situation since it is reasonable to think that the production function is influenced not only by the time, but also by the evaluation of the amount of commodity shipment, namely, the forecasted equilibrium solution. The presence of production excesses may be well justified in periods of economic crisis, so it is possible that some of the amounts of the commodity available are sold out whereas for a part of the products, an excess of production can occur. In the last decade a lot of problems considering a feasible set depending on equilibrium solutions have been studied (see, e.g., [2][3][4]). It is well known that the equilibrium models with fixed constraint sets may be expressed in terms of evolutionary variational inequalities, while models with elastic constraint sets are expressed by evolutionary quasi-variational inequalities. Moreover, the dependence on time leads to considering variational and quasi-variational inequalities in an infinite dimensional setting, for example, a Lebesgue space.
Let us remember that a dynamic oligopolistic market equilibrium problem is the problem of finding a trade equilibrium in a supply-demand market between a finite number of spatially separated firms producing homogeneous goods in a fixed time interval. Moreover, the firms act in a noncooperative behavior. This problem has its origin with Cournot [5]. He considered only two firms and for this reason it was called the duopoly problem. Later, Nash [6,7] extended Cournot's duopoly problem to agents. A more complete and efficient study was done by Nagurney et al. in [8][9][10][11], but the problem was still faced in a static case through a finite dimensional variational approach. Finally, in [12] the time dependence was considered in the model. It allows to explore the change of behavior of equilibrium states for the oligopolistic market models over a finite time interval of interest. As Beckmann and Wallace stressed, for the first time, in [13], "the timedependent formulation of equilibrium problems allows one to explore the dynamics of adjustment processes in which a delay on time response is operating. " Of course a delay on time response always happens because the processes do not have an infinite speed. Usually, such adjustment processes can be represented by means of a memory term which depends on previous equilibrium solutions according to the Volterra operator (see, e.g., [14,15]). Furthermore, in [16] the authors describe the behavior of the market by using the Lagrange multipliers of the infinite dimensional duality theory developed in [17][18][19][20][21]. Such results make use of the notion of tangent cone, normal cone, and quasi-relative interior of sets (see [22,23]), important tools to overcome the difficulty of the emptiness of the topological interior of the ordering cone which defines constraints of several infinite dimensional problems (see [24,25]). Moreover, a sensitivity result has been obtained which states that, under additional assumptions, small changes of the solution happen in correspondence with small changes of the profit function.
Lately, in [1,26], the model presented in [12] has been improved with the addition of production excesses and both production and demand excesses, respectively. Another important question is to find some regularity properties for the solution. In [1,26], the continuity of solution is proved under suitable assumptions, and it results to be very helpful in order to introduce numerical schemes to compute equilibrium solutions (see [27,28]).
In [29,30], the authors abandon the study of the problem from a producer's point of view whose purpose is to maximize his own profit and focus their attention on the policy-maker's perspective whose aim is to control the commodity exportations by means of the imposition of taxes or incentives and formulate the resulting optimization problem as an inverse variational inequality. This paper is structured as follows. In Section 2 we present the dynamic oligopolistic market equilibrium problem with the elastic production function and after that we give the definition of equilibrium according to the generalized Cournot-Nash principle. Moreover, we prove the equivalence with a suitable evolutionary quasi-variational inequality. Section 3 is devoted to prove a result of existence of the solution, while in Section 4 Kuratowski's set convergence will be a preliminary property in order to prove the continuity of the equilibrium solution. In Section 5 we establish a sensitivity result that shows how the equilibrium solution can change if the data have been perturbed. In Section 6 a numerical example is provided to make the theoretical model presented in the previous sections more clearer.

Quasi-Variational Inequalities in Dynamic Oligopolistic Markets
Let us consider firms , = 1, . . . , , that produce a homogeneous commodity and demand markets , = 1, . . . , , that are generally spatially separated. Assume that the homogeneous commodity, produced by the firms and consumed by the markets, is involved during a time interval [0, ], > 0.
We consider a formulation of equilibrium problems where the dependence of the production on the unknown solution * is in the average sense with respect to the time; namely, the following feasibility condition holds: Hence, condition (3) states that the average of the quantity produced by each firm , in the time interval [0, ], must be equal to the commodity shipments from that firm to all the demand markets plus the production excess, at the time ∈ [0, ]. In fact, the production is supposed to depend on the firms' evaluation of the commodity shipments. So one can expect the producers not to evaluate the market practicability instantly, but by an average with respect to the whole time interval.
Abstract and Applied Analysis Moreover, let us associate with each firm a production cost * , = 1, . . . , , and assume that the production cost of a firm may depend upon the entire production pattern; namely, * = * ( , ( ) , ( )) .
Similarly, let us associate with each demand market a demand price for unity of the commodity , = 1, . . . , , and assume that the demand price of a demand market may depend, in general, upon the entire consumption pattern; namely, = ( , ( )) .
For the reader's convenience, we recall that a function V, continuously differentiable, is called pseudoconcave with respect to , = 1, . . . , , a.e. in [0, ] (see [31]), if the following holds a.e. in [0, ]: Moreover, we recall that in the Hilbert space 2 ([0, ], R ), we define the canonical bilinear form on 2 ([0, ], where ∈ ( 2 ([0, ], R )) * = 2 ([0, ], R ), ∈ 2 ([0, ], R ), and Now, let us consider the dynamic oligopolistic market, in which the firms supply the commodity in a noncooperative fashion, each one trying to maximize its own profit at the time ∈ [0, ]. We seek to determine a nonnegative commodity distribution matrix function * for which the firms and the demand markets will be in a state of equilibrium according to the dynamic Cournot-Nash principle.
Proof. First of all, let us prove that the evolutionary quasivariational inequality (20), that we can write as follows: is equivalent to the following point-to-point quasi-variational inequality: where In fact, let us suppose by absurdum that (22) does not hold, namely, ∃ ( ) ∈ K( * ), ∃ ⊆ [0, ] with ( ) > 0 such that Let us choose, now, Abstract and Applied Analysis 5 Hence, let us consider that is a contradiction. The vice versa is immediate. So the equivalence between the evolutionary quasivariational inequalities (20) and (22) is proved.
Let us prove, now, the equivalence between the dynamic Cournot-Nash principle and the evolutionary quasi-variational inequality (20).
Since the profit function V ( , ( )) is pseudoconcave with respect to , we get If we choose ∈ K( * ) such that so we get the contradiction.

An Existence Theorem for Equilibrium Solutions
Now, we prove an existence result for the equilibrium solution to the dynamic elastic oligopolistic market equilibrium problem. To this aim, we recall a general existence result for solutions to quasi-variational inequalities in topological linear locally convex Hausdorff spaces due to Tan [32].
Now, we are able to prove our main result.
It is easy to verify that if ≤ ], for (47), ∈ K( ).   As a consequence, we have that the sequence { } ∈N converges to . It is easy to show that K( ) is a closed, bounded, and convex subset of and since the space is compact, K( ), for all ∈ , is compact too. As a consequence, all the hypotheses of Theorem 3 are satisfied and the existence of at least one solution is guaranteed.

Regularity Results for Equilibrium Solutions
In this section, we study the assumptions under which the continuity of solutions to evolutionary quasi-variational inequality, which expresses the equilibrium condition for the dynamic elastic oligopolistic market equilibrium problem in presence of production excesses, is ensured.

Set Convergence.
First of all, we recall the notion of Kuratowski's set convergence that has an important role in order to establish regularity results. The classical notion of convergence for subsets of a given metric space ( , ) is introduced in the 1950s by Kuratowski (see [33]; see also [34,35] where eventually means that there exists ∈ N such that ∈ K for any ≥ and frequently means that there exists an infinite subset ⊆ N such that ∈ K for any ∈ (in this last case, according to the notation given above, we also write that there exists a subsequence { } ∈N ⊆ { } ∈N such that ∈ K ). In the following, we recall Kuratowski's set convergence.
Definition 5. We say that {K } converges to some subset K ⊆ in Kuratowski's sense and we briefly write K → K, if − lim K = −lim K = K. Thus, in order to verify that K → K, it suffices to check that (i) − lim K ⊆ K, that is, for any sequence { } ∈N frequently in K such that → for some ∈ , then ∈ K; (ii) K ⊂ − lim K , that is, for any ∈ K there exists a sequence { } ∈N eventually in K such that → .
The below lemma establishes that the feasible set K of the dynamic elastic oligopolistic market equilibrium problem in the presence of production excesses satisfies the property of Kuratowski's set convergence.
Since (72) Passing to the limit in (71) as → +∞ and taking into account the continuity assumption on the functions , , and , we have Now, passing to the limit for → +∞ in the left-hand side of (72), we have Then, from (74) and (59), we obtain As a consequence, and, hence, condition (K2) is achieved.

Continuity of Solutions to Weighted Quasi-Variational
Inequalities. In [2,[36][37][38][39] some continuity results for variational and quasi-variational inequalities in infinite dimensional spaces have been obtained. It is worth remarking that similar results have been proved for weighted variational and quasi-variational inequalities in nonpivot Hilbert spaces (see [4,40]). Now, we show a continuity result for equilibrium solutions to the dynamic elastic oligopolistic market equilibrium problem in presence of production excesses. Theorem 7. Let , ∈ 0 ([0, ], R + ), and let ∈ 0 ([0, ] × R + , R + ) be such that Moreover, let V ∈ 1 ([0, ] × R + , R + ) be a vector function satisfying assumption (iii) and such that Then the dynamic elastic market equilibrium distribution in presence of production excesses * ∈ K( * ) is continuous in Proof. The existence of equilibrium solution is ensured by Theorem 4. Moreover, by applying Theorem 8 in [2] and taking into account Lemma 6, we obtain the continuity of * ∈ K( * ) in [0, ].

A Sensitivity Result
In this section a theorem about the sensitivity of solution is presented. The following result establishes that a small change in profit function produces a small change in equilibrium distribution.
Theorem 8. Assume that the profit function changes from V(⋅) to the perturbed functionṼ(⋅) and denote by * andt he correspondent solutions of the following quasi-variational inequalities: Let ∇ V( , ) be a strongly monotone function of constant , namely, for all , ∈ K( * ), ∃ > 0 such that Moreover, let ∇ V be a Carathéodory function such that By adding and subtracting −∇ V(̃) in (84), we have Moreover, by using the strong monotonicity, inequality (85), and the Cauchy-Schwartz inequality, we get (87)

A Numerical Example
This section is devoted to provide a numerical example of the theoretical achievements presented. Let us consider two firms and two demand markets, as in Figure 1. Let , ∈ 2 ([0, 1], R 4 ) be the capacity constraints such that, a.e. in Let ∈ 1 ([0, 1] × , R 2 ) be the production function, such that, a.e. in [0, 1], Abstract and Applied Analysis

11
The dynamic oligopolistic market equilibrium distribution in presence of the excesses is the solution to the evolutionary quasi-variational inequality: Let us observe that all the hypotheses of Theorem 4 are satisfied; hence, evolutionary quasi-variational inequality (94) admits solutions.
In order to compute a solution to (94) we make use of the direct method (see [41]). We consider the following system: ) . (96) Let us observe that the solution satisfies all the constraints; in particular, if we compute we are able to obtain the related production excesses: ( ) = ( 47 + 7 10 25 + 18 18 ) .

Conclusions
In [1] the dynamic oligopolistic market equilibrium problem was studied by introducing production excesses, and the dynamic Cournot-Nash equilibrium was characterized as a solution to a suitable evolutionary variational inequality. In this paper, in order to have a model closer to reality, it was supposed that the production function depends on the equilibrium commodity shipment. Hence, an elastic formulation was introduced that leads to an equivalent formulation by means of a suitable evolutionary quasi-variational inequality.
By means of this mathematical formulation, results of existence and regularity of solutions were proved. Furthermore, a sensitivity analysis is provided. At last a numerical example was provided in order to clarify the theoretical results. In future work, it is possible to consider also demand excesses and elastic demand function, in order to have a more complete and realistic model.

Introduction
Let K = {V ∈ | V ≥ 0, V = 0} be a polyhedral cone in for matrices ∈ × , ∈ × , and let K ∘ be its dual cone; that is, For continuous mappings , : → , the generalized nonlinear complementarity problem, abbreviated as GNCP, is to find vector * ∈ such that Throughout this paper, the solution set of the GNCP, denoted by * , is assumed to be nonempty. The GNCP is a direct generalization of the classical nonlinear complementarity problem and a special case of the general variational inequalities problem [1]. The GNCP was deeply discussed [2][3][4][5] after the work in [6]. The GNCP plays a significant role in economics, operation research, nonlinear analysis, and so forth (see [7,8]). For example, the classical Walrasian law of competitive equilibria of exchange economies can be formulated as a generalized nonlinear complementarity problem in the price and excess demand variables (see [8]).
For the GNCP, the solution existence and the numerical solution methods for the GNCP were discussed [2,3,6]. As an important tool for a mathematical problem, the global error bound estimation for GNCP with the mapping being -strongly monotone and Hölder continuous was discussed in [5], and a global error bound for the GNCP for the linear and monotonic case was established in [4].
In this paper, we will establish a global error bound for the problem (2) without the Hölder continuity of the underlying mapping. To this end, we first develop some new equivalent reformulations of the GNCP under weaker conditions and then establish a sharper global error bound for the GNCP in terms of some easier computed residual functions. The results obtained in this paper can be taken as an improvement of the existing results for GNCP and variational inequalities problem [4,5,[9][10][11].
To end this section, we give some notations used in this paper. Vectors considered in this paper are taken in the Euclidean space equipped with the usual inner product, and the Euclidean 2-norm and 1-norm of vector in are, respectively, denoted by ‖ ⋅ ‖ and ‖ ⋅ ‖ 1 . We use + to denote the nonnegative orthant in and use + and − to denote the vectors composed by elements ( + ) := max{ , 0}, ( − ) := max{− , 0}, 1 ≤ ≤ , respectively. For simplicity, we use ( ; ) to denote vector ( ⊤ , ⊤ ) ⊤ , use to denote the identity matrix with appropriate dimension, use ≥ 0 to denote 2 Abstract and Applied Analysis a nonnegative vector ∈ , and use dist( , * ) to denote the distance from point to the solution set * .

Global Error Bound for the GNCP
First, we give some concepts used in the subsequent.
(ii) -strongly -monotone with respect to : → if there are constants 1 > 0, > 0 such that Remark 2. Based on this definition, -strongly -monotone implies monotonicity, and if ( ) = + , ( ) = + with , ∈ × , , ∈ , then the above Definition 1(i) is equivalent to that the matrix ⊤ is positive semidefinite. Now, we give some assumptions for our analysis based on Definition 1.

Assumption 3.
For mappings , and matrix involved in the GNCP, we assume that (A1) mapping is monotone with respect to mapping ; (A2) matrix ⊤ has full-column rank.
In the following, we will establish a new equivalent reformulation to the GNCP. First, we give the following conclusion established in [2].
is a solution of the GNCP if and only if there exist * 1 ∈ , * 2 ∈ , such that ( * ) ≥ 0, From Theorem 5, under Assumption 3(A2), we can transform the system into a new system in which neither 1 nor 2 is involved. To this end, we need the following conclusion [12].

Lemma 6. If the linear system
= is consistent, then = + is the solution with the minimum 2-norm, where + is the pesudo-inverse of . Lemma 7. Suppose that Assumption 3(A2) holds. Then, for any ∈ , the following statements are equivalent.

Lemma 9.
A point * = ( * , ] * ) ∈ 2 is a solution of (20) if and only if * is a global optimal solution with the objective vanishing of (24).
In the following, we give the error bound for a polyhedral cone from [13] and error bound for a convex optimization from [14] to reach our aims.

Lemma 10. For polyhedral cone
Lemma 11. Let be a convex polyhedron in , and let be a convex quadratic function defined on . Let be the nonempty set of globally optimal solutions of the programming: with being the optimal value of on . There exists a scalar Before proceeding, we present the following definition introduced in [15].
By Lemma 8, ( ) is a convex function and the feasible set Ω is a polyhedral. Combining this with Lemmas 10 and 11, we immediately obtain the following conclusion.

Remark 14. It is clear that if is -strongly -monotone and
is strongly nonexpanding, then Moreover, the conditions which both and are Hölder continuous (or both and are Lipschitz continuous) in Theorem 13 are removed. Thus, Theorem 13 is stronger than Theorem 2.5 in [5]. Furthermore, by Theorem 2.1 in [5], the GNCP can be reformulated as general variational inequalities problem, and the conditions in Theorem 13 are also weaker than those in Theorem 3.1 in [15], Theorem 3.1 in [11], Theorem 3.1 in [10], and Theorem 2 in [9], respectively.
On the other hand, the condition that is -stronglymonotone and is strongly nonexpanding in Theorem 13 is extended compared with the condition that is strongly monotone with respect to (i.e., = 1) in Theorems 3.4 and 3.6 in [15], and it is also extended than compared with the condition is strongly monotone with respect to (i.e., = 1) in Theorem 3.1 in [11], and compared with the condition that ( ) = , ( ) is strongly monotone (i.e., = 1) in Theorem 3.1 in [10].
Using the following Definition 15 developed from the complementarity conditions in (22), we can further detect the error bound of the GNCP.

Proof.
Since by Assumption 3(A1), for any ∈ Ω, we have that is, Abstract and Applied Analysis To prove the assertion, we only need to show that the solution set Ω * is equal to the set For anỹ∈ Ω * , combining Lemma 9 with (20) yields that Letting =̃in (36) yields that Sincẽ, 0 ∈ Ω, using the similar technique to that of (21), we can obtain Since , 0 ∈ Ω, using the similar arguments to that of (21), one has Combining this with (41) yields that From (32), we deduce that Thus, using (21), one has Hence, ∈ Ω * .
Based on Lemma 16, we obtain the following conclusion.
In the following, we give an error bound of the Hölderian type [14].
Based on (18) and (21), the GNCP can also be written as From Lemma 19, we can establish the following global error bound for GNCP.

Theorem 20. Suppose that the hypotheses of Theorem 13 hold, and there exists point̂∈
such that Then, there exists constant 3 > 0 such that Proof.
. By Lemma 8, we have ( ) is a convex quadratic function. Combining this with (51), using Lemma 19 with = 0, this yields the following result where 6 is a positive constant. Obviously, 1 is a closed convex set. Thus, for any ∈ 2 , there exists a vector ∈ 1 such that For convenience, we also let From (50), we have Ω * = Ω ⋂ 1 , where Ω is defined in (24), so for any ∈ 1 , combining Lemma 10, one has where 7 is a positive constant, and the second and third inequalities follow from the fact that where the second equality follows from the fact that and the first inequality is by nonexpanding property of projection operator. Thus, Combining (56) with (59), for any ∈ 2 , we have where the second inequality follows from (56) with constant = 7 √ 2 + 2 + 2 , the third inequality uses (59), the fifth inequality follows from (53), the sixth inequality follows from the fact that the seventh and ninth inequalities follow from the fact that and the last inequality follows by letting 8 = max{√ , √ , √ , 1}.

Remark 21.
When is strongly monotone with respect to , that is, = 1, without the requirement of nondegenerate solution, the square root term in the error bound estimation is removed as stated in Theorem 20. Hence, the error estimation becomes more practical than that in Theorem 4.1 in [4].
Abstract and Applied Analysis 9

Global Error Bound for the GLCP
In this section, we consider the linear case of the GCP such that mappings and are both linear; that is, ( ) = + , ( ) = + with , ∈ × , , ∈ : where For problem (64), combining (18) with (23) and using a similar discussion in Lemmas 8 and 9, we also have the following conclusion.

Lemma 23. * ∈
is a solution of the GLCP if and only if * is global optimal solution with the objective vanishing of (64).
Based on (64), using the argument similar to that of Theorem 13, we can obtain the following conclusion.

Theorem 24. Under Assumptions 3(A1) and 3(A2), and that mappings and are both linear, there exists constant
Proof. For any ∈ , a direct computation yields that dist ( , * ) where the first inequality follows from Lemma 11 with constant 9 > 0 and Lemma 23, and the second inequality uses Lemma 10 with constant 10 > 0. By (67) and letting 4 = 9 max{ 10 , 1}, the desired result follows.
The following result further estimates the error bound for the GLCP.

Theorem 26. Suppose that the hypotheses of Theorem 24 hold, and the GLCP has a nondegenerate solution. Then, there exists
where 0 is a nondegenerate solution of GLCP, and is defined in (64). For any ∈ , a direct computation yields that where the first inequality follows from Lemma 10 with constant 11 > 0, and the second inequality uses (36). Letting Proof.
. By Lemma 22, ( ) is a convex quadratic function, and 2 is a closed convex set. For any ∈ , there exists a vector ∈ 2 such that Combining (51) and applying Lemma 19 yield the following result: where 12 is a positive constant. For convenience, we let From (50), we have * = ⋂ 2 , where is defined in (64). So for any ∈ 2 , combining Lemma 10 and using the similar technique to that of (56), one has where 13 is a positive constant. Using the fact that and using the similar technique to that of (57), one has where the second inequality is by nonexpanding property of projection operator. Thus, Combining (75) with (78), we know that for any ∈ , it holds that where the second inequalities follows from (75) with constant 1 = 13 √ 2 + 2 + 2 , the third inequality follows from (78), the fifth inequality follows from (73), the sixth inequality follows by letting 1 = max{ 1 , [ 1 (‖ ‖ + 2‖ ‖ + ‖ ‖ + 2‖ ‖) + 1] 12 }, and the seventh and ninth inequality follow from the fact that By (79) and letting 6 = 1 max{√ , √ , √ , 1}, the desired result follows.
Remark 29. In Theorem 28, without the requirement of nondegenerate solution, the square root term in the error bound estimation is removed. Hence, the error estimation becomes more practical than that in Theorem 4.1 in [4].

Comparison with Existing Error Bound
In the end of this paper, we will present an example to compare Theorem 13 and Theorem 2.5 in [5]. Furthermore, we will present two examples to show the conclusion in Theorem 13 can provide a global error bound for the GNCP, while the conclusion in Theorem 2.5 in [5] cannot do.

Example 31. For mappings ,
: + → involved in problem (81), we set It is easy to see that the solution set * = {0}, and one has where the first inequality follows from the fact that + +2 ≥ √ + √ . In fact, we consider the following four cases.
Example 32. For mappings , : → involved in problem (81), we set It is easy to see that the solution set * = {0}. Without loss of generality, we let > , and one has where the inequality follows from the fact that In fact, we consider the following four cases. Case 3 ( ≥ 0, < 0 and + ≥ 0). Then, and the desired result follows.

Conclusion
In this paper, we established some global error bounds on the generalized nonlinear complementarity problems over a polyhedral cone, which improves the result obtained for variational inequalities and the GNCP [4,5,[9][10][11] by weakening the assumptions. Surely, under milder conditions, we may establish global error bounds for GNCP and use the error bounds estimation to establish quick convergence rate of the methods for the GNCP. This is a topic for future research.

Introduction
Due to the large number of practical applications in various fields such as marketing, urban planning, supply chain management, and transportation, the continuous facility location problem has aroused the attention of many researchers ever since the pioneering work [1,2]. For a comprehensive review on this topic, see, for example, [3,4]. More specifically, the continuous facility location problem in a space is to seek the optimal locations for facilities serving customers (also called demand points), with certain objectives such as minimizing the sum of distances between facilities and customers.
The vast majority of literature treats the locations of facilities and customers as points in a space. Hence, the distances between facilities and customers are just the point-to-point distance without any ambiguity. In practical applications, however, regional demand arises frequently in such scenarios as uncertain demand, mobile demand, or cumbersome discrete situation whose number of demand points is extremely large. For such scenarios, many authors (e.g., [5][6][7][8][9][10][11]) promote that the regional customer, that is, the locations of customers are geometrically connected regions rather than points, should be considered. Therefore, in this paper, we consider the case of regional customer.
The question to be emphasized here, however, is how to measure the proximity from a regional customer to a facility. In the literature, different kinds of distances have been used to measure the required proximity. For example, the average distance evaluated by the proximity between the facility and some mean point in the interior of a regional customer (e.g., [10][11][12][13]) and the farthest distance measured by the proximity between the facility and the farthest point on the boundary of a regional customer (e.g., [8,9,14,15]). Definitely, the regional customers are treated by the average fashion when the average distance is considered, while the farthest distance realizes the worst-off nature in the sense that the regional customers are represented by their respective farthest points. In some real-life applications, however, the best-off nature is of importance in facility location; see, for example, [6,16,17]. To realize precisely the desired best-off nature, we need to consider the closest distance; that is, the proximity from a regional customer to a facility is evaluated by the closest distance to this facility. Thus, in this paper, we consider the closest distance to measure such a proximity. Note that the three kinds of distances have been well justified in [14,15] and the reader can be referred to them for more details.
Focusing on the real-life applications with vast eyes, the regional customers and the closest distances are highly essential to be considered. An example is about the communication network where several central servers are required to be sited. The demand points on the network are partially connected forming groups, each containing a large number of demand points. The points in the same group are connected to one another, and thus, each group becomes a regional customer in the plane. The servers need to be connected to the closest points in each regional customer, and the rest of the regional customer will be connected to the servers through that connection. Hence, the regional customers and the closest distances need to be used in the locations of these central servers.
As a matter of fact, in the literature of facility location, the distance is usually measured by norms such as -norms and block norms. References [3,18] discuss the approximations about the weighted -norms based on statistical evidence. References [19,20] investigate the use of block norms to obtain good approximation of actual distances. As it is known, the symmetry property of the norm assures that the distance from one point to another is always equal to the distance back. Nevertheless, as one of the first in-depth studies of mathematical location problems, [21] highlights the fact that in numerous real situations the symmetry of the distance is violated, for example, transportation in rush-hour traffic, flight in the presence of wind, navigation in the presence of currents, transportation on an inclined terrain, and so forth. For about two decades after the work of [21], however, no further research on this topic seems to be published, and only in the recent twenty years have the asymmetric distance problems started to attract the researchers' interest, for example, [22][23][24][25][26][27][28].
In this paper, we are interested in the locations of multiple facilities in the space with the aim of minimizing the sum of weighted distance between these facilities and regional customers, where the distance between a facility and a regional customer is evaluated by the closest distance. In addition, we formulate this problem in a quite general way with the aim of enhancing its practical applicability. First, we recognize that, usually, not all the points in the space are possible locations; that is, new facilities are often allowed to be sited within the confines of the restricted areas. Therefore, we introduce locational constraints so that both the unconstrained and constrained problems are taken into consideration in this paper. Second, since the distance in many practical situations is not necessarily symmetric, the gauge is used to measure the distance instead of the widely used norms. With the more general distance measuring function used, both the symmetric and asymmetric distances can be considered in our problem.
The rest of this paper is organized as follows. Section 2 states the formulation of our problem, which is shown to be nonconvex and NP-hard. The spirit of the well-known location-allocation heuristic algorithm, which consists of a location phase and an allocation phase in each iteration, is also discussed in this section. In Section 3, the subproblems arising in location phase and allocation phase are solved. More specifically, for the subproblem arising in allocation phase, the regional customers are allocated to facilities according to the nearest reclassification (NCR) heuristic, whereas for the single-source subproblem arising in location phase, the relationship between the subproblem and a monotone linear variational inequality (LVI) is firstly built, and then, a projection-contraction (PC) method is adopted to find the optimal location of the facility. A new locationallocation heuristic algorithm is proposed in Section 4 for solving our targeted problem, and its convergence is proved in Section 5. Preliminary numerical results are reported in Section 6 to verify the efficiency of the proposed algorithm. Finally, some conclusions are drawn in Section 7.

Model Description
This paper focuses on finding locations in the space for facilities, with the objective to minimize the sum of weighted closest distances between these facilities and regional customers. Note that the minisum single-source models with closest distance have been well justified in [6,16] where the distances are particularly measured by 1 -norm and 2 -norm, respectively. In our multi-source model, however, the gauge is used to measure the distances between facilities and regional customers, and thus, both the symmetric distance (including the 1 -norm used in [16] and the 2 -norm used in [6]) and the asymmetric distance are considered.
Let be a compact convex set in , and let the interior of contain the origin; then, the gauge of is defined by is called the unit ball of (⋅) due to This way to define a gauge from a compact convex set was first introduced by Minkowski [29]. The gauge (⋅) satisfies the following properties: (1) ( ) ≥ 0 and ( ) = 0 ⇔ = 0; (2) ( ) = ( ) for any ≥ 0; (3) ( + ) ≤ ( ) + ( ) for any and ∈ .
It follows from (2) and (3) that any gauge ( ) is a convex function of . The distance measuring function can be derived from a gauge by 3 When is symmetric around the origin, according to the definition of (1), we have Combining (3) and (4), it follows that which means that the distance measured by (⋅) is symmetric. On the contrary, when is not symmetric around the origin, (4) does not hold any more, and thus, the distance measured by the gauge (⋅) is asymmetric. Thus, when different compact convex sets are used as unit balls, different gauges (symmetric or asymmetric) can be generated and employed to measure distances in location problems, which depends on the requirements of practical applications.
Let Λ = { ⊂ | = 1, . . . , } denote the set of regional customers, and let ( = 1, . . . , ) be the location of the th facility. Each regional customer is simply assumed to be closed and convex. We denote the closest point in to the facility by Then, the closest distance between the facility and the regional customer , denoted by ( ), can be represented by When the gauge is used to measure distances, we have the following proposition for ( ) and ( ) which is similar to that in [16]. (6) is well defined, and the closest distance ( ) in (7) is a convex function of .

Proposition 1. Let be the location of a facility; then, the closest point ( ) in
Proof. Since is a convex set and (⋅) is a convex function, (6) is a convex problem, and thus, ( ) is well defined. Now, we prove that ( ) is a convex function of as follows. Let and be two points in and ∈ [0, 1]; then, due to ( ) and ( ) in and the convexity of , it follows that ( ) + (1 − ) ( ) ∈ , and thus, we have Therefore, ( ) is convex with respect to , and the proof is complete.
Based on the notations introduced above, now the constrained multi-source location problem (abbreviated as CMLP) we consider in this paper takes the following formulation: where ≥ 0 is the given demand required by the th customer, = ( 1 , . . . , ) is the variable of the locations of facilities to be determined, = ( ) × is the undetermined variable of which denotes the unknown allocation from the th facility to the th customer, and Π is the locational constraint for the th facility which is assumed to be a convex and closed set in .
More explanations are required for our model (9). First, the locational constraint Π in (9) can also be chosen as , and if all Π ( = 1, . . . , ) are , the CMLP (9) is a unconstrained problem, and thus, both the constrained and unconstrained problems are considered in our model. Second, mark that the minisum models analyzed in [6,16] are two special cases of CMLP (9), where = 1, Π = , and ( ) is particularly 1 -norm in [16] and Euclidean norm in [6].
It is noted that, with the presupposition that each facility is capable of providing sufficient services for the targeted customers, each customer is ultimately served only by the nearest facility in order to minimize the total sum of weighted distances. Therefore, the mathematical model CMLP also has the following form: where Π is the cartesian product of locational constraints; that is, When all are points not regions and all Π are , CMLP (9) reduces to the well-known multi-source Weber problem (MWP) which has wide applications in operations research, marketing, urban planning, and so forth; see, for example, [3,30,31]. Recall the fact that the multi-source Weber problem is nonconvex [32] and NP-hard [33], and therefore, heuristics algorithms are extremely popular and highly appreciated for overcoming the difficulty caused by its nonconvexity and NP-hardness; see, for example, [3,30,[34][35][36][37]. In particular, the classical location-allocation heuristic, also called the Cooper algorithm, has received much attention ever since it was presented originally by Cooper in [34] for MWP, whose attractive characteristic is that each iteration consists of a location phase and an allocation phase. Now, as a more general problem of MWP, CMLP is also nonconvex and NP-hard. Hence, in this paper, we are interested in applying the location-allocation heuristic to solve the CMLP (9). Accordingly, some location subproblems and allocation subproblems occur. To clarify it, let M = {1, 2, . . . , }, N = {1, 2, . . . , }, and then , and each Λ −1 ( = 1, . . . , ) in the partition is called one cluster. Then, at the th iteration, the location phase finds the candidates of locations of facilities (denoted by 1 , 2 , . . . , ) by solving the following constrained singlesource location problems (CSLP) with the closest distance under gauge for each cluster Λ −1 , = 1, . . . , : After the location phase, the allocation phase then revises the current partition of Λ to generate a new disjoint partition of Λ = {Λ 1 , Λ 2 , . . . , Λ } by the following nearest center reclassification (NCR) heuristic (see [38]): for some customer ∈ Λ −1 ℎ ( ∈ {1, 2, . . . , } and ℎ ∈ {1, 2, . . . , }), if ( ̸ = ℎ) is the nearest point for among all computed by (11), then (11) is the nearest facility for each regional customer in Λ −1 for any ∈ {1, 2, . . . , }, then ( = 1, 2, . . . , ) are the desirable locations of facilities and stop. Otherwise, we set = + 1 and repeat the iterations.

The Subproblems in Location and Allocation Phases
In this section, we will discuss the subproblems arising in the location phase and allocation phase. The allocation phase will partition the customers to clusters by the nearest center reclassification (NCR) heuristic, and the location phase will find the optimal location for each cluster by solving CSLPs (11).

Nearest Center Reclassification for Allocation of Customers.
The implementation of NCR heuristic to allocate regional customers can be executed by the following framework; see [38] for more details about this heuristic.
Algorithm 2 (the implementation of NCR). Given an initial Step 1. Set = 0 ( stores the number of reassignments); Step 2. Compute the facility of Λ −1 by solving CSLP (11), for = 1, 2, . . . , ; Step 3. For = 1, 2, . . . , do: Step 4. If = 0, then the iteration terminates with { 1 , . . . , } being the desirable locations for facilities and the customers in Λ −1 being served by ( = 1, . . . , ). (11). According to the spirit of location-allocation heuristic algorithm, our central task for the CMLP is to solve CSLP (11) in location phase by an efficient means. Recall that the CSLP is a generalized problem of the minisum models discussed in [6,16], where Π = and the gauge (⋅) are the particular 1norm in [16] and 2 -norm in [6]. For the model under 1 -norm in [16], by taking advantage of the piecewise linearity of the objective function, this model can be reduced to a standard minisum problem which can be easily solved by obtaining a median point for each coordinate separately. For the minisum model under 2 -norm in [6], an efficient Weiszfeld-type method is proposed, and the convergence of this method is analyzed. Similar to Weiszfeld procedure [2], one problem of the proposed method is that the singular case, that is, the current iterate happens to be within some location of regional customers, may occur during its implementation. Due to the use of the gradient of objective function in the iteration, this method will terminate unexpectedly once the singular case occurs. In order to tackle the undesirable singular case and make the Weiszfeld-type method computational effective, the authors suggest to ignore the gradient of ‖ − ( )‖ if ∈ and then add an extra descent check and a boundary check to the iteration. As pointed out by Theorem 1 in [6], however, the sequence generated by the proposed Weiszfeld-type method is possible to be convergent to a nonoptimal point which is on the boundary of the regional customer.

The Variational Inequality Approach for CSLP
In this section, a variational inequality approach is proposed to solve the general CSLP (11), where the locational constraints are imposed to the facility and the gauge is used as distance measuring function. Note that the study of variational inequality has received much attention due to its various applications arising in engineering, operations research, economics, transportation, and so forth; see, for example, [39][40][41][42][43][44][45]. Specifically, the CSLP (11) considered in this paper is first reformulated into an equivalent linear variational inequality (LVI), and then, a projection-contraction (PC) method is adopted to solve the LVI. Consequently, a sequence will be generated by the variational inequality approach, which is shown to be convergent to the optimal location of the facility of (11) even in the singular case. In addition, the closest points to the facility and the dual vectors with respect to the gauge can also be obtained from the generated sequence.
For convenience and succinctness, with the assumption that Λ −1 contains customers, throughout this section, we ignore some superscripts and subscripts in (11) and consider the simplified model of (11) without confusion: According to Proposition 1, it follows that CSLP (11), or equivalently MCSLP (12), is convex problem of .

LVI Reformulation of MCSLP.
For the gauge (⋅) in (12) which is defined by (1), there exists a dual gauge, (⋅), defined by Let be the unit ball of the dual gauge (⋅), which is also convex and compact and exactly the polar set of . The dual gauge of (⋅) is again (⋅), that is, which can also be rewritten as For more details about gauge and dual gauge, as well as their unit balls, the readers can be referred to [27]. According to (15), MCSLP (12) is equivalent to the following min-max problem: where each is a vector in = { ∈ | ( ) ≤ }. Since ( ) is the closest point to in , we can introduce ∈ to replace ( ). Hence, (16) is equivalent to Denote = ( 1 , . . . , ) , =( 1 , . . . , ) , and let ( * , * , * ) ∈ Π× ×B be the solution of (17); then, it follows that ( * , * , * ) is the saddle point of the objective function ∑ =1 ( − ); that is, Thus, ( * , * , * ) is the solution of the following linear variational inequality: A compact form of (20) is where = ( , , ) , Ω = Π × × B , Note that in (22) is a skew-symmetric matrix, then it is positive semidefinite, and thus, the linear variational inequality (21)-(22) is monotone.
Note that ( * , * , * ) is the saddle point of ( , , ) if and only if ( * , * , * ) ∈ Ω and which implies that On the other hand, let be any vector in B ; then, we have We choose in (25) Combining (24) and (26), it follows that all terms in (24) are equal, and therefore, which implies that ( * , * , * ) is the solution of (17).

Remark 4.
It is worth pointing out that the equivalence between the MCSLP (12) and the linear variational inequality (21)- (22) can also be obtained by the duality theory and the variable and the set ( = 1, . . . , ) in (16) are, respectively, the dual vector and dual ball in the space which satisfy ∈ .
The norms especially 1 , 2 , and ∞ are frequently used to measure distances in the literature; see, for example, [18,19]. It should be noted that the gauge used in this paper is an extension of norms which include 1 , 2 , and ∞ . When the gauge (⋅) is chosen as the 1 , 2 , and ∞ -norm, the dual gauge (⋅) will be the ∞ , 2 , and 1 -norm, respectively. Let where ‖ ⋅ ‖ 2 is the Euclidean norm; then, as a particular case of our single-source location problem (11), the problem under 2 -norm analyzed in [6] can be reformulated into the LVI (21)- (22) in which B is equal to B 2 and Π = .

A Projection-Contraction Method for LVI (21)-(22).
Among numerous effective numerical algorithms for solving VI, especially LVI, one famous one is the projectioncontraction (PC) method which was originally proposed by Uzawa [46]. The attractive characteristics of the PC method, for example, simpleness and effectiveness, have motivated further development on VI especially in computational aspects; see, for example, [39,[47][48][49]. In this section, we will summarize some concepts and results about linear variational inequalities and then adopt the projection-contraction method in [48] for solving LVI (21)- (22). More details about the proposed PC method can be referred to [48]. Let be a nonempty closed convex set of Q . For a given V ∈ Q , the projection of V onto denoted by (V) is the unique solution of the following problem: A basic proposition of the projection mapping on a closed convex set is It is well known (see, e.g., [50]) that for any > 0, * is the solution of LVI(Ω, , ) if and only if In and Ω * be the set of solutions of LVI(Ω, , ); then, for the positive semidefinite (not necessarily symmetric) matrix , the following theorem can be obtained.
Step 2. Calculate ( ) and set ( ) as Step 3. Calculate +1 as Step 4. Adjust as follows and set Let = + 1 and go to Step 1.

Remark 7. In
Step 4 of Algorithm 6, the parameter is selfadaptive during the iterations according to the value of ( ).
Note that is skew-symmetric, and thus, ( ) can also be rewritten as which is shown to be in [0, 1]. It follows from (39)-(40) that the two terms ‖ ( )‖ 2 and ‖ ( )‖ 2 in the denominator of ( ) are balanced by the self-adaptive parameter .

A Location-Allocation Heuristic
for CMLP (9) Recall the fact that for the well-known multi-source Weber problem (MWP), heuristics algorithms are extremely popular and frequently used for overcoming its nonconvexity and NP-hardness. In particular, the location-allocation heuristic algorithm has drawn much attention ever since its presentation by Cooper [34]. Note that the targeted CMLP (9) is an extension of the MWP and it is harder than MWP, and thus, in this paper, we also focus on applying the locationallocation heuristic algorithm for solving the CMLP in the spirit of Cooper's work.
Our previous analysis indicates that each iteration of the location-allocation heuristic algorithm to be presented consists of an allocation phase and a location phase. The allocation task generates a new disjoint partition of all the regional customers according to the principle of NCR as in the Cooper algorithm, and the location phase identifies the optimal locations for the current partition of customers via implementing the variational inequality approach for solving CSLPs.
Mark that the CMLP (9) differs from MWP mainly in that the customers are represented by regions rather than points. Consequently, the CSLPs involved in the location phase are constrained location problems with regional demand and closest distances under gauge. No doubt that the numerical implementation of the heuristic algorithm to be presented is expected to be more complicated than the location-allocation algorithms for MWP. Therefore, how to accelerate the convergence of the proposed heuristic deserves further consideration. To achieve this objective, we here consider a particular strategy for the initial partition of regional customers or the initial locations of facilities. In practical implementation, we suggest to choose the solution of the following constrained multi-source Weber problem (CMWP) as the initial locations of facilities for CMLP: CMWP: min where 's are geometric centers of the regional customers. Then, we apply the NCR to determine an initial partition of regional customers according to the solution of (42). For solving the constrained multi-source Weber problem (42), we 8 Abstract and Applied Analysis employ the location-allocation heuristic algorithm in [35]. As we will show by numerical experiments, this initialization strategy can accelerate the convergence of the proposed algorithm greatly. In the spirit of Cooper's work, the new heuristic algorithm is ready to be presented for solving the targeted CMLP (9), and its iterative framework can be elaborated as follows.
Step 2 (allocation phase). Update the partition of regional customers Λ from Λ −1 to Λ based on the spirit of NCR heuristic.
Step 3 If ‖ − −1 ‖ < , the current locations and partition are heuristic locations of facilities and heuristic partition of customers. Otherwise, set = + 1 and go to Step 1.

Remark 10.
At the ( + 1)th iteration, it is recommended to use (and the corresponding and ) in Step 1 as the initial iterate in the variational inequality approach for solving +1 ( = 1, . . . , ), considering the fact that Λ usually differs from Λ −1 slightly in practical implementation.
Remark 11. Compared to the main body of the proposed location-allocation heuristic algorithm, the workload of the initialization is relatively less. However, this initialization strategy can reduce its number of iterations and computing time, which will be verified by the numerical experiments to be reported in Section 6.2. Hence, the convergence of the proposed algorithm is accelerated greatly by this initialization strategy.
During the implementation of Algorithm 9, we denote the map L : Π × → Π × as the location operation in Step 1 and the map A : Π× → Π × as the allocation operation in the Step 2. It follows that where and +1 are, respectively, the locations of facilities in the th and ( + 1)th iteration, is the variable of the closest points to in the partition of Λ ,̃is the closest points to the new +1 in the partition of Λ , and +1 is the closest points to the new +1 in the new partition of Λ +1 . Then, the iterate scheme of the location-allocation heuristic is Let ( 0 , 0 ) denote the iterative sequence generated by the location-allocation heuristic for CMLP with the initial iterate ( 0 , 0 ). During the implementation of the proposed heuristic algorithm, we choose the initial iterate in location phase for solving CSLP as Remark 10 indicates. We first give the following proposition which reveals the monotonicity of the generated sequence ( 0 , 0 ).
Proof. Since +1 ̸ = , there exists at least one ∈ M such that +1 ̸ = . For such 's according to the following convex optimization problem we have Based on Remark 10, we know that if is the solution of (46), then +1 will be equal to . Therefore, +1 ̸ = implies that is not the solution of (46), and thus, +1 ( +1 ) < ( ). It follows that On the other hand, based on the principle of NCR in the allocation phase of Algorithm 9, we also have Abstract and Applied Analysis 9 By Combining (48) and (49), it follows that The proof is complete.
Based on the monotonicity of the generated sequence ( 0 , 0 ), the following theorem can be proved. Proof. After a finite number of iterations, if +1 = , then the iterates after will be constant, and thus, { } is convergent to ( , ) ∈ Π × , and the two assertions are both true.
Below, we will discuss the case that +1 ̸ = for any ∈ . First, we prove the first assertion. Since ∈ Π × and Π × is a compact space, it follows from the Bolzano-Weierstrass theorem that there exists a subsequence of { } which, say { } , converges to an element ∈ Π × ; that is, Note that ( , ) is a continuous function according to (43); then, Due to that { ( )} is a monotone sequence (Proposition 12) and has lower bound, then { ( )} is convergent. Thus, any subsequence of { ( )} will be convergent to the same value. Note that { ({ } )} is a subsequence of { ( )} and it is convergent to ( ); then, it follows that The second assertion can easily be proved. Let be an accumulation point of { }; then, there exists one subsequence { } which converges to , and due to the continuity of ( , ), we have The first assertion has shown that { ( )} is convergent, and note that { ({ } )} is a subsequence of { ( )}; then, it follows that Thus, all accumulation points of { } have the same objective functional values equal to lim → ∞ ( ).
Proof. Note that the CSLP (11) is a convex problem, then are continuous. Since Π × is a compact space and also a Hausdorff space and every continuous map from a compact space to a Hausdorff space is closed, it follows that L is closed over Π × .
Lemma 15. Let 0 be a given vector in Π × and Δ : Proof. It is known that every closed subset of a compact space is also compact, and therefore, it is enough to prove that Δ is a closed set.
Due to the continuity of , it follows that On the other hand, { } ∈ Δ implies By combining (58), (59), and the continuity of , the following inequality is obtained: and accordingly, ∈ Δ. This means that Δ is closed, and the proof is complete. Now, we are ready to prove the convergence of the proposed location-allocation heuristic (Algorithm 9). Let Ξ ⊆ Π × be the nonempty local solution set of CMLP (9). Recall that in the location-allocation Cooper algorithm for MWP, if +1 = occurs after a finite number of iterations, the iterates after will be constant. Then, no further improvement is possible for MWP, and it follows from [32,34] that the is a local solution of MWP. Similarly, in the proposed location-allocation heuristic algorithm for CMLP, if +1 = occurs, the iterates after = ( , ) will also be constant. Then, exactly as in the location-allocation Cooper algorithm for MWP, no further improvement is possible for CMWP, and a local solution, namely, , is obtained. Hence in this case, { } is convergent to the ∈ Ξ. However, it is not assured that +1 = always occurs during the implementation of Algorithm 9, and therefore, we assume that +1 ̸ = for any ∈ and prove the convergence in this case.
So, it is enough to prove ∈ Ξ. We prove this by contradiction. Assume that ∉ Ξ; that is, is not a solution, and we consider the subsequence { +1 } . Denote Δ = { ∈ Π × | ( ) ≤ ( 0 )}. Due to Proposition 12, it follows that for all ∈ we have ( +1 ,̃) ∈ Δ and ( +1 , +1 ) ∈ Δ. According to Lemma 15, Δ is compact, and thus, there exists ⊂ such that According to Lemma 14, the map L is closed at ∈ Π × ; then, it follows that 1 = L( ). Further, due to ∉ Ξ, V 1 will be not equal to . Otherwise, we can choose as the initial iterate of the location-allocation heuristic algorithm, then, the sequence generated by the algorithm will be constant. It follows from the first case (i.e., +1 = ) that will be a local solution, which contradicts with ∉ Ξ. Therefore, we can obtain V 1 ̸ = . Together with the monotone proposition of L (48), we have the inequality On the other hand, note that A( +1 ,̃) = ( +1 , +1 ); then, by the monotonicity of A (49), it follows that and thus, by taking the limit for (64) and by the continuity of , we know ( 2 ) ≤ ( 1 ). Combining this with (63), we obtain However, note that 2 and are two accumulation points of { }, and according to the second assertion of Theorem 13, ( 2 ) = ( ), which will contradict with (65). Therefore, our assumption is wrong, and thus, ∈ Ξ.
As a result, we have the following convergence theorem for the sequence generated by the proposed locationallocation algorithm.
Theorem 17. The sequence ( 0 , 0 ) generated by the proposed location-allocation heuristic algorithm either converges to a point in Ξ or all accumulation points of ( 0 , 0 ) belong to Ξ.

Numerical Results
This section reports some preliminary numerical results to verify the theoretical assertions proved in previous sections. Section 3.1, reports some numerical results of the proposed variational inequality approach for the CSLP (11) (or equivalently (12)) which includes (1) the results of the comparison between our approach and the Weiszfeld-type method by solving the example in [6] and some randomly generated unconstrained examples under Euclidean distances and (2) the results of our approach for solving some randomly generated constrained examples under a gauge. These numerical results demonstrate the efficiency of the proposed variational inequality approach for CSLP. In the second subsection, we apply the proposed location-allocation heuristic algorithm to solve some randomly generated examples of the CMLP (9). In particular, the effectiveness of the initialization strategy adopted in this heuristic for accelerating convergence will be justified. All the programming codes are written by Matlab 2012b and were run on an ASUS notebook (Intel Core2 Duo T6670 2.20 GHz).

Numerical Results of Variational Inequality Approach for CSLP.
When applying the variational inequality approach for solving CSLP and MCSLP (12), theoretically, the initial iteration 0 in Algorithm 6 can be chosen arbitrarily in Ω. In practical implementation, however, we choose 0 judiciously similar to the initialization strategy in the location-allocation heuristic: let ( = 1, . . . , ) be the centers of regional customers, solve the following single-source Weber problem (SWP): by the projection-contraction method in [35], and then use its solution as the initial iterate for Algorithm 6. We call this the initialization strategy of variational inequality approach. In addition, throughout our experiments of VI approach, the 1 and 2 in Algorithm 6 are chosen as 1 and 0, respectively. We first solve the example given in [16] by the proposed variational inequality approach and the Weiszfeld-type method in [16].
In order to clarify the comparison of two methods, we choose the same stopping criterion as ‖ +1 − ‖ ≤ 10 −4 (throughout this section, ‖ ⋅ ‖ is the ∞ -norm). We test this example for 100 times with the same initial iterate for the two methods which is randomly generated in [0, 5] × [0, 3], and the numerical results including the location of new facility, the closest points to the facility, number of iterations, and computing time in units of second are reported in Table 1. According to Table 1, it follows that both methods can get the optimal location of the new facility. Though our variational inequality approach needs smaller iteration numbers and less computing time, both methods are efficient for solving Example 1 given in [6].
Recall that the sequence generated by the Weiszfeld-type method in [16] is possible to be convergent to a nonoptimal point on the boundary of the regional customer, as indicated in [16]. The variational inequality approach, however, can obtain the optimal location of the new facility, which is guaranteed by the theoretical analysis in Section 3.2. To illustrate this, we compare two methods by solving the following particular example.
Example 19. Similar to Example 18, the number of customers = 5, and all customers are unit squares whose edges are parallel to the axes. Let 1 be the distance between any two neighboring customers, and we set 1 equal to 1, 0.1, 0.01, and 0.001, respectively. The weights ( = 1, . . . , 5) are randomly generated in the area of regional customers.
Note that in Example 19 the parameter 1 reflects how close the customers are away from its neighborhoods. We test this example for a large number of times with the stopping criterion ‖ +1 − ‖ ≤ 10 −4 , and the average numerical results are reported in Table 2. In this table, each row reports the average results by testing Example 19 for one hundred times. The column of "No. of Iter. 0 " gives the average iteration times of initialization in the variational inequality approach, and the column of "No. of Iter. " reports the average iteration times of both methods. The two columns of "CPU" give the average computing time in units of second for variational inequality approach (including the computing time for initialization) and Weiszfeld-type method, respectively. The columns of "Obj. " give the average objective functional value obtained by the two methods, and the column of "Impro. Percent" gives the improvement percentage in objective functional values of the VI approach to Weiszfeld-type method. Remark that the convergence of VI approach to the optimal location of new facility is guaranteed, and then the column of "Freq. Num. " reports the frequency among one hundred times that the Weiszfeld-type method can get the same solution as VI approach; that is, it does not converge to the nonoptimal solution on the boundary of the customer.
According to Table 2, it follows that both the VI approach and the Weiszfeld-type method are efficient for solving this particular example, and both of them need a small number of iterations and little computing time. In comparison of two methods, we can find that VI approach needs more iteration times and computing time than Weiszfeld-type method. From the column of "Impro. Percent, " however, we can conclude that VI approach can obtain a better solution (in fact, the solution obtained by VI approach is the optimal location of new facility) than Weiszfeld-type method. In addition, according to the last column, we find that when 1 decreases, that is, the customers become closer and closer, the frequency that Weiszfeld-type method obtains the optimal solution gets smaller and smaller. When 1 is 0.001, which implies that the customers are quite close to the neighborhoods, this frequency is totally less than 20. In other words, for this particular example with 1 = 0.001, the sequence generated by Weiszfeld-type method has a great possibility (more than 80%) to be convergent to a nonoptimal solution which is on the boundary of the customer. On the contrary, when 1 is 1, which means that the customers are enough far away from one another, Weiszfeld-type method can obtain the optimal location of facility in most cases, and the frequency even exceeds 90.
Since our main effort in this paper is to solve the general location problem under gauge and locational constraint, it is necessary to apply the variational inequality approach to solve some CSLPs. In particular, we test a large number of randomly generated CSLPs with the number of customers from 10 to 2000. In the experiments, all regional customers are assumed to be square units, and their edges are parallel to the coordinate axes. The geometric centers of all regional customers are randomly generated in [−100, 100] 2 ; the weights of the regional customer are all randomly chosen in (1,5); the locational constraint is ‖ − ‖ ≤ , where the center is randomly generated in [−100, 100] 2 , and the radius is randomly generated in (1,5); the stopping criterion of VI approach is chosen as and the initial iterate is randomly generated in [−100, 100] 2 ; the gauge (⋅) is generated with the unit ball set as For each , we test one hundred randomly generated CSLPs, and the average numerical results are reported in Table 3.
To illustrate the effect of the initialization strategy of VI approach, we also report the results of VI approach without the initialization strategy. The columns of "VI approach with Initial. " and "VI approach without Initial., " respectively, report the average number of iterations and average computing time of variational inequality approach with and without the initialization strategy.  According to Table 3, it is easy to conclude that the variational inequality approach is effective for solving CSLP under gauge considering the difficulty of this problem. In addition, the number of iterations of "VI approach with Initial. " is less than that of "VI approach without Initial., " which shows that the initialization strategy can accelerate the convergence of the variational inequality approach. This strategy, however, does not necessarily reduce the computing time of VI approach, especially when is small, for example, = 10, 20, 50, and 100, which can be explained as follows. When the number of customers is small, the variational inequality problem is small scale, and thus, it can be solved in a short time. In this case, the computational workload of initialization plays an important role in the total workload, and therefore, the computing time of "VI approach with Initial. " is greater than that of "VI approach without Initial. " due to the computational iterations for initialization. With increasing, the scale of VI problem as well as the number of iterations becomes larger. Then, in the comparison of the workload of VI approach, the workload of initialization can almost be ignored, and therefore, the computing time of "VI approach with Initial. " will be smaller than that of "VI approach without Initial. " As a matter of fact, Table 3 reveals the computational necessity of the initialization strategy of variational inequality approach for large-scale CSLP; for example, the iteration number and computing time are reduced about 1/5 by the initialization strategy when = 2000.

Numerical
Results of Heuristic Algorithm for CMLP. This subsection applies the proposed location-allocation heuristic algorithm (Algorithm 9) to solve a large number of CMLPs  (9) and also verifies the necessity of the initialization strategy in Algorithm 9. In the experiments, we again generate a large number of CMLPs with unit square customers and assume that the edges of these regions are parallel to the coordinate axes. The geometric centers of all the demand regions are randomly generated in [−100 100] 2 , and all the demands, ( = 1, 2, . . . , ), are randomly generated in [ 1 10 ]. The gauge (⋅) is also defined with the unit ball set as (68). We test the scenario with = 100, 200, 500, and 1000 and = 2, 4, 6, 8, and 10; the locational constraints are ‖ − ‖ ≤ ( = 1, . . . , ), where the radius is randomly generated in [1,10], and the center is given in advance as follows: The initial locations of facilities are randomly generated in [−100 100] 2 , and the stopping criterion used in Algorithms 6 and 9 is chosen as To show the significance of the initialization strategy, we compare the numerical performance of the locationallocation heuristic algorithm with initialization strategy (denoted by "Algorithm 9 with Initial. ") and without this initialization strategy (denoted by "Algorithm 9 without Initial. "). In the initialization step of Algorithm 9, the location-allocation algorithm in [35] is adopted to solve the corresponding CMWP (42), where a PC method is proposed to solve the subproblems in location phase, and the numbers of iterations of the algorithm in [35] (denoted by "Iter 0 ") and the average iteration numbers of PC method in one iteration of the algorithm (denoted by "Iter 0 -PC") are reported. The columns of "Iter. " and "CPU, " respectively, report the number of iterations and computing time of Algorithm 9 with and without initialization strategy. Since the efficiency of Algorithm 9 is mainly determined by the number of iterations of the variational inequality approach, we also report the average number of iterations of Algorithm 6 in one iteration of Algorithm 9 (denoted by "Iter.-PC"). For each given pair ( , ), we test the CMLP for 100 times, and the computational performance is reported in Table 4.
It follows from Table 4 that the proposed locationallocation heuristic, with or without the initialization strategy, is capable of tackling the CMLP (9) efficiently, even for large-scale cases. Also, the necessity of the initialization strategy is evident. In fact, this strategy reduces both the number of iterations and the computing time by about 50%.
Another interesting fact obtained from Table 4 deserves further illustration. Recall that with the number of new facilities ( ) increasing, the number of subproblems (CSLPs) in location phase increases too. According to Table 4, however, we find that for fixed number of customers ( ), with increasing, the number of iterations for solving CSLPs in one iteration of Algorithm 9 does not increase but almost decreases with , especially for large-scale CMLP. This can be illustrated roughly as follows. For fixed , when increases, the average number of customers in each Λ becomes smaller, which implies that the scale of the involved CSLP (11) in location phase is smaller. According to Table 3, it follows that we need smaller number of iterations for small-scale CSLP, and thus, the total number of iterations for solving CSLPs decreases. Similarly, due to the same reason, the computing time for solving CMLP also decreases with increasing, as reported in the column of "CPU" in Table 4.

Conclusion
In this paper, we are interested in the locations of multiple facilities in the space with regional demands, where the closest distance is used to measure the proximities between facilities and customers. With locational constraints introduced for the locations of new facilities and with the gauge used as the distance measuring function, the problem considered in this paper has much more applications in practice. Due to its nonconvexity and NP-hardness, a new location-allocation heuristic algorithm is proposed to solve this problem, and its convergence is proved under mild assumptions. Some preliminary numerical experiments are reported to verify the computational efficiency of the proposed algorithm.

Introduction
In various fields of applied mathematics and engineering, many problems can be equivalently formulated as a separable convex optimization problem with two blocks; that is, given two closed convex functions : R → R ∪ {+∞}, = 1, 2, to find a solution pair ( * 1 , * 2 ) of the following problem: where is a matrix in R × , = 1, 2, and is a vector in R . The classical alternating direction method of multipliers (ADMM) [1,2] applied to problem (1) yields the following scheme: where is a Lagrangian multiplier and > 0 is a penalty parameter. Possibly due to its simplicity and effectiveness, the ADMM with two blocks has received continuous attention both in theoretical and application domains. We refer to [3][4][5][6][7][8] for theoretical results on ADMM with two blocks and [9][10][11][12][13] for its efficient applications in high-dimensional statistics, compressive sensing, finance, image processing, and engineering, to name just a few.
In this paper, we concentrate on the linearly constrained convex programming problem with three blocks: 2 Abstract and Applied Analysis where 3 : R 3 → R ∪ {+∞} is a closed convex function and 3 is a matrix in R × 3 . For solving (3), a nature idea is to extend the ADMM with two blocks to the ADMM with three blocks in which the next iteration ( +1 2 , +1 3 , +1 ) is updated by Similar to the ADMM with two blocks, the ADMM with three blocks has found numerous applications in a broad spectrum of areas, such as doubly nonnegative cone programming [14], high-dimensional statistics [15,16], imaging science [17], and engineering [18]. Even though its numerical efficiency is clearly seen from those applications, the theoretical treatment of ADMM with three blocks is challenging and the convergence of the ADMM is still open given only the convex assumptions of the objective function. To alleviate this difficulty, the authors of [19,20] proposed predictioncorrection type methods to solve the general separable convex programming; however, numerical results show that the direct ADMM outperforms its variants substantially. Therefore, it is of great significance to investigate the theoretical performance of the ADMM with three blocks even only to provide sufficient conditions to guarantee the convergence.
To the best of our knowledge, there exist only two works aiming to attack the convergence problem of the direct ADMM with three blocks. By using an error bound analysis method, Hong and Luo [21] proved the linear convergence of the ADMM with blocks for sufficiently small subject to some technical conditions. However, the sufficiently small requirement on makes the algorithm difficult to implement. In [22], Han and Yuan employed a contractive analysis method to establish the convergence of ADMM under the strongly convex assumptions of and the parameter less than a threshold depending on all the strongly convex moduli. In this paper, we firstly prove the convergence of ADMM with three blocks under two conditions weaker than those of [22]. In our conditions, the threshold on the parameter only relies on the strongly convex moduli of 2 and 3 , and furthermore 1 is not necessarily strongly convex in one of our conditions. Also, the restricted range of in this paper is shown to be at least three times as big as that of [22].
In order to accelerate the ADMM with three blocks, we also propose a relaxed ADMM with three blocks which involves an additional computation of optimal step size. Specifically, with the triple ( 2 , 3 , ), we first generate a predictor (̃2,̃3,̃) according to (5) and then obtain ( +1 2 , +1 3 , +1 ) in the next iteration by where ∈ (0, 2) and * is special step size defined in (43). The convergence of the relaxed ADMM is also established under mild conditions. We should mention that it is possible to modify the analyses given in this paper to be problems with more than three blocks of separability. But this is not the focus of this paper.
The remaining parts of this paper are organized as follows. In Section 2, we list some preliminaries on the strongly convex function, subdifferential, and the ADMM with three blocks. In Section 3, we first show the contractive property of the distance between the sequence generated by ADMM with three blocks and the solution set and then prove the convergence of ADMM under certain conditions. In Section 4, we extend the direct ADMM with three blocks to the relaxed ADMM with an optimal step size and establish its convergence under suitable conditions. We conclude our paper in Section 5.
Abstract and Applied Analysis 3

Preliminaries
Throughout this paper, we assume , = 1, 2, 3, are strongly convex functions with modulus ≥ 0; that is for each . Note that is a strongly convex function with modulus 0 being equivalent to the convexity of . Let be a point of dom( ); the subdifferential of at is defined by From Proposition 6 in [23], we know that, for each , is strongly monotone with modulus which means The next lemma introduced in [22] plays a key role in the convergence analysis of the ADMM and the relaxed ADMM with three blocks.

The ADMM with Three Blocks
In this section, we first investigate the contractive property of the distance between the sequence generated by ADMM with three blocks and the solution set under the condition that 0 < ≤ min{ 2 /‖ 2 ‖ 2 2 , 3 /‖ 3 ‖ 2 2 }.
By (14) and the monotonicity of 3 (⋅) (11), it is easily seen that Then for each , where the last "≥" follows from the elementary inequality Since by direct computations, we further obtain that ), implies Note that We complete the proof of this lemma.
Proof. By the inequality (13), it follows that the sequence Hence { 1 1 } is also bounded. Moreover, from (13) we see immediately that According to the condition that 0 < It therefore holds that Therefore, the sequence 3 3 , }, implies that { 2 , 3 , } is bounded, and { 1 } is bounded given the condition 1 > 0 or 1 is of full column rank. Moreover, since by the first equality in (25) and the third equality in (26), it holds that Abstract and Applied Analysis 5 We proceed to prove the convergence of ADMM by considering the following two cases.

The Relaxed ADMM with Three Blocks
For the ADMM with two blocks, Ye and Yuan [25] developed a variant of alternating direction method with an optimal step size. Numerical results demonstrated that an additional computation on the optimal step size would improve the efficiency of the new variant of ADMM. In this section, by adopting the essential idea of Ye and Yuan [25], we propose 6 Abstract and Applied Analysis a relaxed ADMM with three blocks to accelerate the ADMM via an optimal step size. For notational simplicity, we write With = ( 2 , 3 , ), the new iterate of extended ADMM is produced by wherẽis the solution of (5) and * is defined by * := Φ ( ,̃) Lemma 5. Let the sequence { } be generated by the relaxed ADMM with three blocks. Then, if 0 < ≤ 3 /‖ 3 ‖ 2 2 , the following statements are valid: Proof. By direct computations to Φ( ,̃), we know that where the second inequality follows Cauchy inequality. It therefore holds that which completes the proof of the first part. By Lemma 1 and the elementary inequality (17), it can be easily verified that and then This, together with the fact that * ≥ 1/6, completes the proof.
Based on the above inequality, we are able to prove the following convergence result of the relaxed ADMM with three blocks. Since the proof is in line with that of Theorem 3, we omit it.

Conclusion Remarks
In this paper, we take a step to investigate the ADMM for separable convex programming problems with three blocks. Based on the contractive analysis of the distance between the sequence and the solution set, we establish theoretical results to guarantee the global convergence of ADMM with three blocks under weaker conditions than those employed in [22]. By adopting the essential idea of [25], we also present a relaxed ADMM with an optimal step size to accelerate the ADMM and prove its convergence under mild assumptions.

Introduction
Let be a real Hilbert space with inner product ⟨⋅, ⋅⟩ and norm ‖ ⋅ ‖, let be a nonempty closed convex subset of , and let : → be an operator. Then the variational inequality problem VI( , ) [1] is to find * ∈ such that * ∈ , ⟨ − * , ( * )⟩ ≥ 0, ∀ ∈ . (1) The literature contains many methods for solving variational inequality problems; see  and references therein. According to the relationship between the variational inequality problems and a fixed point problem, we can obtain * is the solution of VI ( , ) where the projection operator is the projection from onto , that is, In this paper, : → is an operator with : -Lipschtz and -strongly monotone; that is, satisfies the following conditions: If is small enough, then is a contraction. Naturally, the convergence of Picard iterates generated by the right-hand side of (2) is obtained by Banach's fixed point theorem. Such a method is called the projection method or more results about the projection method see [6,8,20] and so forth.
In fact, the projection in the contraction methods may not be easy to compute, and a great effort is to compute the projection in each iteration. Yamada and Deutsch have provided a hybrid steepest-descent method for solving the VI( , ) [2,3] in order to reduce the difficulty and complexity of computing the projection . Subsequently, the convergence of hybrid steepest-descent methods was given out by Xu and Kim [4] and Zeng et al. [5]. Naturally, by analyzing several three-step iterative methods in each iteration by the fixed pointed equation, we can obtain the Noor iterations. Recently, Ding et al. [7] proposed a threestep relaxed hybrid steepest-descent method for variational inequalities, and the simple proof of three-step relaxed hybrid steepest-descent methods under different conditions was introduced by Yao et al. [24]. The literature [14,16] described a modified and relaxed hybrid steepest-descent (MRHSD) method and the different convergence of the MRHSD method under the different conditions. A set of practical numerical experiments in the literature [16] demonstrated that the MRHSD method has different efficiency under different conditions. Subsequently, the prediction-correction and relaxed hybrid steepest-descent method (PRH method) [15] makes more use of the history information and less decreases the loss of information than the methods [7,14]. The PRH method introduced more descent directions than the MRHSD method [14,16], and computing these descent directions only needs the history information.
In this paper, we will prove the strong convergence of PRH method under different and suitable restrictions imposed on parameters (Condition 12), which differs from that of [15]. Moreover, the proof of strong convergence is different from the previous proof in [15], which is not similar to that in [7] in Step 2. And more importantly, numerical experiments verify that the PRH method under Condition 12 is more efficient than that under Condition 10, and the PRH method under some descent directions is more slightly efficient than that of the MRHSD method [14,16]. Furthermore, it is easy to obtain these descent directions.
The remainder of the paper is organized as follows. In Section 2, we review several lemmas and preliminaries. We prove the convergence theorem under Condition 12 in Section 3. In Section 4, we give out a series of numerical experiments, which demonstrated that the PRH method under Condition 12 is more efficient than under Condition 10. Section 5 concludes the paper.

Preliminaries
In order to proof the later convergence theorem, we introduce several lemmas and the main results in the following.

Lemma 1. In a real Hilbert space H, there holds the inequality
The lemma is a basic result of a Hilbert space with the inner product.

Lemma 2 (demiclosedness principle). Assume that is a nonexpansive self-mapping on a nonempty closed convex subset of a Hilbert space . If has a fixed point, then ( − ) is demiclosed. That is, whenever is a sequence in weakly
converging to some ∈ and the sequence ( − ) strongly converges to some ∈ , it follows that ( − ) = . Here is the identity operator of .
The following lemma is an immediate result of a projection mapping onto a closed convex subset of a Hilbert space.

Convergence Theorem
Before analyzing the convergence theorem, we first review the PRH method and related results [15].
We obtain the strong convergence theorem of PRH method for variational inequalities under different assumptions. Proof. We divide the proof into several steps.
6 Abstract and Applied Analysis By (32), we immediately obtain By a series of computations, we can get Hence, by (28), (33), and (34), we also obtain Using Steps 2 and 3, it is easy to obtain the following corollary.
By * being the unique solution of VI( , ), we can obtain Since ‖̂− ‖ ≤ ‖̂− ‖ → 0, we immediately conclude that Step 5. By Step 1 and Lemma 1, we have Abstract and Applied Analysis and 0 ≪ < ∞. Denote We can rewrite (42) as In fact, , satisfies Lemma 5; according to we obtain Moreover, by Step 4, we also obtain Furthermore, by (43), (47), and (48), it is easy to obtain Consequently apply Lemma 5 to obtain

Numerical Experiments
The problem considered in this section is where ‖ ⋅ ‖ is the matrix Fröbenis norm; that is, Note that the matrix Fröbenis norm is induced by the inner product The problems arise from finance and statistics, and we form the test problems similarly as in [9,21].
Let , be given × symmetric matrices, and asymmetric which differs from previous approaches [9,21], and it is to be noted that the extended contraction method (EC method) [9] has much difficulty in computing the examples when is asymmetric, where ≤ in element wise:
Then (51) is equivalent to the following variational inequality: So we get According to Condition 10, we take the following parameter sequences, and let Condition 10 denote the parameter sequences: According to Condition 12, we take the following parameter sequences, and let Condition 12 denote the parameter sequences: Obviously, we have much difficulty in computing the projection of [ ], for all ∈ . In order to reduce the difficulty and complexity of computing the projection , we define by where ( ) = min ( , max ( , )) , which can be computed without difficulty and the fixed point set of Fix( ) = . According to Theorems 11 and 13, the sequences generated by Algorithm 8 under Conditions 10 and 12 are convergent. The computation begins with ones ( , ) in MATLAB and stops as soon as ‖ +1 − ‖ ≤ 10 −3 or 10 −2 . All codes were implemented in MATLAB 7.1 and ran at a Pentium R 1.70G processor, 2G Acer note computer.
We test the problems with = 100, 200, 300, 400, 500, 1000, and 2000. The test results with the PRH method under   Tables 1, 2, 3, and 4. And the CPU time is in seconds. It is to be noted that the results of extended contraction method are only given out when the iteration step (It) is less than or equal to 100.
Test Examples 1. In this example we generate the data in a similar manner as in [9]. The entries of diagonal elements of are randomly generated in the interval (0, 2); the entries of off-diagonal elements of are randomly generated in the interval (−1, 1) (Algorithm 1): When ≥ 1000 and tolerance 10 −4 , the computation time of the proposed method is too long, so the results of the PRH method give out approximate solution with ≥ 1000 and tolerance 10 −3 in the following. And the extended contraction method (EC method) has much difficulty in computing the examples when is asymmetric. Furthermore, by introducing auxiliary variable, the certain projection method or relaxed-PPA method [10] can be implemented by these tests.
From Tables 1 to 3, we found that the iteration numbers and CPU time of PRH under Condition 12 are more efficient than that under Condition 10. In Table 4 of our method, the tests' results give out that the PRH method under some descent directions is more slightly efficient than those of the MRHSD method [14,16], and it is easy to obtain these descent directions. Furthermore, it is important to find , , and by Tables 2 and 4.

Conclusions
We have proved the strong convergence of PRH method under Condition 12, which differs from Condition 10. The result can be considered as an improvement and refinement of the previous results [14]. And more importantly, numerical experiments demonstrated that the PRH method under Condition 12 is more efficient than that under Condition 10, and the PRH method under some descent directions is more slightly efficient than that of the MRHSD method. How to select parameters of the PRH method for solving variational inequalities is worthy of further investigations in the future.

Introduction
The variational inequality (VI (Ω, )) in the finite-dimensional space is to determine a vector ∈ Ω such that where Ω ∈ R is a nonempty closed convex subset and is a continuous mapping from R 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: where ∈ R × 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-knownextragradient type methods (EGM) as special cases. Assume that is Lipschitz continuous; that is, there is ∈ (0, 1), such that Then at each iteration EGM takes the following general form: In this paper, we consider the following variational inequalities: find a vector ∈ D such that with := ( ) , ( ) := ( ( ) ( ) ) , where D ∈ R + is a nonempty closed convex subset and : X → R and : Y → R are monotone operators. Problem (5) is referred to as a structured variational inequality (SVI) [6].

Abstract and Applied Analysis
By attaching a Lagrange multiplier vector ∈ R to the linear constraints + = , the VI problem (5) is converted into the following form: where The compact form is For the purpose of parallel computing, the proximal alternating directions method (PADM) generates̃= (̃,̃, ) ∈ Ω as follows [7,8]: first find añ∈ X such that Then find añ∈ Y such that Finally, updatẽviã 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][9][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]) where > 0 is a chosen step size. The PADM (11)-(13) is often easy to implement under the assumption that the decomposed subproblems have closedform 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/ ) ) −1 ( ) and ( + (1/ ) ) −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 ( + ) −1 ( ). 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-quadraticproximal-(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.

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) X, Y 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, (A2) The mapping is point-to-point, monotone, and continuous.
The following well-known properties of the projection operator will be used in the coming analysis.

Lemma 1.
Let Ω ∈ R be a nonempty closed convex set; let Ω (⋅) be the projection operator onto Ω under the -norm. Then For any arbitrary positive scalar and ∈ Ω, let ( , ) denote the residual function associated with the mapping ; that is,

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: Obviously, is a symmetric positive definite matrix whenever > 0, > 0, and > 0, and we also have = .

Description of the Algorithm
Algorithm 4. It is a prediction-correction-based algorithm for the SVI (Ω, ).
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.

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 (Ω, ). (22)- (24), and let the matrix be given in (20). Then one has
The following lemma plays a key role in proving the convergence of the algorithm.
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, Abstract and Applied Analysis 5 Then one has where Proof. One can see that On the other hand, since = , using the monotonicity of and Lemma 9, we have Combining (42)-(43) together, we have Thus, ( ) is a lower bound of ( ) for any > 0.
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 Proof. According to Lemma 11, Moreover, it follows from (26) that the step size 6 Abstract and Applied Analysis Based on the conditions (33), we have Substituting (50) into (47), we have 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 contractiontype methods. It follows from (46) that { } is bounded, and we have that Consequently, since (see (22) and (23) Becausẽ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 Consequently, 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 Thus, the sequence { } converges to ∞ , which is a solution of SVI (Ω, ).

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 MAT-LAB Version 2009b. We consider the following optimization problem: where ∈ R × , ∈ R × are symmetric positive semidefinite matrixes, ∈ R × , ∈ R × , and ∈ R .
Using the KKT condition, the problem (59) can be converted into the following variational inequality: find * = ( * , * , * ) ∈ R + + such that 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: (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: 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: 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 ) = (0, 0, 0), and the stopping criterion is 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.
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.

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 8 Abstract and Applied Analysis 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.
For the special case of (1) with = 2, the problem has been studied extensively. Among lots of numerical methods, one of the most popular methods is the alternating direction method of multipliers (ADMM) which was presented originally in [1,2]. The iterative scheme of ADMM for (2) is as follows: where is Lagrange multiplier associated with the linear constraints and > 0 is the penalty parameter. The convergence of ADMM for (2) was also established under the condition that the involved functions are convex and the constrained sets are convex too.
While there are diversified applications whose objective function is separable into ≥ 3 individual convex functions without coupled variables, such as traffic problems, the problem of recovering the low-rank, sparse components of matrices from incomplete and noisy observation in [3], the constrained total-variation image restoration and reconstruction problem in [4,5], and the minimal surface PDE problem in [6], it is thus natural to extend ADMM from 2 blocks to blocks, resulting in the iterative scheme: . . .
Unfortunately, the convergence of the natural extension is still open under convex assumption, and the recent convergence results [7] are under the assumption that all the functions involved in the objective functions are strongly convex. This lack of convergence has inspired some ADM-based methods, for example, prediction-correction type method [3,[8][9][10][11], that is, the iterate +1 1 , +1 2 , . . . , +1 is regarded as a prediction, and the next iterate is a correction for it. However, the numerical results show that the algorithm (4) always performs better than these variants. Recently, Han and Yuan [7] show that the global convergence of the extension of ADMM for ≥ 3 is valid if the involved functions are further assumed to be strongly convex. This result does not answer the open problem regarding the convergence of the extension of ADMM under the convex assumption, but it makes a key progress towards this objective.
In the following, we abuse a little the notation and still write with ; that is, the problem under consideration is where : R → R ∪ {+∞} ( = 1, 2, . . . , ) are closed proper strongly convex function with the modulus (not necessarily smooth).
The rest of the paper is organized as follows. In the next section, we list some necessary preliminary results that will be used in the rest of the paper. We then describe the algorithm formally and analyze its global convergence under reasonable conditions in Section 3. We complete the paper with some conclusions in Section 4.

Preliminaries
In this section, we summarize some basic concepts and their properties that will be useful for further discussion.
Let ‖ ⋅ ‖ denote the standard definition of the -norm, and particularly, let ‖ ⋅ ‖ = ‖ ⋅ ‖ 2 denote the Euclidean norm. For a symmetric and positive definite matrix , we denote ‖ ⋅ ‖ the -norm, that is, ‖ ‖ = √ . If is the product of a positive parameter and the identity matrix , that is, = , we use the simpler notation: If the domain of denoted by dom = { ∈ R | ( ) < +∞} is not empty, then is said to be proper. If for any ∈ R and ∈ R , we have then is said to be convex. Furthermore, is said to be strongly convex with the modulus > 0 if and only if Abstract and Applied Analysis 3 A set-valued operator defined on R is said to be monotone if and only if and is said to be strongly monotone with modulus > 0 if and only if Let Γ 0 (R ) denote the set of closed proper convex functions from R to R ∪ {+∞}. For any ∈ Γ 0 (R ), the subdifferential of which is the set-valued operator, defined by is monotone. Moreover, if is strongly convex function with the modulus , is strongly monotone with the modulus . Let be a mapping from a set Ω ⊂ R → R . Then is said to be co-coercive on Ω with modulus > 0, if Throughout the paper, we make the following assumptions.
Remark 2. Assumption 1 is a little restrictive. However, some problems can satisfy it. A remarkable one is the following route-based traffic assignment problem.
Consider a transportation network (N, ), where N is the set of nodes. We denote the set of links by A, and the number of the element of A by A , respectively. Let RS denote the set of origin-destination (O-D) pairs. For an O-D pair rs ∈ RS, let rs be its traffic demand; let rs be the set of routes connecting rs, and ∈ rs ; N rs denotes the number of the routes connecting rs; let ℎ rs be the route flow on . The feasible route flow vector ℎ = ( ∈ rs | rs ∈ RS) is thus given by ℎ rs = rs , ℎ rs ≥ 0, ∀ ∈ rs , rs ∈ RS } } } = {ℎ | (ℎ rs 1 , ℎ rs 2 , . . . , ℎ rs rs ) = rs , ℎ rs ≥ 0, ∀ ∈ rs , rs ∈ RS} .
Define as the link-route incidence matrix such that Then, link flow can be written as By denoting the link cost function as ( ) and for the additive case, the route cost function as (ℎ), they can be related by The user equilibrium traffic assignment problem can be formulated as a VI: find * ∈ such that or equivalently, find ℎ * ∈ such that where = { } is the vector of the link cost function. In general, it is easy to show that is a row of and is not a full column rank (if is, then the above variational inequality is strongly monotone).
For simplicity, in the following, we only consider the case for = 3. Notice that for ≥ 3, it can be proved similarly following the processing of = 3.

Abstract and Applied Analysis
The iterative scheme of ADMM for problem (19) is as follows: where is the Lagrangian multiplier associated with the linear constraints and > 0 is the penalty parameter.

Convergence
In this section, we prove the convergence of the extended ADMM for problem (19). As the assumptions aforementioned, by invoking the first-order necessary and sufficient condition for convex programming, we easily see that the problem (19) under the condition is characterized by the following variational inequality (VI): find * ∈ U and * ∈ ( * ) such that where :=( We denote the VI (21)-(22) by MVI(U, ).
Similarly, in [7], we propose an easily implementable stopping criterion for executing (20): and its rationale can be seen in the following lemma.
By using the notion of ‖ +1 − * ‖ 2 , from (38) we have The proof is complete.

Proof. From Lemma 5, we have
where From Assumption 1, it follows that Consequently, From (45), we have which means that the generated sequence { } is bounded. Furthermore, it follows that which means that Therefore, we have Since ‖ ‖ is nonzero and bounded, from (48) Since { } is bounded, { } has at least one cluster point, say . Let { } be the corresponding subsequence that converges 8 Abstract and Applied Analysis to . Taking a limit along this subsequence in (25) and (49), we obtain * ∈ ( * ), which follows that is an optimal Lagrange multiplier. Since * is arbitrary, we can set * = in (46) and conclude that the whole generated sequence converges to a solution of MVI(U, ).

Conclusions
In this paper, we extend the convergence analysis of the ADMM for the separable convex optimization problem with strongly convex functions to the case in which the individual functions are composites of strongly convex functions with a linear transform. Under further assumptions, we established the global convergence of the algorithm.
It should be admitted that although some problems arising from applications such as traffic assignment fall into our analysis, the problems considered here are too special. Thus, it is far away to solve the open problem of convergence of the ADMM with more than three blocks.

Assumption 1.
The solution set of (1) is denoted by X * , and it is assumed to be nonempty.

Assumption 2.
The objective function is simple. This means that, for a given constant , the following proximal problem admits a closed-form solution or can be solved efficiently with high precision: where is any given vector. At first sight, this assumption seems to be quite restrictive, but this is indeed for many problems in practice. For example, nuclear norm function in matrix completion problem, 1 -norm function in basis pursuit problem, and so forth.
Many fundamental methods have been developed over the past decades to solve problem (1). Proximal point algorithm (PPA) is one of the leading approaches for solving convex optimization problems. It is earlier used for regularized linear equations and has been applied to convex optimization by Martinet [11]. There are some significant theoretical achievements [12][13][14][15][16][17][18][19] in the field of PPA for convex optimization and monotone variational inequalities (VIs). Nowadays, it is still the object of intensive investigation [20] and leads to a variety of primal and dual methods. Common to PPA and its variants is the difficulty of their subproblems; this restricts the practical interest. Augmented Lagrangian method (ALM) [21] is a powerful method for linearly constrained problems. It can be regarded as a variant of PPA applied to the dual problem of (1). However, with the additional regularized term ‖ − ‖ 2 , its subproblems require an inverse operator of the form ( +(1/ ) ) −1 which is hard to implement in some cases. Particularly, is general or large scale, so the computation of inverse operator may fail. Hence, ALM is not sufficiently competitive when the objective function ( ) is "simple. " Extragradient method (EGM) [22] is a practical method for (1) which employs the information of current iteration. In fact, EGM is an explicit type method and requires two calls to the gradient per iteration; therefore, it might suffer slow convergence. Recently, He and Yuan [23] proposed a contraction method based on PPA (PPA-CM) to solve (1), which is elegant and simple. Inspired by PPA-CM, a Lagrangian PPA-based (LPPA) contraction method is presented in [24] which employs an asymmetrical proximal term [25]. These two PPA-based methods have nice convergence properties that are similar in many ways to PPA, but they both require a quite restrictive condition for convergence and are sensitive to the initial choice of stepsizes.
In this paper, we focus on development of PPA-type method for solving (1). Based on LPPA, we propose a general self-adaptive relaxed-PPA method (SRPPA) which is simple yet efficient. Our algorithm capitalizes certain features of PPA, hence, we term it relaxed-PPA. The proposed algorithm has several nice fronts. First, our method is a PPA-type method with asymmetrical linear term, which is clearly a different nature to classical PPA. It relaxes the jointly structure of subproblem to a tractable one. Second, we provide two simple search directions for new iterate. Third, the stepsize choice is flexible, and the condition for convergence guarantee is weaker than both PPA-CM and LPPA. Finally, simple adaptive strategies are employed to choose stepsize, and this appealing aspect is significantly important in practice. We also demonstrate that our method is relevant for various applications whose practical success is made possible by our efficient algorithm. This paper is organized as follows. In Section 2, we provide some notations and preliminaries which are useful for subsequent analysis. In Section 3, we review some related works. The general relaxed-PPA and its variant are formally presented in Section 4. Self-adaptive strategies to choose stepsize are also described. Next, in Section 5, the convergence analysis is provided. In Section 6, we present some concrete applications of (1) and elaborate on the implementation of our method; preliminary numerical results are also reported to verify the efficiency of our proposed method. Finally, in Section 7, we conclude the paper with a discussion about the future research directions.

Preliminaries
In this section, we first establish some important notations used throughout this paper. Then, we describe the variational inequality formulation of (1) which is convenient for the convergence analysis.
Recalling the monotonicity of , it is easy to get that VI(Ω, ) (9) is monotone. Since the solution set of (1) is assumed to be nonempty, the solution set of VI(Ω, ), denoted by Ω * , is also nonempty. Our analysis will be built upon this equivalent VI formulation.

The Existing Related Methods
There are basically two lines of research for VI(Ω, ) (9), either deal with it by implicit methods that are in general computationally intractable or concentrate on relaxing it with explicit methods. In this section, we first briefly review the well-known classical PPA and EGM. And then, PPA-CM [23] and LPPA [24] are discussed, which will provide motivation for our general self-adaptive relaxed-PPA.

Classical PPA for the Equivalent Variational Inequality.
PPA and its variants are implicit methods which have fast asymptotical convergence rate and produce highly accurate solutions. At each iteration, the subproblem of classical PPA consists of a regularized term, which can be expressed as Then, the update step is taken as follows: PPA has a nice convergence property Although classical PPA is conceptually appealing, subproblem (10) presents certain computational challenges. More specifically, primal variable and dual variable are tied together, and their subproblems are treated as a joint problem.
In most cases, this joint subproblem may be as difficult as the original problem (9). As a result, PPA is "conceptual" rather than implementable. It lacks capability in exploiting potential decomposable/specific structure of subproblem. Variants of classical PPA have been explored in the literature, in order to solve the proximal subproblem (10), inexactly, see, for example, [14,15,17,19]. Unfortunately, inexact variants take expensive computation for obtaining approximative solutions.

The Methods Based on the Simplest Relaxation.
To overcome the drawbacks of the classical PPA, it is instinctive to relax subproblem (10) to a solvable one. The most straightforward and simplest relaxation is to replace (̃) with ( ) in the proximal subproblem (10), which amounts to the following subproblem: The solution of the relaxed problem (13) is given bỹ= . It is clear that methods with such relaxation are explicit type methods. However,̃cannot be accepted directly as the new iterate. Using the terminology "predictor-corrector, " such point can be viewed as a predictor. Here, we list two simple methods which employ predictort o obtain corrector as the new iterate.
(i) The extragradient method (EGM) updates the new iterate (corrector) by (ii) The projection and contraction methods (PCM) [28][29][30] perform update as follows: The sequence { } generated by the above mentioned EGM or PCM satisfies which is similar to PPA. Both EGM and PCM use the simplest relaxation to obtaiñin th iteration, hence are computationally practical. These methods have appealing practical aspects; however, such simplest relaxation does not exploit the inner structure of VI(Ω, ) (9). This observation prompts the need for dedicated relaxations.

PPA-Type Contraction
Method. The algorithms that are closely related to ours are PPA-CM [23] and LPPA [24]. The PPA-CM obtains the predictor̃by solving the following subproblem: find (̃,̃) ∈ Ω such that And perform the update The framework of LPPA is as follows: where = ( 0 ) .
And the new iterate is defined by Both procedures are simple and can solve subproblem efficiently; but their nice convergence results require a quite restrictive condition, that is; > ‖ ‖ in PPA-CM and > (1/2)‖ ‖ in LPPA, respectively. The stepsizes , are directly determined by such condition; hence, it is important to estimate ‖ ‖. Overestimation may lead to poor stepsizes and slow convergence, while underestimation may result in divergence. In addition, they are both quite sensitive to the choice of , . To overcome those drawbacks, we propose a general self-adaptive relaxed-PPA, and as mentioned earlier, it can provide improved guarantee for convergence and has potential progress in the choice of stepsize. Furthermore, selfadaptive strategies are designed to improve performance. The relaxed-PPA described here involves two steps. First, we solve the relaxed subproblem ( * ), ( * * ) to obtain predictor, which is nice and efficient for the nature of the problem under study. Note that the -predictor step ( * ) involves minimizing plus a convex quadratic function, and under Assumption 2, it can be efficiently solved or it admits a closed form solution. And then, -predictor step ( * * ) is just a projection onto Λ which is tractable and computationally efficient. It is clear that the prediction step employs a Gauss-Seidel manner to update information efficiently. The correction step ($) only involves matrix-vector multiplication which is very simple and straightforward.
Remark 3. We first make some insight into the correction step in Algorithm 1. The obtained̃plays no direct role as the new iterate. Instead, we need some kind of "corrector" defined in ($). Although matrix in ($) is just a required positive symmetry definite, our goal here is to fully integrate the information of and̃to construct effective, simple search direction −1 ( −̃) for the corrector. Based on this consideration, we elaborately provide two simple choices of . . This case is a little less intuitive, but it can lead to a simple update form as well as Case 1. The underlying derivation is a little more complicate. Applying in ($), we get Recalling is a lower triangular matrix, by the fact that its inverse is also a lower triangular, we have Plugging the previous relationship to (31), we have ) . (30) In fact, this is a scheme of Gaussian back substitution.
Both cases only involve one matrix-vector multiplication which makes the update form simple. And the computational cost is usually inexpensive.
and its compact form We observe that subproblem (32) is similar to (10) in PPA, except for the construction of asymmetry matrix . As mentioned before, (32) is the same as the prediction subproblem in [24]. Even though they are closely related, the stepsize choice here is quite different. We provide more specific and weaker condition for stepsize , . It is clear that condition (#) does not need prior knowledge of matrix . Furthermore, it only involves matrix-vector multiplication, and so, it is easy to verify, and it is amenable to large-scale . If , fail to meet this convergence condition (#), one should appropriately increase , . In the following subsection, we will elaborate on the self-adaptive strategies to increase the stepsizes. At this point, condition (#) may be seen somewhat unmotivated. Some insight into this will be provided later, as we proceed with the convergence analysis. The convergence condition in [24] has a quite different feature: , satisfy It is stronger than condition (#). The following lemma is devoted to the proof of this result.
Proof. Note that Recall that matrix described in Algorithm 1 can be designed in two different cases.
Abstract and Applied Analysis According to Cauchy-Schwarz inequality, we get The last inequality follows directly from condition (33).
, by the definition of ( ,̃), we obtain Then, we get The first inequality follows from the Cauchy-Schwarz inequality, and the last one follows directly from condition (33).
Condition (33) is not only stronger than Condition (#), but it also requires that matrix is positive semidefinite, while condition (#) does not. Furthermore, condition (33) may require the explicit expression of or knowledge of ‖ ‖. Despite these drawbacks, condition (33) is appealing to the problems in which ‖ ‖ is known beforehand or easy to compute/obtain. For instance, is small scale, an identity matrix or a projection operator, and so forth. It is clear that both condition (#) and (33) are more flexible than the one in PPA-CM [23]. The most aggressive condition (#) may lead to further improvement in stepsize choice. Moreover, it is worthwhile to notice that condition (#) is elegantly designed and provides ( ,̃) with favourable property. In fact, for general matrix , condition (#) also can guarantee convergence.
Remark 6. The update stepsize plays an important role here. In fact, it can be regarded as an optimal stepsize which will be illustrated in the following section.

Remark 7.
We should restrict the adjustment in Step 5 of Algorithm 1 within a limited number to avoid divergence.
In Algorithm 1, we carry out the -predictor before performing -predictor. The roles of and are symmetric; hence, sweeping the order will not break the Gauss-Seidel structure. We switch and and obtain a variant of relaxed-PPA with the order of the -predictor step and -predictor step reversed. This variant is illustrated in Algorithm 2. However, there is no a priori information to know which algorithm is superior. Here, we let

Adaptive Enhancements.
Both PPA-CM and LPPA employ fixed stepsizes , . Experiments reveal that they will suffer slow convergence when stepsizes , are chosen inappropriately. A natural question is, how to choose the proper initial stepsizes , . Here, we propose self-adaptive strategies with the goal of improving the convergence in practice, as well as making performance less dependent on the initial choice of stepsizes. Our strategies dynamically incorporate the information of the current iteration to perform more informative choice of stepsizes for the next iteration [31]. When doing so, the algorithm will be adaptive and free from the initial choice. Denote and then, (31) can be rewritten as Under -norm, the quantity can be viewed as a residual for the dual feasibility condition, and can be viewed as a primal residual. These two residuals converge to zero as relaxed-PPA proceeds. Note that Abstract and Applied Analysis And this implies that small values of tend to reduce the primal residual but have potential to enlarge violations of dual feasibility and tend to produce larger dual residual. This observation motivates us to balance primal and dual residuals. When condition (#) fails, we increase stepsizes , properly according to the adaptation shown in Algorithm 3.
Here, 1 > 1, 2 > 1. This adaptation strategy increases when the dual residual ‖ ‖ 2 appears large compared to the primal residual ‖ ‖ 2 and increases when the dual residual ‖ ‖ 2 seems too small relative to the primal residual ‖ ‖ 2 .
As mentioned, condition (33) is stronger than condition (#). If one chooses condition (33), our RPPA also converges. It must have predetermined stepsizes satisfying = ‖ ‖ (here, ≥ 0.5). However, there is no priority knowledge of the choice of individual or . Here, we can also adjust , automatically when choosing condition (33). Intuitively, we consider expansion of the entire residual under -norm: there are three terms involving or , and we intend to balance these terms. For fixed , take = ( / )‖ ‖; then Applying (44) into (43), clearly we have Now, we consider adjusting stepsize to balance ‖ −̃‖ 2 and (1/ )( ‖ ‖‖ −̃‖ 2 + ‖ ( −̃)‖ 2 ) and obtain another adaptation strategy (see Algorithm 4). It is worth noting that too many adjustments of stepsizes by Algorithm 4 might cause the algorithm to diverge in practice, and we therefore restrict these adaptations within a limited number of iterations. If one chooses Algorithm 4, there is no need to carry out Step 5 in Algorithm 1 (or Algorithm 2). These techniques embedded into relaxed-PPA automatically choose a "better" stepsize for the next iteration.

Convergence Analysis
In this section, we analyze convergence of our primal-dual relaxed-PPA. The convergence analysis of dual-primal scheme can follow a similar procedure.
Combining Lemmas 8 and 9, we immediately obtain the following convergence theorem. Proof. First, for each ∈ Ω, we have It follows from (57) and consequently, This means that ∞ is a solution of VI(Ω, ). Note that the inequality (57) is true for all solution points of VI(Ω, ), and hence, we have Sincẽ→ ∞ ( → ∞) and −̃→ 0 ( → ∞), for any given > 0, there exists an integer > 0 such that Therefore, for any ≥ , it follows from (63) and (64) that This implies that the sequence { } converges to ∞ , which is a solution of VI(Ω, ) (9).

Applications and Preliminary Numerical Experiments
The general self-adaptive relaxed-PPA (SRPPA) offers a flexible framework for solving many interesting problems. We illustrate our algorithm with different applications: basis pursuit problem, nearest correlation matrix problem. In this section, we describe the results of experiments whose goal is to demonstrate the efficiency of general relaxed-PPA (RPPA) and its self-adaptive version. To that end, we compare RPPA with certain state-of-the-art algorithms on different problems. Our experiments focus on efficiency and speed of convergence and evaluate the methods in terms of their number of iterations and computational times. All the codes were written by Matlab R2009b version, and all the numerical experiments were performed on a Lenovo desktop computer with Intel (R) Core (TM) i5 CPU with 3.2 GHz and 3.5 GB RAM.
And here, ‖ ⋅ ‖ 1 denotes the 1 norm defined as ‖ ‖ 1 := ∑ =1 | |. BP is a fundamental problem in image processing and modern statistical signal processing, particularly the theory of compressed sensing; see, for example, [1][2][3][4] for intensive study. We now discuss our approach to BP problem of over-complete representations. Our experiments in this subsection use synthetic data which were mainly designed to illustrate the nice performance of our RPPA. The synthetic problem that we test here is similar to the one employed in [32]. We generate the data as follows: matrix is a random × matrix, with Gaussian i.i.d. entries of zero mean and variance 1, with = /2. original ∈ is the original sparse signal, constructed with /5 nonzero values, randomly from standard normal distribution. We use original to generate the measurements as = original . It is desirable to use test problems that have a precisely known solution. In fact, when original is very sparse, it is the solution to the minimization problem (66). Hence, in our synthetic problem, original is exactly the solution.
In our first experiment, we compared general RPPA using two different 's mentioned in Section 4.1. For BP problem, we use condition (#) and Algorithm 3. Since constructed here is a general random matrix, and when is large scale, ‖ ‖ might be obtained costly. A simple stopping criterion was used in this experiment, and the stopping tolerance Tol was set to 10 −10 . In all the tests, initial stepsizes were set as = 10, = 1, the primal variable 0 was initialized as zeros( , 1), and the dual 0 was ones( , 1) in Matlab. Table 1  are Gaussian back substitute form methods and perform a little slower than Gaussian type methods. We also plot a figure to graphically illustrate the performance of four SRPPAs. Figure 1 shows the results from the test with = 1000 and = 6000, depicting error versus CPU time. Qualitywise, SPDRPPAG1-I was on par with SDPRPPAG1-I.
In the second experiment, we compare the performance of SPDRPPAG1-I with TFOCS (source code can be found at http://cvxr.com/tfocs/) [32], ADMM (source code can be found at http://www.stanford.edu/∼boyd/papers/admm/) [33], and PPA-CM. To make the comparison independent of the stopping criterion for each algorithm, we first run TFOCS to get its solution TFOCS and set a benchmark error benchmark err . = TFOCS − original 2 and then run other algorithms until they obtain smaller errors than this benchmark. TFOCS was stopped upon 10 Abstract and Applied Analysis  Since we found that Tol = 10 −12 is small enough to guarantee very high accuracy, we set Tol = 10 −12 in TFOCS. The parameters of TFOCS and ADMM were taken with their defaults. To guarantee the convergence, fixed stepsizes , were set to = 100, = 1.01 * ‖ ‖/ for PPA-CM. In SPDRPPAG1-I, we also choose the same convergence condition (#) and initial step size = 10, = 1 as the previous experiment. We varied the size of from = 512 ( = /2) to = 8500. The results of this experiment are summarized in Table 2. There, we report the run time in seconds, the number of iterations, and the error of the recovery solution. In Table 2, "-" means "out of memory. " We observe from Table 2 that four algorithms reach high accuracy around 10 −9 . SPDRPPAG1-I is about two times faster than the first-order method implemented in the TFOCS package, and moreover, it usually outperforms TFOCS in terms of iterations. For medium size problems, SPDRPPAG1-I is clearly faster than ADMM. Even for small size problems, SPDRPPAG1-I shows its superior performance. The main reason lies in that ADMM computed ( + ) −1 to solve its subproblem exactly which would take expensive computational cost. Not surprisingly, the general SPDRPPAG1-I performs better than the primary PPA-CM. Here, the total iterations of SPDRPPAG1-I are less than 50% of PPA-CM. As we have mentioned, "optimal" update stepsize and more flexible condition for convergence may provide SPDRPPAG1-I improved performance. SPDRPPAG1-I is faster than PPA-CM in terms of CPU times. However, the superiority of CPU time is not so significant as iteration number. For the cases = 4300, it is just about 62% of PPA-CM. This is not particularly surprising; compared to PPA-CM, SPDRPPAG1-I has to take extra computation for convergence condition and "optimal" in each iteration.

Nearest Correlation Matrix
Problem. The nearest correlation matrix problem is solving the problem where ∈ R is the vector whose entries are all 1 , + denotes the cone of positive definite symmetric matrices, diag( ) is the vector of diagonal elements of , and ‖ ⋅ ‖ denotes the matrix Fröbenius norm ‖ ‖ = trace ( ) 1/2 . Here, we apply PPA-CM, LPPA, and SDPRPPA1-II for solving (70). The standard Matlab Mex interface mexsvd is used to conduct the eigenvalue decomposition. We constructed test data sets and stopping criterion like those of [24]. As mentioned in the prequel, we expect our SRPPA to produce robust performance. To assess the effectiveness of the adaptive strategies proposed in Section 6, we now move on to the description of experiments that focus on the consequences of the initial stepsizes. For investigating, we used dimensions ∈ {500, 1000} and varied from 0.05 to 100, and initial points were set to 0 in all cases. Note that = ; we fixed = 1.01/ for PPA-CM, = 0.65/ for LPPA and chose = 0.65/ as initial start for SDPRPPAG1-II. Since the experiments with other values of give qualitatively similar results, we therefore do not plot those results to avoid clutter in the figures. The respective numerical results are plotted in Figure 2.
It is clear that, for PPA-CM and LPPA, the convergence performance was a result of the stepsize selection. They are both fairly sensitive to initial stepsize (or ). The results confirm that, with inappropriate stepsizes, both PPA-CM and LPPA become significantly slow. SDPRPPAG1-II yields significantly robust performance with adaptive strategy. And it is independent of the initial stepsizes and illustrates its superior performance. Furthermore, SDPRPPAG1-II yields competitive results even when PPA-CM and LPPA chose the "good" initial stepsize. This underlies the importance of adaptive strategy in producing good performance. Of course, care should be taken. For instance, the cost of computing optimal stepsize here is negligible, compared to the computation of SVD; when they are more costly, general LPPA will be expected to perform slower than PPA-CM.

Conclusions
In this paper, we proposed an efficient general self-adaptive relaxed-PPA method for linearly constrained convex programming and provided theoretical convergence analysis for this method. The stepsizes choice condition is flexible and simple. Self-adaptive strategies are provided to make our method more efficient and robust. Experiments of the method in comparison to other state-of-art methods are provided to confirm the efficiency of the proposed method. Our numerical results suggest that SRPPA is effective and simple to implement. There are a few directions for further research, but we list here only two. The first is the question of whether we may modify the algorithm to work with more general constrained convex problems. Second, we aim to provide various relaxations of the subproblem for the practical purpose.

Introduction
During the last few decades, the celebrated Banach contraction principle, also known as the Banach fixed point theorem [1], has become one of the core topics of applied mathematical analysis. As a consequence, a number of generalizations, extensions, and improvement of the praiseworthy Banach contraction principle in various direction have been explored and reported by various authors; see, for example,  and the references therein. In parallel with the Banach contraction principle, Kannan [5] and Chatterjea [6] created, respectively, different type, fixed point theorems as follows. Then has a unique fixed point in .
The characterization of the renowned Banach fixed point theorem in the setting of multivalued maps is one of the most outstanding ideas of research in fixed point theory. The remarkable examples in this trend were given by Nadler [2], Mizoguchi and Takahashi [3], and M. Berinde and V. Berinde [4]. On the other hand, investigation of the existence of a fixed point of non-self-maps under certain condition is an interesting research subject of metric fixed point theory, see, for example, [19][20][21][22][23][24][25][26][27], and references therein.
The following attractive result was reported by M. Berinde and V. Berinde [4] in 2007.
Then has a fixed point in .
If we take = 0 in Theorem 3, then we conclude the remarkable result of Mizoguchi and Takahashi [3] which is a partial answer of problem 9 in [8]. Then has a fixed point in .
Recently, Du [12] established the following theorem which is an extension of Theorem 3 and hence Theorem 4.
Then has a fixed point in .
The basic objective of this paper is to investigate the existence of coincidence and fixed points of multivalued nonself-maps under the certain conditions in the setting of metric spaces. The presented results generalize, improve, and extend several crucial and notable results that examine the existence of the coincidence/fixed point of well-known maps such as Kannan type, Chatterjea type, Mizoguchi-Takahashi type, Berinde-Berinde type, Du type, and other types in the context of complete metric spaces.

Preliminaries
is said to be the Hausdorff metric on CB( ) induced by the metric on . It is also known that (CB( ), H) is a complete metric space whenever ( , ) is a complete metric space. Let be a nonempty subset of , : → a singlevalued map, and : → N( ) a multivalued map. A point in is a coincidence point of and if ∈ . If = is the identity map, then = ∈ and call a fixed point of . The set of fixed points of and the set of coincidence point of and are denoted by F ( ) and COP ( , ), respectively. In particular, if ≡ , we use F( ) and COP( , ) instead of F ( ) and COP ( , ), respectively. Throughout this paper, we denote by N, and R, the set of positive integers and real numbers, respectively.
In what follows that, we recall some characterizations of MT -functions proved first by Du [12].

Existence Theorems of Coincidence Points and Fixed Points for Multivalued Non-Self-Maps of Kannan Type and Chatterjea Type
In this section, we prove the existence of coincidence points and fixed points of multivalued non-self-maps of Kannan type and Chatterjea type. For this purpose, we first established a new intersection theorem of COP ( , ) and F ( ) for multivalued non-self-maps in complete metric spaces.
Proof. It is obvious that any of these conditions (K1)-(K7) implies condition (D3) as in Theorem 8. So the desired conclusion follows from Theorem 8 immediately.
The following fixed point theorem for multivalued nonself-maps of generalized Kannan type can be established immediately from Theorem 9 for ≡ (the identity mapping).
Then F ( ) ̸ = 0. As a consequence of Theorem 10, we obtain the following generalized Kannan type fixed point theorems for multivalued maps.
Applying Theorem 14, we can prove the following fixed point theorems for multivalued maps of generalized Chatterjea type.

New Coincidence and Fixed Point Results
for Various Multivalued Non-Self-Maps:
Then COP ( , ) ∩ F ( ) ̸ = 0. As a direct consequence of Theorems 19 and 20, we obtain the following fixed point result for multivalued non-selfmaps of Du type in complete metric spaces.  Step 1. For given { }, let the sequence {V } be generated iteratively by where proj C is the metric projection and { } is a real number sequence.
Then the sequence { } generated by Algorithm 8 converges strongly to * ∈ Υ.

Introduction
The notion of coupled fixed point was introduced by Guo and Lakshmikantham [1] in 1987. In a recent paper, Gnana Bhaskar and Lakshmikantham [2] introduced the concept mixed monotone property for contractive operators of the form : × → , where is a partially ordered metric space, and then established some coupled fixed-point theorems. After that, many results appeared on coupled fixedpoint theory in different contexts (see, e.g., [3][4][5][6]). Later, Berinde and Borcut [7] introduced the concept of tripled fixed point and proved tripled fixed-point theorems using mixed monotone mappings (see also [8][9][10]).
Very recently, Roldán et al. [11] proposed the notion of coincidence point between mappings in any number of variables and showed some existence and uniqueness theorems that extended the mentioned previous results for this kind of nonlinear mappings, not necessarily permuted or ordered, in the framework of partially ordered complete metric spaces, using a weaker contraction condition, that also generalized other works by Berzig and Samet [12], Karapınar and Berinde [13].
Partial metric spaces were firstly introduced by Matthews in [14] as an attempt to generalize the metric spaces by establishing the condition that the distance between a point to itself (which is not necessarily zero) is less or equal than the distance between that point and another point of the space. In the mentioned papers, Matthews studied topological properties of partial metric spaces and stated a modified version of a Banach contraction mapping principle on this kind of spaces. After Matthews' pioneering work, the theory of partial metric spaces and particularly the field of fixed-point theorems have expansively been developed due to the increasing interest in this area and motivated by its possible applications (see [15,16] and references therein).
In this paper, our main aim is to study a weaker contractivity condition for nonlinear mappings of any number of arguments. This condition can be particularized in a variety of forms that let us extend the previously mentioned results and other recent ones in this field (see [2,5,7,9,11,12,[16][17][18][19][20]). We also notice that our results cannot be obtained by the very recent paper of Haghi et al. [21] (for more details see Remark 26).
From these properties, we can easily deduce that ( , ) ≥ 0 and ( , ) = ( , ) for all , ∈ . The last requirement is called the triangle inequality. If is a metric on , we say that ( , ) is a metric space (for short, an MS).
Definition 1 (see [22]). A triple ( , , ≤) is called a partially ordered metric space if ( , ) is a MS and ≤ is a partial order on .
Definition 2 (see [2]). An ordered MS ( , , ≤) is said to have the sequential -monotone property if it verifies (i) if { } is a nondecreasing sequence and { } → , then ≤ for all ; (ii) if { } is a nonincreasing sequence and { } → , then ≥ for all .
If is the identity mapping, then is said to have the sequential monotone property.
If ( , ≤) is a partially ordered space, , ∈ , and ∈ Λ , we will use the following notation: Let : → and : → be two mappings.
In this case, ( , ) is a partial metric space (for short, a PMS).
define partial metrics on , where : → R + 0 is an arbitrary function and ≥ 0.
Obviously, if ( , ) is a MS and we define = , then ( , ) is a PMS. Indeed, a partial metric on verifies for all , ∈ , are (usual) metrics on . On a PMS, the concepts of convergence, Cauchy sequences, completeness, and continuity are defined as follows.
We have used the previous notation because we need to distinguish between -convergence and -convergence on and usual convergence for real sequences.
Proof. We know that

Remark 25.
In the previous theorem, if the image Im of the metric is not the whole set [0, ∞[, then and can only be defined on Im , and we can consider a wider range of mappings since it is only necessary to impose that they are continuous and nondecreasing on Im .
Step 2. We claim that ( 1 , 2 , . . . , ) is the unique Υcoincidence point of and such that = for all . It is similar to Step 2 in Theorem 11 in [11].
It is natural to say that is injective on the set of all Υcoincidence points of and when = for all implies = for all when ( 1 , 2 , . . . , ), ( 1 , 2 , . . . , ) ∈ are two Υ-coincidence points of and . For example, this is true is is injective on .

Corollary 31.
In addition to the hypotheses of Theorem 30, suppose that is injective on the set of all Υ-coincidence points of and . Then, and have a unique Υ-coincidence point.
Consider : → and : → defined by It is not difficult to prove the following statements.
(1) and are -continuous mappings (since generates the discrete topology on ).

Introduction
Throughout this paper, we assume that is a real Banach space with the norm ‖ ⋅ ‖, * the dual space of , ⟨⋅, ⋅⟩ the duality between and * , and a nonempty closed convex subset of . R + denotes the set of nonnegative real numbers and N the natural number set. The mapping : → 2 * with ( ) = { * ∈ * : ⟨ , * ⟩ = ‖ ‖ 2 , * = ‖ ‖} , ∈ , is called the normalized duality mapping. Let : → be a nonlinear mapping. ( ) denotes the set of the fixed points of .
As we know, a mapping : → is said to be pseudocontractive, if, for all , ∈ , there exists ( − ) ∈ ( − ), such that Variational inequalities introduced by Stampacchia in the early sixties have had a great impact and influence on the development of almost all branches of pure and applied sciences and have witnessed an explosive growth in theoretical advances, algorithmic development, and so forth.
Recently, some authors also studied the problem of finding the solution set of variational inequalities and the common element of the fixed point set for generalized nonexpansive mappings in the framework of real Hilbert spaces and Banach spaces. As is known to all, the variational inequality problem, nonlinear optimization problem, and fixed point problem are equivalent to each other under certain conditions.
In 2012, Chang et al. [1] introduced a more general class of pseudocontractive mappings and studied the methods for approximation of the split common fixed points.
(II) A mapping : → is said to be ( , )asymptotically strictly pseudocontractive, if there exist