# Monte Carlo Sampling Method for a Class of Box-Constrained Stochastic Variational Inequality Problems

Authors: Pei-Yu Li
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

This paper uses a merit function derived from the Fishcher–Burmeister function and formulates box-constrained stochastic variational inequality problems as an optimization problem that minimizes this merit function. A sufficient condition for the existence of a solution to the optimization problem is suggested. Finally, this paper proposes a Monte Carlo sampling method for solving the problem. Under some moderate conditions, comprehensive convergence analysis is included as well.


---

## Body

## 1. Introduction

Letl and u be two n-dimensional vectors with components li∈ℝ∪−∞ and ui∈ℝ∪+∞ satisfying li<ui, and denote, by S, the nonempty and possibly infinite box l,u≔x∈ℝn|li≤xi≤ui,i=1,…,n. Then, the box-constrained variational inequality problem (BVIP, for short) is to find a vector x∗∈S such that(1)x−x∗TFx∗≥0,∀x∈S,where F:ℝn⟶ℝn is a given function. This problem is also called the mixed complementarity problem [1]. Both linear complementarity problems [2] and nonlinear complementarity problems [3] have played an important role in studying economic equilibrium and engineering problems [4, 5].Much effort has been made to derive merit functions for the BVIP and then, using these functions, to develop solution methods. Recently, Kebaili and Benterki [6] propose a penalty approach for a box-constrained variational inequality problem (BVIP). It is replaced by a sequence of nonlinear equations containing a penalty term. A homotopy method for solving mathematical programs with box-constrained variational inequalities is presented [7]. A reformulation of the BVIP, based on the Fishcher–Burmeister function, is described in paper [8] by Sun and Womersley. Further reformulations can be obtained by replacing the Fishcher–Burmeister function by the introduced function from [9] which seems to have somewhat stronger theoretical properties and a better numerical behavior.Stochastic variational inequality problem (SVIP, for short) model is a natural extension of deterministic variational inequality models. Over the past few decades, deterministic variational inequality has been extensively studied for its extensive application in engineering, economics, game theory, and networks; see the book on the topic by Facchinei and Pang [1]. While many practical problems only involve deterministic data, there are some important instances where problem data contain some uncertain factors and consequently SVIP models are needed to reflect uncertainties. Gürkan et al. [10] have shown how to extend a simulation-based method, sample-path optimization, to solve SVIP, and Robinson [11] has provided a mathematical justification for sample-path optimization, while Shapiro et al. [12] suggested Monte Carlo sampling methods for SVIP. SVIP can be found in pricing game [13] and inventory competition [14] among several firms that provide substitutable goods or services. Some stochastic dynamic noncooperative games [15] and competitive Markov decision processes [16] can be formulated as examples of SVIP. Jiang and Xu [17] proposed a stochastic approximation method for numerical solution of SVIP. The method is an iterative scheme where, at each iterate, a correction is made and the correction is obtained by sampling or other stochastic approximation. Xu [18] applied the well-known sample average approximation (SAA, for short) method to solve the same class of stochastic variational inequality problems (SVIP).Aimed at a practical treatment of the SVIP, box-constrained stochastic variational inequality problem (BSVIP, for short) is meaningful and interesting to study [19]. Motivated by Sun and Womersley [8], Luo and Lin [20] formulated a class of BSVIP as an optimization problem that minimizes the expected residual of the merit function derived on the Fishcher–Burmeister function and proposed a Monte Carlo sampling method for solving the problem.In this paper, we consider the following BSVIPF,S to find a vector x∗∈S such that(2)x−x∗TEFx∗,ξω≥0,∀x∈Sω∈Ωa.s.,where ξ:Ω⟶Ξ⊂ℝm is a random vector defined on probability space Ω,ℱ,P, F:S×ℝm⟶ℝn is a mapping, E denotes the mathematical expectation with respect to the distribution of ξω, and “a.s.” is the abbreviation for “almost surely” under the given probability measure. To ease the notation, we will write ξω as ξ, and this should be distinguished from ξ being a deterministic vectors of Ξ in the context.We are concerned with the numerical solution of BSVIP. If we are able to obtain a closed form ofEFx∗,ξ, then BSVIP becomes a deterministic BVIP, and the existing numerical methods for the latter [8] can be applied directly. However, in practice, obtaining a closed form of EFx∗,ξ or computing the value of it numerically is often difficult either due to the unavailability of the distribution of ξ or because it involves multiple integration.Motivated by the above work, we make use of a merit function derived from the Fishcher–Burmeister function (FB-function, for short) and formulate BSVIP as an optimization problem that minimizes this merit function. We study a sufficient condition for the existence of a solution to the optimization problem. Finally, we propose a Monte Carlo sampling method for solving the problem. Under some moderate conditions, comprehensive convergence analysis is included as well.The organization of this paper is as follows. In Section2, we study some preliminary knowledge. Section 3 shows under what conditions the level sets of the merit function are bounded. In Section 4, we make use of a Monte Carlo sampling method to handle the expectation; moreover, we establish convergence of global optimal solutions of approximation problems generated by the proposed method. Preliminary numerical results are reported in Section 5. Finally, we give some conclusions.


