# The Primal-Dual Active Set Method for a Class of Nonlinear Problems withT-Monotone Operators

Authors: Xiahui He; Peng Yang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The family of primal-dual active set methods is drawing more attention in scientific and engineering applications due to its effectiveness and robustness for variational inequality problems. In this work, we introduce and study a primal-dual active set method for the solution of the variational inequality problems withT-monotone operators. We show that the sequence generated by the proposed method globally and monotonously converges to the unique solution of the variational inequality problem. Moreover, the convergence rate of the proposed scheme is analyzed under the framework of the algebraic setting; i.e., the established convergence results show that the iteration number of the methods is bounded by the number of the unknowns. Finally, numerical results show that the efficiency can be achieved by the primal-dual active set method.


---

## Body

## 1. Introduction

The variational inequality problem associated withT-monotone operators has many applications, e.g., the diffusion problem involving Michaelis-Menten or second-order irreversible reactions; see, for example, [1–7] and the references therein for details. Hence, there is a growing interest in finding robust and efficient methods for solving this kind of complementarity problems, which reflects in an increasing number of proposals of numerical methods for its solution in recent years. It is well known that the projected relaxation method [8–10] is a popular solution technique for this class of complementarity problems. A great advantage of this approach is that it is easy to implement and can be convergent for both problems with symmetric and nonsymmetric operator. However, the convergence of this kind of methods depends crucially on the choice of the relaxation parameter. Another popular approach for the solution of variational inequality problems is the Schwarz algorithm [1, 7, 11–14], which is based on the framework of domain decomposition methods [15–17]. A major advantage of the Schwarz methods is amenable to implement and its convergence rate will not be deteriorated with the refinement of the mesh size when it is applied to the system arising from the discretization of partial differential equations (PDEs). The theory of monotone and global convergences for the classical Schwarz algorithms is also obtained. However, this kind of methods depends on the shape of the computational domain. To fix these issues, the family of active set strategies [18–24] can be used to solve variational inequality problems in an efficient way, which is the focus of this work.The active set method consists of two major steps: in the first phase, an index set is decomposed into active and inactive parts with respect to the solution vector, based on a criterion specifying a certain active set method; and then in the second phase, a reduced system associated with the active and inactive sets is solved. We briefly mention a few related publications that partially motivated our current work. In [22], Kanzow shows that the primal-dual active set algorithm is an efficient and accurate method for large-scale linear complementarity problems. In [25], Puterman and Brumelle analyze the convergence properties for the primal-dual active set method with continuous state and control spaces. They showed that the algorithm is equivalent to a Newton’s method, but under very restrictive assumptions which were not easily verifiable. In [24, 26], an active-set method with nonlinear elimination is proposed for the fully implicit simulation of two-phase flows. In the proposed algorithm, a variational inequality formulation of two-phase flow problems is used to avoid nonphysical undershoot or overshoot of the saturation fractions, and employ a class of active-set reduced-space algorithms to solve the resultant nonlinear complementarity system arising at each implicit time step. In particular, [18] builds the relationship between the primal-dual active set method and the semismooth Newton method and shows that the proposed method is a specific case of the Newton-type method under reasonable assumptions.In general, there are two ways to study convergence of primal-dual active set methods for solving variational inequality problems [18, 27]. The first one is to prove that the method produces a monotone sequence, in which the initial guess often starts from a super-solution or a lower-solution of the problem. Convergence theorems established in this way are often based on the assumption that the matrix from the linear system is an M-matrix, and the corresponding variational inequality is usually restricted to the unilateral case. The other way is to prove that the active set method generates a minimizing iterative sequence under some objective function, which is usually used to the variational inequality with the bilateral case. In the later case, the matrix is often supposed to be symmetric and positive definite. However, when the primal-dual active set method is used to solve the variational inequality problem with T-monotone operators, there is not a general convergence theorem, owing to the fact that it is difficult to find such a merit function. In this paper, we get the equivalent relation between the primal-dual active set method and Howard’s algorithm introduced in [9, 28, 29] for the variational inequality problem with T-monotone operators, then show the convergence theorem of Howard’s algorithm, and thus obtain the convergence theorem of the primal-dual active set method. Using the equivalence of the proposed algorithms for variational inequality problems with T-monotone operators, we give a simple proof for the global monotone convergence of the primal-dual active set method and conclude that the primal-dual active set method converges in no more than n+1 iterations, where n is the size of the solution vector. To the best of the authors’ knowledge, this is the first attempt to apply the primal-dual active set method for the variational inequality problem associated with T-monotone operators.The rest of the paper is organized as follows. In Section2, we present some notations and model problem and give some preliminaries which are used throughout the paper. In Section 3, the primal-dual active set method for the bilateral variational inequality problem is proposed. We establish the equivalence between Howard’s algorithm and the primal-dual active set method in Section 4. Finally, in Section 5 we report some numerical results for the proposed methods, and the paper is concluded in Section 6.


