# A Kind of Stochastic Eigenvalue Complementarity Problems

Authors: Ying-xiao Wang; Shou-qiang Du
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

With the development of computer science, computational electromagnetics have also been widely used. Electromagnetic phenomena are closely related to eigenvalue problems. On the other hand, in order to solve the uncertainty of input data, the stochastic eigenvalue complementarity problem, which is a general formulation for the eigenvalue complementarity problem, has aroused interest in research. So, in this paper, we propose a new kind of stochastic eigenvalue complementarity problem. We reformulate the given stochastic eigenvalue complementarity problem as a system of nonsmooth equations with nonnegative constraints. Then, a projected smoothing Newton method is presented to solve it. The global and local convergence properties of the given method for solving the proposed stochastic eigenvalue complementarity problem are also given. Finally, the related numerical results show that the proposed method is efficient.


---

## Body

## 1. Introduction

Computational electromagnetics is a science, which spans many subjects. It is an organic combination of mathematical theory, electromagnetic theory, and computer science. Electromagnetism is the classical dynamics theory. Electromagnetic phenomena are also closely related to eigenvalue problems and essentially can be deduced by Maxwell eigenvalue equation. The typical eigenvalue problem of the electromagnetic field is the resonant problem of the cavity and the guide wave problem of the waveguide. No matter the resonance problem of closed cavity, or the propagation problem of uniform guided waves, the wave equation is homogeneous without considering any source or field excitation process. For simplification, the conductors that constitute the cavity are idealized. So the field satisfies the homogeneous boundary condition on the boundary. There are many methods for solving computational electromagnetics, such as finite difference method and finite element method; see literatures [1–7] for details. In recent years, various methods in computational electromagnetics have been continuously improved. However, the input data are usually affected by many uncertain factors, such as the environment, which often results in the difference between calculation and measurement. So, we consider establishing a new model to solve some related electromagnetics problems. This is also the motivation of this paper.Among the optimization model, as we all know, the eigenvalue complementarity problem arises from several important applications in engineering and physics, such as the study of the resonance frequency and the stability of dynamic systems. In the last few years, the eigenvalue complementarity problem has drawn increasing attention, in many literature systems, such as [8–13] and the references therein. Among them, in [8], the authors study an eigenvalue complementarity problem and find its origins in the solution of a contract problem in mechanics. In [9], the eigenvalue complementarity problems with symmetric real matrices are considered. The authors transform this problem into a differentiable optimization program involving the Rayleigh quotient on a simplex and find its stationary point by the spectral projected gradient algorithm. In [10–13], many methods are proposed to solve the eigenvalue complementarity problems, such as Levenberg-Marquardt method and the derivative-free projection method. In [14], the stability of dynamic system is studied. The model of eigenvalue complementarity problem is established and the finite element method is used to solve the model.On the other hand, the stochastic linear complementarity problems, stochastic nonlinear complementarity problems, and some related stochastic optimization problems attract many researchers’ attention. And a lot of related theories and algorithms are given, such as [15–20] and their references. The expected value (EV) to reformulate the problems is usually used. Many methods are effectual for solving the stochastic complementarity problems. The global and local convergence of these methods are also given. Besides, the expected residual minimization (ERM) also can reformulate the stochastic complementarity problems, such as those in [21–23]. And the stochastic complementarity problems are eventually transformed into a certain minimal problem with constraints.In order to reduce the difference between calculation and measurement in computational electromagnetics, we consider getting the approximate distribution, by statistics, of a large number of experimental data and establishing a new stochastic model. So, in this paper, we propose a new kind of stochastic eigenvalue complementarity problem and the smoothing Newton method is also proposed to solve the given problem.Denote(Ω,F,P) as a probability space with Ω⊆Rn, where the probability distribution P is known. The new kind of stochastic eigenvalue complementarity problem is to find a scalar λ>0 and a vector x∈Rn∖{0}, such that(1)u=λBω-Aωx,x≥0,u≥0,xTu=0,where Aω and Bω are n×n symmetric positive definite stochastic matrices, Ω⊆Rn is the probability space, and ω∈Ω={ω1,ω2,…,ωm} is a random vector given in probability distribution P.The rest of this paper is organized as follows. In Section2, we reformulate problem (1) as a nonsmooth system firstly and give the related preliminary results. In Section 3, we present the projected smoothing Newton method and give the convergence results. In Section 4, some preliminary numerical results are proposed. The last section is the conclusion.