## 2. Preliminaries

We first consider a merit function derived from the NCP-function [21] for BVIP (1). It is easy to see that BVIP is equivalent to its Karush–Kuhn–Tucker (KKT, for short) system:(3)Fx+∑ui≠+∞λiei−∑lj≠−∞μiej=0,0≤λi⊥xi−ui≤0,ui≠+∞,i=1,…,n,0≤μj⊥xj−lj≥0,lj≠−∞,j=1,…,n,where ei is a n-dimensional unit vector. If x∈Rn solves BVIP, then x,λ,μ∈Rn×Rn×Rn solve the KKT system (3). Conversely, if x,λ,μ∈Rn×Rn×Rn solve the KKT system (3), then x solves BVIP (1).Letϕ be a NCP-function. Then, the KKT system (3) can be written as(i)
Ifli=−∞,ui=+∞,1≤i≤n, then Fix=0.(ii)
Ifli≠−∞,ui=+∞,1≤i≤n, then μi=Fix,0≤Fix⊥xi−li≥0, and we have(4)ϕFix,xi−li=0.(iii)
Ifli=−∞,ui≠+∞,1≤i≤n, then λi=−Fix,0≥Fix⊥xi−ui≤0, and we have(5)ϕ−Fix,ui−xi=0.(iv)
Ifli<xi<ui,1≤i≤n, then(6)Fix+λi−μi=0,0≤λi⊥xi−ui≤0,0≤μi⊥xi−li≥0,⇔0≤λi⊥xi−ui≤0,0≤Fix+λi⊥xi−li≥0,and we have(7)ϕλi,ui−xi=0,ϕFix+λi,xi−li=0,⇔ϕλi,ui−xi+ϕFix+λi,xi−li=0.Based onϕ, we let Φix as(8)Φix≔Fix,ξ,ifli=−∞,ui=+∞,ϕFix,ξ,xi−li,ifli≠−∞,ui=+∞,ϕ−Fix,ξ,ui−xi,ifli=−∞,ui≠+∞,ϕλi,ui−xi+ϕFix+λi,xi−li,ifli≠−∞,ui≠+∞.Then, it is easy to see that(9)x∗solves BVIP⇔x∗solvesΦx=0.In turn, corresponding to this reformulation of the BVIP, we may define the merit function:(10)Ψx≔12ΦxTΦx,Therefore,Ψx is an unconstrained differentiable merit function for BVIP (1). In order to define this merit function, we may consider the FB-function, which is given by(11)ϕa,b=a+b−a2+b2,and let us introduce a partition of the index set I:(12)If≔i∈I|−∞=li<ui=+∞,Il≔i∈I|−∞<li<ui=+∞,Iu≔i∈I|−∞=li<ui<+∞,IF+≔i∈I|−∞<li<ui<+∞,EFix,ξ≥0,IF−≔i∈I|−∞<li<ui<+∞,EFix,ξ<0.Consequently, we give the merit function for BSVIP (2) as follows:(13)θx=12∑i∈IfEFix,ξ2+∑i∈Il∪IF+ϕ2EFix,ξ,xi−li+∑i∈Iu∪IF−ϕ2−EFix,ξ,ui−xi.That is to say, we formulate BSVIP as the following optimization problem:(14)min,θx,s.t.,x∈S.Throughout, we assume thatFx,ξ is continuous with respect to x for any ξ∈Ξ. Furthermore, suppose that the sample space Ξ is nonempty, and for every x∈S,(15)EFx,ξ2<+∞,E∇xFx,ξ<+∞,EFx,ξ∇xFx,ξ<+∞,where ⋅ means the Euclidean norm. By the above assumption and Th16.8 in [22], we get that fx≔EFx,ξ is continuous with respect to x.