## 2. Preliminaries

In this section, we present some notations and the model problem and give some preliminaries that are used throughout the paper. First of all, we introduce some notations. LetN={1,2,…,n} denote an index set. For any index sets I,J⊆N, AIJ is denoted as the submatrix of a matrix A∈Rn×n that consists of aij(i∈I,j∈J), and rI is defined as the subvector of a vector r∈Rn consisting of ri(i∈I). Let K be a subset of Rn and f an operator from K to Rn. Let any element v∈K be expressed by v=v++v- with v+=max⁡{v,0}, v-=min⁡v,0. Then the notion of T-monotone is defined as follows.Definition 1.
A functionf is called T-monotone over the subspce K⊆Rn, if it satisfies(1)fv-fw,v-w+≥0,∀v,w∈K,where ·,· denotes the inner product.Definition 2.
A functionf is called strictly T-monotone over the subspce K⊆Rn, if, for all v,w∈K, f(v)-f(w),(v-w)+=0 is equivalent to (v-w)+=0.Based on this definition, in this study we consider variational inequality problems associated withT-monotone operators as follows.Problem 3.
LetK={v∈Rn:ϕ≤v≤ψ}, ϕ,ψ be given vectors in Rn satisfying ϕ≤ψ, and f is a continuous and strictly T-monotone operator. Then the variational inequality problem associated with T-monotone operators is defined by(2)findx∈K,suchthatfx,v-x≥0,∀v∈K.T-monotone operators are a kind of important operators, which include many linear and nonlinear elliptic operators; see, e.g., [1, 3, 7]. One of the advantages of T-monotone operators is that the relevant algorithms often possess monotone convergence property. We would like to point out that problem (2) has a unique solution; see, e.g., [27, 30]. Moreover, if all components of vector ϕ become -∞ and ψ become +∞, then problem (2) reduces to the system of nonlinear equations(3)fx=0.Hence, system (3) also has a unique solution when f is a continuous and strictly T-monotone operator. Moreover, the T-monotone operator has the following properties.Lemma 4.
LetI, J be the subsets of the index set N={1,2,…,n} that satisfies J=N∖I. For any vectors y,z∈K, if the subvector yI=zI with yJ≥zJ, and the function f is a continuous T-monotone operator over the subspace K, then fI(y)≤fI(z).Proof.
Let the index setI^ be defined as(4)I^=i∈I:fiy>fizand J^=N∖I^. Without loss of generality, the set I^ is denoted as {1,2,…,k} and J^={k+1,…,n}.
In the following, we prove it by contradiction based on the assumption thatI^ is not empty. The set w1 is defined as(5)w1=z1+δ1,…,zk+δk,zk+1,…,znand w2=y with δi being a positive integer for all i∈I^. If δi is small enough, then we have fi(w1)<fi(w2) for all i∈I^ by the continuity of the operator f. In addition to that, since we have (w1-w2)J^+=0, then(6)0≤fw1-fw2,w1-w2+=∑i=1kδifiw1-fiw2<0,which is a contradiction. Therefore, I^=∅ and we can get the conclusion fI(y)≤fI(z).Lemma 5.
Let the functionf be a continuous strictly T-monotone operator over K and I, J be the subsets of the index set N={1,2,…,n} that satisfies J=N∖I. For any vectors y,z∈K, if yI≤zI and fJ(y)≤fJ(z), then y≤z.Proof.
LetI^={i∈N:yi≤zi} and J^=N∖I^. Similarly, without loss of generality, let I^={1,2,…,k} and J^=N∖I^.
In the following, we also prove it by contradiction thatJ^ is not empty. By the definition of I^, we have that I⊂I^ and J^⊂J. Moreover, we have(7)fJ^y≤fJ^z,yJ^>zJ^,yI^≤zI^.Moreover, we know that(8)fJ^zI^,yJ^≤fJ^yI^,yJ^=fJ^y≤fJ^z=fJ^zI^,zJ^,where the first inequality comes from (7) and Lemma 4.
Letw1=(zI^,yJ^) and w2=z, by the definition of the T-monotone operator; then we obtain that(9)0≤fw1-fw2,w1-w2+=fJ^w1-fJ^w2,yJ^-zJ^+≤0,where the second inequality comes from fJ^(w1)≤fJ^(w2) and (yJ^-zJ^)+>0. Hence(10)fJ^w1-fJ^w2,yJ^-zJ^+=0,which means yJ^-zJ^≤0, since f is a strictly T-monotone operator. This is a contradiction to (7), and it means that y≤z.