## 2. Preliminary Results

In this section, we give the reformulation of problem (1) and some related preliminaries.Firstly, we consider using the expected value model to solve (1) as the following equations.Denote(2)B¯=∑i=1mpiBωi,A¯=∑i=1mpiAωi,where pi=Pωi∈Ω≥0,i=1,…,m. So (1) is equivalent to the problem, which is defined as(3)x≥0,λB¯-A¯x≥0,xTλB¯-A¯x=0,(4)λBωi-Aωix≥0,i=1,…,m.As (3) is a linear complementarity problem, we can use the complementarity functions to transform it. The complementarity function has many different forms, such as the functions given in [24–29]. Among them, the Fischer-Burmeister (FB) function(5)ϕx=a2+b2-a+bis used widely. By the FB function transformation, (3) is equivalent to Φ(x)=0, where(6)Φx=ϕx1,λB¯-A¯1x⋮ϕxn,λB¯-A¯nx.In this paper, we consider using the smoothing Fischer-Burmeister function, which is differentiable at any point. This smoothing method, which is an important method to solve nonsmooth problems, has been widely used in recent years, such as in [30–33]. The idea of a smoothing method is to approximate nonsmooth functions by a sequence of smooth functions. There are many forms of smooth functions. Qi and Chen proposed a smooth approximation method for nonsmooth functions in [34]. By introducing the smoothing approximation function of Φμ:Rn→Rn, μ>0, smoothing methods are proposed by Chen et al. So far, a lot of progress has been made in the study of smoothing methods. In this paper, similar to [28], we use the smoothing Fischer-Burmeister function as(7)ϕμa,b=a2+b2+μ-a+b,where μ>0.Therefore, (3) is equivalent to Φμ(x)=0, where(8)Φμx=ϕx1,λB¯-A¯1x⋮ϕxn,λB¯-A¯nx.Then, (3) and (4) are equivalent to the following problem:(9)Hx,λ,μ,y=0,y≥0,λ>0,μ>0,where(10)Hx,λ,μ,y=ΦμxλBω1-Aω1x-y1⋮λBωm-Aωmx-ymxTx-12μ-1=0and y=(y1T,…,ymT)T∈Rmn. Let z=(xT,λ,μ,yT)T∈R(m+1)n+2 and define a merit function of (9) as f(z)=(1/2)H(z)2.If (1) has a solution, then solving (9) is equivalent to finding a global solution of the following optimization problem:(11)min⁡fzs.t.z∈Z,where Z={xT,λ,μ,yT)T∣x≥0,λ>0,μ>0,y≥0.Next, we give some related definitions and propositions.Definition 1.
Assume thatF:Rm→Rn is a locally Lipschitzian function. The B-subdifferential of F at x is (12)∂BFx=V∈Rm×n∣∃xk⊆DF:xk⟶x,F′xk⟶V,where DF is the differentiable points set and F′xk is the Jacobian of F at a point x∈Rn.
Here, for anyz=(xT,λ,μ,yT)T∈R(m+1)n+2, we have(13)∂BHz=VΦμVΦμ′VΦμ″0⋯0λBω1-Aω1Bω10-I⋯0⋮⋮⋮⋮⋮⋮λBωm-AωmBωm00⋯-I2x000000002μln⁡200,where VΦμ is the B-subdifferential of Φμ to x, VΦμ′ is the B-subdifferential of Φμ to λ, and VΦμ″ is the B-subdifferential of Φμ to μ. I is the n×n identity matrix and λ∈R+,μ∈R+.Definition 2 (see [18]).
F is said to be semismooth at x if(14)limV∈∂BFx+th′,h′→h,t↓0⁡Vh′exist for any h∈Rn.Definition 3 (see [18]).
F is said to be strongly semismooth at x if F is semismooth at x and for any V∈∂F(z+h),h→0,(15)limh→0,V∈∂BFx+h⁡Fx+h-Fx-Vhd2<∞.Definition 4 (see [35]).
F is said to be BD-regular at a point x if all the elements in ∂BF(x) are nonsingular.