## 3. Boundedness of Level Sets

In this section, we discuss conditions for boundedness of the level sets of the merit function (13). Consider the level set defined by(16)LθSc≔x∈S|θx≤c,where c≥0 is a given scalar.Definition 1.
A functionf:ℝn⟶ℝ is a uniform P-function if there exists a positive constant μ such that, for every x and y in ℝn,(17)max1≤i≤nxi−yifix−fiy≥μx−y2.
Now let us focus on the properties of FB-functionϕa,b.Lemma 1.
For givena,b∈R, we have FB-function ϕa,b satisfying(18)23+22min2a,b≤ϕ2a,b≤6+42min2a,b.Proof.
From Tseng [23], for any two numbers a,b∈R, we have(19)22+2mina,b≤ϕa,b≤2+2mina,b.
Then,(20)23+22min2a,b≤ϕ2a,b≤6+42min2a,b.
Then, we study the boundedness of the level sets.Theorem 1.
Suppose thatfx≔EFx,ξ is a uniform P-function. Then, for any c≥0, LθSc is bounded.Proof.
Suppose that there is a nonnegative numberc¯ such that LθSc¯ is unbounded. This implies that there exists a sequence xk⊆LθSc such that limk⟶∞xk=+∞. We first define the index set I⊆1,…,n by(21)I≔i|limk⟶∞xik=+∞.
By assumption, we haveI≠∅. We now define a new sequence yk⊆ℝn as follows:(22)yik≔0,ifi∈I,xik,ifi∉I.
From the definition ofyk and the fact that f is a uniform P-function, we obtain that, for some constant μ>0,(23)μ∑i∈Ixik2=μxk−yk2≤max1≤i≤nxik−yikfixk−fiyk,=maxi∈Ixikfixk−fiyk≤∑i∈Ixik2∑i∈Ifixk−fiyk,from which we obtain(24)μ∑i∈Ixik2≤∑i∈Ifixk−fiyk.
By definition, the sequenceyk remains bounded. From the continuity of f, it follows that the sequencefiyk is also bounded for every i=1,…,n. Hence, we deduce from (24) that there is at least one index i0∈I such that(25)xi0k⟶∞,fi0xk⟶∞,k⟶+∞.
In what follows, we will show that(26)θxk=12∑i∈IfEFixk,ξ2+∑i∈Il∪IF+ϕ2EFixk,ξ,xi−li+∑i∈Iu∪IF−ϕ2−EFixk,ξ,ui−xi⟶+∞,k⟶+∞.
Now, we consider above three cases, respectively.Case 1
(i∈If): obviously, EFi0xk,ξ2=fi02xk⟶+∞,k⟶+∞.Case 2
(i∈Il∪IF+): by Lemma 1, we have(27)ϕ2EFi0xk,ξ,xi0k−li0≥23+22min2EFi0xk,ξ,xi0k−li0.Case 3
: ifminEFi0xk,ξ,ξ,xi0k−li0=EFi0xk,ξ, it follows from (25) that(28)ϕ2EFi0xk,ξ,xi0k−li0≥23+22EFi0xk,ξ2,=23+22fi02xk⟶+∞,k⟶+∞.Case 4
: ifminEFi0xk,ξ,ξ,xi0k−li0=xi0k−li0, then, by (25), we have(29)ϕ2EFi0xk,ξ,xi0k−li0≥23+22xi0k−li02⟶+∞,k⟶+∞.Case 5
(i∈Iu∪IF−): by Lemma 1, we have(30)ϕ2−EFi0xk,ξ,ui0−xi0k≥23+22min2−EFi0xk,ξ,ui0−xi0k.Case 6
: ifmin−EFi0xk,ξ,ui0−xi0k=−EFi0xk,ξ, it follows from (25) that(31)ϕ2−EFi0xk,ξ,ui0−xi0k≥23+22EFi0xk,ξ2,=23+22fi02xk⟶+∞,k⟶+∞.Case 7
: ifmin−EFi0xk,ξ,ui0−xi0k=ui0−xi0k, then, by (25), we have(32)ϕ2−EFi0xk,ξ,ui0−xi0k≥23+22xi0k−ui02⟶+∞,k⟶+∞.
Based on the above fact, we have(33)c¯≥θxk⟶+∞,k⟶+∞.
This is a contradiction, and hence,LθSc is bounded for any c≥0.
SinceS is nonempty, we see from Theorem 1 that problem (14) has at least one optimal solution when f is a uniform P-function on S.
Note that the expectation function of problem (14) is generally difficult to evaluate exactly. In what follows, we employ a Monte Carlo sampling method for numerical integration to address this question.