## 3. Primal-Dual Active Set Method for the Variational Inequality Problem

In this section, we use the family of primal-dual active set methods [18–23] for solving the variational inequality problem (2) and then establish the equivalence between the primal-dual active set method and Howard’s algorithm [9, 28, 29]. The focus of this study is on the following equivalent problem of the model problem (2):(11)fx-λ=0,Bx,λ=0,where(12)Bx,λ≔max⁡min⁡λ,cx-ϕ,cx-ψ.Here, the max-operation or min-operation is understood componentwise; c>0 is a constant. The primal-dual active set method is based on using (12) as a prediction strategy. In the method, an index set is partitioned into active and inactive parts, based on a criterion specifying a certain active set method; i.e., given a current primal-dual pair (x,λ), the choice for the next active and inactive sets is given by(13)Jϕ=i∈N:λi-cx-ϕi>0,Jψ=i∈N:λi-cx-ψi<0,and I=N∖(Jϕ∪Jψ). Below, we present a high level description of the basic algorithm for a general problem in Algorithm 1.Algorithm 1:The primal-dual active set method.
Step1. Initialize x0 and λ0. Set k:=0. Step2. Determine the active and inactive sets by J ϕ k = i ∈ N : λ i k - c x k - ϕ i > 0, J ψ k = { i ∈ N : λ i k - c x k - ψ i < 0 } . and Ik=N∖(Jϕk∪Jψk). Step  3. Let xk+1 and λk+1 be the solution of the nonlinear system f ( x k + 1 ) - λ k + 1 = 0, x k + 1 = ϕ  o n   J ϕ k, x k + 1 = ψ  o n   J ψ k, λ k + 1 = 0  o n   I k . Step4.If k≥1 and xk=xk-1, then stop. Otherwise set k:=k+1 and return to Step 2.Remark 6.
The convergence theorem of the primal-dual active set method for the variational inequality problem with linear operator is based on a merit function, which is to prove that the proposed method generates a minimizing sequence for the merit function [23]. In this case, the matrix from the discretization of the linear operator is often supposed to be symmetric and positive definite. Theoretically, this condition number estimate cannot be applied immediately to the family of variational inequality problems with T-monotone operators, since these operators do not have these properties and we cannot use the technique of the merit function to get the convergence of Algorithm 1.In this study, we establish the equivalence between and the primal-dual active set method and Howard’s algorithm [9, 28, 29] for the solution of (2), and then we obtain the convergence theorem of Howard’s algorithm. The use of Howard’s algorithm is based on the following equivalent problem of (2):(14)findx∈Rn,suchthatmax⁡min⁡fx,x-ϕ,x-ψ=0.And then we define the following two functions F,G∈Rn by(15)Fx≔min⁡fx,x-ϕ,Gx≔max⁡Fx,x-ψ,respectively. Based on above notations, for ∀x∈Rn (14) is equivalent to(16)findx∈Rn,suchthatGx=0.In the following, we introduce the other equivalent formation of (2). Let A={0,1}n with α=(αi)i∈N∈A; then we can also define the function(17)Fx≔minα∈A⁡fαx=0,where(18)fαx≔fixifαi=0,x-ϕiifαi=1.Moreover, for every β∈A, the function Fβ(x):Rn→Rn is defined by(19)Fβx≔Fixifβi=0,x-ψiifβi=1.In a similar way, the function G(x)=maxβ∈A⁡Fβ(x), and then (14) is equivalent to the following problem:(20)findx∈Rn,suchthatmaxβ∈A⁡Fβx=0.Based on the above notations, we also present a high level description of Howard’s algorithm for solving (20) in Algorithm 2.Algorithm 2:Howard’s algorithm.
Step1. Initialize β0 in A={0,1}n, and set k:=0. Step2. Find xk∈Rn such that Fβk(xk)=0. If k≥1 and xk=xk-1, then stop. Otherwise go to Step 3. Step  3. For every i∈N, take β k + 1 ≔ 0 i f   F i β x k ≥ x k - ψ i , 1 i f   F i β x k < x k - ψ i . Set k:=k+1 and return to Step 2.