## 3. Projected Smoothing Newton Method and Convergence Analysis

In this section, we present the projected smoothing Newton method and establish the global and local convergence results of the given method.Method 1 (projected smoothing Newton method).
Step  0. Choose σ∈(0,1),γ∈(0,1), μ¯∈R+ satisfying γμ¯<1, β∈(0,1).
Step  1. Compute gk=∇fzk. If gk<ϵ, stop.
Step  2. Compute dk by(16)Hzk+Vk+αIdk=βz¯,where Vk∈∂BH(zk),α≥0,z¯=(μ¯,0)∈R((m+1)n+2).
Step  3. Let m be the least nonnegative integer satisfying(17)fPZzk+βmdk≤1-2σ1-γμ¯βmfzk.Let zk+1=PZzk+βmdk.
Step  4. Set k≔k+1; return to Step  1.
In the following, we give the convergence results of Method1.Theorem 5.
Let{zk} be an infinite sequence generated by Method 1, and z∗ is an accumulation point of {zk}. If z∗ is a BD-regular solution of f(zk), z∗ is a solution of H(z)=0.Proof.
From Step  3 of Method1, we have f(zk+1)<f(zk), for all k≥0. Hence, f(zk) is monotonically decreasing. Since (zk)≥0(k≥0), there exists f∗, such that limk→∞fzk=f∗=f(z∗).
We supposef(z∗)>0. By (17), we have(18)fzk+1≤1-2σ1-γμ¯βmfzk.Take the limit of (18) on both sides. We get(19)fz∗≤1-2σ1-γμ¯βmfz∗.Then, we have (2σ(1-γμ¯)βm)f(z∗)≤0.
This is in contradiction toσ∈(0,1),γμ¯<1. Hence, H(z∗)=0. The proof is completed.Theorem 6.
Let{zk} be an infinite sequence generated by Method 1 and z∗ be a BD-regular solution; the rate of convergence is quadratic.Proof.
By Proposition2.1.1(b) in [36] and ∂BHzk being nonempty compact set at any point, we have that Vk∈∂BH(zk) is bounded in bounded set of points. So, there is a constant c1>0 and c3>0, such that, for all zk∈⋃(z∗,δ1), where δ1 is a sufficiently small positive constant, (Vk+αI)-1≤c1 and Vk≤c3 hold. Besides, by Definition 3, there exists a positive constant c2, such that(20)Hzk-Hz∗-Vkzk-z∗≤c2zk-z∗2,for all zk∈⋃(z∗,δ1), where δ2 is a sufficiently small positive constant. Moreover, there exists a positive constant c4, such that βz¯=c4zk-z82 hold.
Denoteδ=min⁡(δ1,δ2), for all zk∈⋃(z∗,δ). We have(21)Vk+αIzk+1-z∗=Vk+αIPZzk+βmdk-z∗=Vk+αIzk+dk-z∗+PZzk+βmdk-zk+dk=Vk+αIzk-z∗+dk+Vk+αIPZzk+βmdk-zk+dk=Vk+αIzk-z∗+Vk+αIdk+Vk+αIPZzk+βmdk-zk+dk=Vkzk-z∗+αzk-z∗+βz¯-Hzk+Vk+αIPZzk+βmdk-zk+dk.
SinceH is a locally Lipschitzian function, we premultiply this equation by (Vk+αI)-1 and take norms on both sides. So we get(22)zk+1-z∗≤Vk+αI-1Vkzk-z∗+αzk-z∗-Hzk+Vk+αIPZzk+βmdk-zk+dk+βz¯≤c1Hz∗-Hzk+Vkzk-z∗+αzk-z∗+PZzk+βmdk-zk+dk+βz¯≤c1c2+c3zk-z∗2+zk+dk-z∗+c4zk-z∗2≤c1c2+c-3+c4zk-z∗2+zk-z∗+Vk+αI-1βz¯-Hzk≤c1c2+c-3+c4zk-z∗2+c1c4zk-z∗2+c2zk-z∗2+c3zk-z∗2=2c1c2+c3+c4zk-z∗2=τzk-z∗2,where τ=2c1(c2+c3+c4). Therefore, similar to the proof of Theorem 2.3 in [37], we know that the rate of convergence is Q-quadratic. The proof is completed.