## 4. Monte Carlo Sampling Method and Convergence Analysis

We can view the generated sampleξ1,ξ2,…, as a sequence of random vectors, each having the same probability distribution as ξ. If the generated random vectors are stochastically independent of each other, we say that the sample is independent identically distributed (iid). With the generated sample ξi∈Ξ,i=1,2,…,N, we associate the sample average function:(34)θ^Nx≔12∑i∈Iff^iNx2+∑i∈Il∪IF+ϕ2f^iNx,xi−li+∑i∈Iu∪IF−ϕ2−f^iNx,ui−xi,∀x∈S,where f^iNx≔1/N∑j=1NFix,ξj, and consider the following approximation problem of (14):(35)min,θ^Nx,s.t.,x∈S.We refer (14) as the true problem and (35) as the Monte Carlo sampling average approximation problem.If the sample is iid, then the law of large numbers (LLN) holds pointwise. We say that the LLN holds, forf^Nx, pointwise if f^Nx converges w.p.1 to fx, as N⟶∞, for any fixed x∈S. See [12, 24, 25], for more details about the Monte Carlo sampling method.Theorem 2.
SupposeFx,ξ is continuous with respect to x for any ξ∈Ξ and EFx,ξ2<+∞,∀x∈S; if xN solves problem (35) for each N and x∗ is an accumulation point of xN as N tends to infinity, then x∗ is an optimal solution of the true problem (14) with probability one.Proof.
LetX⊆S be a compact set which w.p.1 contains a neighborhood of x∗. We first show that θ^Nx converges w.p.1 to θx for any x∈X. Recall (13) and (34), and we have(36)θ^Nx−θx=12∑i∈Il∪IF+ϕ2f^iNx,xi−li−ϕ2fix,xi−li+∑i∈Iu∪IF−ϕ2−f^iNx,ui−xi−ϕ2−fix,ui−xi+∑i∈Iff^iNx2−fix2.
Note that(37)∑i∈Il∪IF+ϕ2f^iNx,xi−li−ϕ2fix,xi−li=∑i∈Il∪IF+ϕf^iNx,xi−li−ϕfix,xi−liϕf^iNx,xi−li+ϕfix,xi−li=∑i∈Il∪IF+f^iNx−fix1+f^iNx+fixf^iNx2+xi−li2+fix2+xi−li2ϕf^iNx,xi−li+ϕfix,xi−li.
SinceX is a compact set and Fx,ξ is continuous with respect to x for any ξ∈Ξ, hence, for every x∈X, there exists a constant 0<M1<+∞ such that(38)1+f^iNx+fixf^iNx2+xi−li2+fix2+xi−li2ϕf^iNx,xi−li+ϕfix,xi−li≤M1,i∈Il∪IF+.
Combining (37) and (38), we have(39)∑i∈Il∪IF+ϕ2f^iNx,xi−li−ϕ2fix,xi−li≤M1∑i∈Il∪IF+f^iNx−fix.
Similarly, we can have that, for everyx∈X, there exists constant 0<M2<+∞ and 0<M3<+∞ such that(40)∑i∈Iu∪IF−ϕ2−f^iNx,ui−xi−ϕ2−fix,ui−xi≤M2∑i∈Iu∪IF−f^iNx−fix,(41)∑i∈Iff^iNx2−fix2≤M3∑i∈Iff^iNx−fix.
LetM=maxM1,M2,M3; by (39), (40), and (41), then(42)θ^Nx−θx≤M2∑i=1nf^iNx−fix.
f^Nx converges w.p.1 to fx, as N⟶∞, for any fixed x∈S; then, for any 2ε/3M>0, we have that w.p.1 for N is large enough:(43)∑i=1nf^iNx−fix≤2ε3M,w.p.1.
Obviously, for everyx∈X, we can have(44)θ^Nx−θx≤ε3,w.p.1.
Forx∗,xN∈X, we have from (44) that(45)θ^Nx∗−θx∗≤ε3,w.p.1,(46)θ^NxN−θxN≤ε3,w.p.1.
SincexN⟶x∗ as N⟶∞ and θ^Nx is continuous, we have(47)θ^NxN−θ^Nx∗≤ε3.
Combining (47) and (45), then there exists Nω>0 such that w.p.1 for N>Nω:(48)θ^NxN−θx∗≤θ^NxN−θ^Nx∗+θ^Nx∗−θx∗≤2ε3.
Note that(49)θxN−θx∗≤θxN−θ^NxN+θ^NxN−θx∗.
Combining (46) and (48), we have that w.p.1 for sufficiently large N:(50)θxN−θx∗≤ε3+2ε3=ε.
This implies that w.p.1 a global minimizer of (35) becomes a global minimizer of (14); hence, it is concluded.
Section5 will demonstrate the proposed approach.