## 4. Convergence Results

In the section, we show that the primal-dual active set method is equivalent to Howard’s algorithm for the model problem (2), and then we obtain the convergence theorem of the primal-dual active set method.To begin, let us focus on Howard’s algorithm (i.e., Algorithm2). Note that, at each step, Howard’s algorithm satisfies(21)λk+1=fxk+1,λik+1=0,i∈Ik,xk+1-ϕi=0,i∈Jϕk,xk+1-ψi=0,i∈Jψk.Then we have the following notation:(22)F~xk+1=0≔fixk+1ifcxk-ψi≤fixk≤cxk-ϕi,cxk+1-ϕiiffixk>cxk-ϕi,cxk+1-ψiiffixk<cxk-ψi.In the following, Howard’s algorithm is used to solve the equivalent problem(23)max⁡min⁡fx,cx-ϕ,cx-ψ=0.For k≥0, if we set(24)βik+1=0,αik+1=0fori∈Ik,βik+1=0,αik+1=1fori∈Jϕk,βik+1=1,αik+1=0fori∈Jψk.Note that αk+1 and βk+1 are defined from the previous step xk, as introduced in Howard’s algorithm, and then we obtain F~=Fβk+1. Therefore, (22) is equivalent to(25)Fβk+1xk+1=0,and xk+1 is defined as in Howard’s algorithm applied to (23).Remark 7.
If we choose the initial guess withλ0=f(x0) in Algorithm 1, then the primal-dual active set method (i.e., Algorithm 1) and Howard’s algorithm (i.e., Algorithm 2) for the variational inequality problem (2) are equivalent.We now state the main convergence result for Howard’s algorithm; i.e., the iteration number of the method is bounded by the number of the unknowns. We start with the following preliminary results to show the monotone convergence of Algorithm2.Lemma 8.
Iff is a continuous strictly T-monotone operator, then the functions F, Fβ, and G are monotone operators with each β∈A={0,1}n.Proof.
Let the vectorsx∈Rn,y∈Rn such that F(x)≤F(y), and let αx∈A be a minimizer that is associated with the function F(x) defined as follows:(26)αx≔0iffix≤x-ϕi1iffix>x-ϕifori∈N.Based on this notation, we obtain(27)fαxx=Fx≤Fy=minα∈A⁡fαy≤fαxy.Hence, we conclude that the inequality x≤y holds, owing to the fact that fαx is a monotone operator.
In what follows, we show that the functionFβ is a monotone operator. Let x,y be in Rn such that Fβ(x)≤Fβ(y). On the one hand, if βi=0 for some i∈N, then we have Fiβ(x)≤Fiβ(y) and therefore(28)fiαxx≤fiαyy≤fiαxy.On the other hand, if βi=1 for some i∈N, then we have (x-ψ)i≤(y-ψ)i.
Letα¯∈A be defined as follows:(29)α¯≔0ifβi=0,αix=01ifotherwisefori∈N.Then we can obtain that if βi=0 and αix=0 then(30)fiα¯x=fiαxx≤fiαxy=fiα¯y,and otherwise,(31)fiα¯x=x-ϕi≤y-ϕi=fiα¯y.Hence, fiα¯(x)≤fiα¯(y), and we also conclude that the inequality x≤y holds by using the monotonicity of the function fα¯.
In a similar way, we show that the functionG is a monotone operator. We assume G(x)≤G(y). By (20), let βy be the index such that G(y)=Fβy(y). Then we have(32)Fβyx≤Gy=Fβyy.Hence, the inequality x≤y is obtained, based on the fact that the function Fβy is a monotone operator.Theorem 9.
Letf be a continuous strictly T-monotone operator and x∗∈Rn be the unique solution of (2). Then the sequence {xk} given by Algorithm 2 satisfies(a)
xk≥xk+1 for all k≥0 and xk≤ψ for all k≥1;(b)
xk→x∗ at most n+1 iterations.Proof.
The uniqueness of the solution is obtained by the monotonocity of the functionG.(a)
Note that(33)Fβk+1xk+1=0=Fβkxk≤Gxk=Fβk+1xk.Then, we havexk≥xk+1, since the function Fβk+1 is monotone by Lemma 8. In the following we prove that(34)xk≤ψforallk≥1.In fact, we need only to prove thatx1-ψ≤0. If βi1=0, then we have (x1-ψ)i=0 by definition of x1. On the other hand if βi1=1, then we obtain Fi(x0)≥(x-ψ)i. Moreover, either Fi(x0) or (x0-ψ)i is zero by definition of x0. Hence (x0-ψ)i≤0 and (x1-ψ)i≤0.(b)
We first show that the sequence{βk}k≥1 is decreasing in A={0,1}n. In fact, if βik=0 for some k≥1, then Fi(xk)=0. Using (xk-ψ)i≤0 we have that βik+1=0. Since {βk}k≥1 is decreasing, the set Ik:={i∈N:Fi(xk)≥(xk-ψ)i} is thus increasing for k≥0. On the other hand, since card⁡(Ik)≤n, there exists a first index k∈N such that Ik=Ik+1, and we have αk=αk+1. Moreover, G(xk+1)=Fαk+2(xk+1)=Fαk+1(xk+1)=0, and thus xk+1 is the solution for some k≤n. This makes at most n+1 iterations.