## 4. Numerical Results

In this section, we give some related numerical experiments. Because problem (1) is proposed for the first time, we give the following examples combined with the examples in literatures [1–7, 38]. In [7], the authors consider the complex shell inside the complicated enclosures of the electronic component and the high frequency radiation. The enclosure is assumed to be large compared to the wavelength. In this case, wavelengths have strong dependence on the properties of waves, including eigenvalues, eigenfunctions, scattering, and impedance matrices under small perturbations. The field fluctuations within the enclosure are described in a statistical sense using random matrix theory. In [38], in the microwave frequency band, the resonant circuit is usually realized by a resonant cavity. The cavity is filled with air or other mediums, and the electromagnetic oscillation in the cavity can be generated in the cavity by the way of probe and small hole. The waveguide cavity is formed by sealing the waveguides along the ends of the propagation direction with a conductor. Example 10 is converted from a two-dimensional rectangular resonant cavity problem in random system, which is proposed in [3, 38]. The program is written in Matlab 7.0. The parameters in Method 1 are taken as(23)δ=0.1,σ=0.1,μ0=1,μ¯=1,γ=0.1.The stopping rules for Method1 are ∇f(zk)≤1e-3 or kmax⁡=2×104.The initial points are randomly selected. In these tables of the numerical results,x∗ denotes the point of the minimization of f(z), λ denotes the eigenvalue, and DIM denotes the dimension of the problem.Example 7.
We consider problem (1), where(24)Aω=3+ω112+ω,Bω=2+ω111+ω,Ω=0,1,pi=Pωi∈Ω=0.5,i=1,2.
The numerical results of Example7 are given in Table 1 and the values of f(z) at each iteration are shown in Figure 1.Table 1
Numerical results of Example7.
λ x ∗ f z ∗ 1.5343 ( 0.3447,0.2414 ) T 0.4328 1.5280 ( 0.3827,0.3481 ) T 0.4193 1.5159 ( 0.4038,0.3369 ) T 0.3759 1.5472 ( 0.3267,0.1698 ) T 0.4327 1.4130 ( 0.6240,0.3500 ) T 0.1902Figure 1
Change curve of numerical value off(z) at every iteration.Example 8.
We consider problem (1), where(25)Aω=2+ω1111112+ω1111112+ω1111112+ω1111112+ω1111112+ω,Bω=3+ω1111113+ω1111113+ω1111113+ω1111113+ω1111113+ω,Ω=0,1,pi=Pωi∈Ω=0.5,i=1,2.
The numerical results of Example8 are shown in Table 2 and the values of f(z) at each iteration are shown in Figure 2.Table 2
Numerical results of Example8.
λ x ∗ f z ∗ 0.8996 ( 0.4856,0.5973,0.3077,0.4007,0.5974,0.5284 ) T 0.5782 0.9013 ( 0.6375,0.2817,0.3071,0.5837,0.3521,0.5110 ) T 0.2802 1.0000 ( 0.2128,0.2264,0.0828,0.2039,0.1965,0.0488 ) T 0.3125 0.9001 ( 0.5624,0.6134,0.4119,0.5258,0.5934,0.5930 ) T 0.4441 1.0000 ( 0.0096,0.0049,0.0080,0.0014,0.0042,0.0091 ) T 0.7799Figure 2
Change curve of numerical value off(z) at every iteration.Example 9.
We consider problem (1), where(26)Aω=2+ω1⋯112+ω⋯1⋮⋮⋱⋮11⋯2+ωn×n,Bω=3+ω1⋯113+ω⋯1⋮⋮⋱⋮11⋯3+ωn×n,Ω=0,1,pi=Pωi∈Ω=0.5,i=1,2.
Whenn takes different values, the numerical results of Example 9 are shown in Table 3 and the values of f(z) at each iteration are shown in Figure 3.
In the following, we consider the application of problem (1) in resonant cavity.
The calculation of the resonator can be described by the field equation and the boundary condition. With passivity in the cavity, the field vector satisfies the homogeneous Helmholtz equation. The vector potentialA and the vector potential Am are used, and an arbitrary field vector a is taken as a leading vector. We have(27)A=aΨ,Am=aΨ∗,where Ψ and Ψ∗ satisfy(28)∇2Ψ+k2Ψ=0,∇2Ψ∗+k2Ψ∗=0.
On the inner wallS of the cavity, the boundary conditions can often be converted to(29)ΨS=0,(30)Ψ∗S=0.
Besides, as we all know, generally speaking, microwave engineering can be divided into two categories: eigenvalue problems and noneigenvalue problems. The problems of the propagation of electromagnetic wave in the waveguide and the distribution of the midfield in the cavity are all eigenvalue problems. In [38], the eigenvalue problem here is as follows:(31)Ax=λx.
The reason for such a conclusion is that if the problem is abstracted, the Laplace operator∇2 acts on the function Ψ and the product of matrix A and vector x can be regarded as a linear operator L acting on an object f.
In [3], the author considers the eigenvalue problem as follows:(32)Kφ-λBφ=0,where K and B are matrices, φ is the unknown quantity, and λ is the eigenvalue.
In a random system, different loop lengths can produce different resonant frequencies, resulting in multiple modes. So we denote the resonant frequencies asλω, where ω is random.
In Example10 of this paper, we transform the above eigenvalue problem into the eigenvalue complementarity problem (1) and consider the form in the stochastic system. Next, we give Example 10.Table 3
Numerical results of Example9.
DIM λ f(z∗) 50 1.4000 0.7196 50 1.4000 0.7248 50 1.4000 0.6907 50 1.4000 0.6868 50 1.4000 0.6811 100 1.4000 0.9427 100 1.4000 0.8941 100 1.4000 0.8774 100 1.4000 0.8967 100 1.4000 0.8902 200 1.6595 1.7401 200 1.6610 1.6844 200 1.6626 1.5796 200 1.6584 1.8135 200 1.6611 1.7051Figure 3
Change curve of numerical value off(z) at every iteration, when DIM takes different values.
(a)
DIM = 50 (b)
DIM = 100 (c)
DIM = 200Example 10.
We consider problem (1), where(33)Bω=12ωω-ω2ωωω14ωω-ωω-ωω16ω-ωω2ω-ω-ω12ω-ωωωω-ω11ω,K=10231121212131111-11219111-1115,Ω=0,1,pi=Pωi∈Ω=0.5,i=1,2.
The numerical results of Example10 are shown in Table 4 and the values of f(z) at each iteration are shown in Figure 4.Table 4
Numerical results of Example10.
λ x ∗ f(z∗) 1.1038 (0.0027,0.0073,0.0074,0.0034,0.0050)T 0.5255 1.1097 (0.0010,0.0004,0.0055,0.0072,0.0075)T 0.5169 1.1004 (0.0013,0.0057,0.0047,0.0001,0.0033)T 0.5172 1.1718 (0.0065,0.0319,0.0129,0.0190,0.0034)T 0.6617 1.1417 (0.0321,0.0181,0.0307,0.0310,0.0236)T 0.9004Figure 4
Change curve of numerical value off(z) at every iteration.


## 5. Conclusion

The numerical experiments show that the projected smoothing Newton method is effective in solving the eigenvalue complementarity problems, even for high dimensional problems. After several experiments, we find that the value ofμ in smoothing function has a great influence on the numerical results. When μ=0.01, the numerical results are optimal. Using the smoothing Newton method, the objective function is polished. So the objective function is guided everywhere. Besides, the projected smoothing Newton method is Q-quadratic. In the application of the projected smoothing Newton method, we consider the two-dimension resonator problem. Therefore, using the new technology of optimization method to further improve the accuracy and efficiency of the given method, the applications in more complicated resonator problems are of great significance and practical value.



---

*Source: 7397592-2018-02-20.xml*