## 5. Numerical Results

In this section, we used the notation and example in Wang et al. [26] to illustrate the model of BSVIP and the formulation. In our experiments, we used the command random in Matlab R2010a to generate pseudorandom sequences and employed fmincon to solve problem (35).Example 1.
(see [26]). Consider the stochastic variational inequality problem (2), in which ξ is uniformly distributed on Ξ=0,1,S=0,4×0,4×0,4, and F:ℝ3×Ξ⟶ℝ3 is given by(51)Fx,ξ=x1−ξx2+3−2ξ−ξx1+2x2+ξx3−2−ξξx2+3x3−3−ξ.
Namely,li=0andui=4,i=1,2,3. This problem has a solution x∗=0,1,1T for each ξ∈Ξ. The numerical results are shown in Table 1.Table 1
The computational results for Example1.
NxN=xN,1,xN,2,xN,3θ^NxN20(1.382e − 2, 1.086204, 1.103972)1.0861e − 2200(7.736e − 3, 0.990252, 1.068744)6.2177e − 31000(5.3081e − 4, 0.996378, 0.987653)3.4724e − 410000(8.2676e − 5, 1.000652, 1.000933)4.5232e − 520000(3.275e − 6, 1.000083, 0.999951)1.2911e − 6Example 2.
(see [26]). Consider the stochastic variational inequality problem (2), in which ξ is uniformly distributed on Ξ=0,1,S=0,4×0,4×0,4 and F:ℝ3×Ξ⟶ℝ3 is given by(52)Fx,ξ=x12−ξx2+3−2ξ−ξx1+2x22+ξx3−2−ξξx2+3x32−3−ξ.
Namely,li=0andui=4,i=1,2,3. It is easy to prove that the function EFx,ξ is strongly monotonous. So, this stochastic variational inequality problem has a unique solution x∗=0,1,1T for each ξ∈Ξ. The numerical results are shown in Table 2.
From the above analysis for Examples1-2, our preliminary numerical results for these examples indicate that the proposed method yields a reasonable and better solution of the stochastic variational inequality problem (2).Table 2
The computational results for Example2.
NxN=xN,1,xN,2,xN,3θ^NxN20(9.0123e − 3, 0.932573, 0.965513)9.4742e − 3200(1.7361e − 3, 1.009566, 1.010213)1.2275e − 31000(8.7204e − 5, 0.995087, 1.008562)7.8327e − 510000(1.3375e − 5, 1.999361, 1.000375)1.0253e − 520000(1.0017e − 6, 1.000027, 0.999987)1.0009e − 6


## 6. Conclusions

In this paper, the well-known Mote Carlo sampling method is applied to solve a class of box-constrained stochastic variational inequality problems. Firstly, a merit function derived from the FB-function is used and BSVIP is formulated as an optimization problem that minimizes this merit function. A sufficient condition is suggested for the existence of a solution to the optimization problem. Finally, this paper proposes a Monte Carlo sampling method for solving the problem. Under some moderate conditions, comprehensive convergence analysis is included as well.



---

*Source: 6618180-2021-06-02.xml*