## 5. Numerical Experiments

In the numerical experiments, we focus on investigating the qualitative properties of the solution algorithms and comparing them to some other known algorithms. In this test, we consider a problem defined in the unit squareΩ=(0,1)×(0,1) as follows: find u∗∈K such that(35)-△u+gu,x,y>0,u=ϕ,-△u+gu,x,y<0,u=ψ,-△u+gu,x,y=0,ϕ≤u≤ψ,where K={u∈H01(Ω):ϕ≤u≤ψa.e.inΩ}, ϕ, and ψ are given functions and g(u,x,y)=u/(1+u)+10x+y-8. We discretize the model problem by using the standard five-point difference scheme with a constant mesh step size: hx=hy=1/(m+1), where m denotes the number of mesh sizes in x- and y-directions (N=m2 is the total number of unknowns). Then the discretized problem (35) belongs to a variational inequality problem with a T-monotone operator [7].In the following tests, we choose the bound obstaclesϕ=0 and ψ=1 for the variational inequality problem and conduct the following experiments:(1)
Compare the primal-dual active set method with project successive overrelaxation (PSOR) method [8–10] and the classical additive Schwarz method [7, 12, 13].(2)
Compare different solution algorithms from the point of view of iteration numbers and the total computing time for the variational inequality problem withT-monotone operators.Throughout this section, for the construction of the classical overlapping additive Schwarz method as introduced in [7, 31], denoted by Schwarz, we partition the computational domian N=N1∪N2 into two equal parts with the overlapping size O(1/10), and the corresponding subproblems are solved by PSOR with the relaxation parameter ω=1.8. In the proposed primal-dual active set method, i.e., Algorithm 1, the corresponding nonlinear systems are solved by nonlinear Gauss-Seidel method. We mainly consider the effect of dimension (denoted by N) on the performance of each algorithm. The initial guess for all the methods is chosen as u0=8A-1(1,…,1)T with A being the coefficient matrix of -△. The tolerance of the three methods is chosen to be equal to 10-6 in the ·2-norm for both inner and outer iterations.In the experiment, we consider PSOR, the classical additive Schwarz method, and Algorithm1 for the numerical solution of (35). Algorithm 1 for c=1 is not convergent for (35), where similar results were discussed in [23]. Hence, we choose c=1000 in Algorithm 1. In the test, we mainly focus on the effect of the computational mesh size (denoted by N) to the performance of each algorithm, as listed in Tables 1 and 2, respectively. Below we list the observations made from the results.(1)
From Table1, we can see that the number of iterations for Schwarz and PSOR grows much with the increase of the number of the computational mesh sizes, while the number of iterations for Algorithm 1 does not grow much with the increase of the number of mesh sizes, which implies that Algorithm 1 is insensitive to the number of mesh sizes. Hence, we can conclude that the performance of the primal-dual active set method is better than the Schwarz or PSOR method in terms of the iteration number.(2)
On the other hand, if we focus on the total computing time, the performance of PSOR is better than that of the Schwarz method or the primal-dual active set method, as shown in Table2. Since the operator in the test problem (35) is semilinear, i.e., almost linear, the PSOR method behaves much better than it usually does, while the Schwarz method and the primal-dual active set method spend much time to solve the related subproblem at each outer iteration step. As a result, the Schwarz method and the primal-dual active set method performed not so good as one expected in terms of computing time, as shown in our numerical results. It is consonant with [31], in which the PSOR and Schwarz methods are used to solve linear complementarity problem.Table 1
Comparison of iteration numbers.
N PSOR Schwarz Algorithm1 100 35 16 8 400 42 23 12 900 63 39 15 1600 120 60 19 2500 187 85 23 3600 267 116 27Table 2
Comparison of the computing time in seconds.
N PSOR Schwarz Algorithm1 100 0.015 0.031 0.015 400 0.5 0.703 0.765 900 4.281 5.796 6.484 1600 23.984 29.312 33.312 2500 89.5 106.89 126.359 3600 289.046 334.64 406.734


## 6. Concluding Remarks

In this work, we have developed the family of primal-dual active set methods to solve the discrete nonlinear systems, arising from the variational inequality problems withT-monotone operators. We build the equivalent relation between the primal-dual active set method and Howard’s algorithm, then show the convergence theorem of Howard’s algorithm, and thus obtain the convergence theorem of the primal-dual active set method. The numerical experiments confirm the advantage of the primal-dual active set method, when compared to the traditional PSOR and Schwarz methods.



---

*Source: 2912301-2019-03-17.xml*