# A Continuous Trust-Region-Type Method for Solving Nonlinear Semidefinite Complementarity Problem

Authors: Ying Ji; Tienan Wang; Yijun Li
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We propose a new method to solve nonlinear
semidefinite complementarity problem by
combining a continuous method and a trust-region-type method. At every iteration, we need to calculate a second-order cone subproblem. We show the
well-definedness of the method. The global convergent result is
established.


---

## Body

## 1. Introduction

This paper deals with the semidefinite complementarity problem (SDCP) with respect to a mappingF:S→S, denoted by SDCP(F), to find an X∈S such that
(1)(X,F(X))∈S+×S+,〈X,F(X)〉=0,
where S⊂Rn×n is a set comprising those X∈Rn×n that are real symmetric. SDCP is the generalization of linear complementarity problems (LCPs) and semidefinite programs (SDPs) which has wide applications in engineering and economics [1]. The study of this problem can be dated back to the work of Shibata et al. [2]. Since then much attention has been attracted to SDCPs and various reformulations of SDCPs to minimization problem based on merit functions have been presented [3–6]. In general, there are two ways to derive the global convergence of an algorithm: trust-region methods and line search methods. The above methods proposed for solving SDCP(F) are all based on a line search strategy; methods based on trust-region technique are relatively fewer. Despite having been studied by many researchers [7, 8], trust-region methods are robust, can be applied to ill-conditioned problems, and have strong global convergence properties. Therefore, different from the above methods, we propose a new algorithm based on trust-region method to solve SDCPs.A functionF:S→S is said to be monotone if
(2)〈X-Y,F(X)-F(Y)〉≥0
for any X,Y∈S. An SDCP(F) is called a monotone SDCP(F) if the involved function is a monotone function. The Frobenius norm of a matrix X is defined by
(3)∥X∥F≔∥X∥≔〈X,X〉.
Let DF(X):S→S be a linear operator satisfying
(4)lim⁡ΔX→0∥F(X+ΔX)-F(X)-DF(X)ΔX∥∥ΔX∥=0;
then, F is said to be Fréchet differentiable at X and DF(X) is the Fréchet derivative of F at X. The function F is said to be differentiable if it is differentiable at each X∈S and to be continuously differentiable if also DF(X) is continuous at each X∈S. In this paper, we suppose that F:S→S is a continuously differentiable monotone function.Recently, there has been much interest in SDCP(F). A few methods have been developed to solve this problem, such as interior point methods, merit function methods, and noninterior point continuation/smoothing methods [3–5].Our new algorithm is based on the following smoothed Fischer-Burmeister function:(5)Φμ(X)=X+F(X)-X2+F(X)2+2μ2I,
where (μ,X,Y)∈R×S×S and I is the n×n identity matrix. This smoothing function was introduced by Kanzow [9] in the case of the NCP based on the Fischer-Burmeister function. Let
(6)Hμ(X)≔∥Φμ(X)∥2.
From Lemma 1 of [3], we know that if μ→0, then
(7)Hμ(X)⟶H0(X*)≔[X*-[X*-F(X*)]+,0]T,
where [X*-Y*]+ denotes the orthogonal projection of X*-Y* at S+, whereas by Lemma 2.1 of [4],
(8)H0(X*)=0⟺X*sloves  SDCP(F).
Thus, we can solve SDCP(F) by using the following approach: reformulate SDCP(F) as a system of nonsmooth equation H0(X)=0 and then approximate nonsmooth equations by parameterized smooth equations Hμ(X)=0; we solve the smooth equations at each iteration and make ∥Hμ(X)∥ decrease gradually by reducing the smoothing parameter μ to zero. In practice, however, it is usually impossible to solve the equation Hμ(X)=0 exactly for μ>0.In this paper, we present a continuous and approximate method to solve SDCP(F). At each iteration, the method solves a quadratic semidefinite program, which can be converted to a linear semidefinite program with a second-order cone constraint. A subproblem of this kind can be solved quite efficiently by using some recent software for semidefinite and second-order cone programs. The method is shown to be globally convergent under certain assumption.The rest of this paper is organized as follows. Section2 gives the algorithm and discusses the well-posedness for the algorithm; Section 3 analyzes the global convergence for the new algorithm; Section 4 presents the numerical results for the new algorithm; Section 5 concludes this paper.


## 2. The Algorithm

In this section, we will propose a smoothing trust-region-type algorithm for solving SDCP(F) and prove that the proposed algorithm is well-defined.Define(9)Qμk(ΔX)≔∥Φμk(Xk)+DΦμk(Xk)ΔX∥2.We begin with a formal statement of the algorithm which is in the spirit of [10–12].Algorithm 1 (continuous trust-region-type method).
(S0)
(Initialization) Choose0<ρ1<ρ2<1, 0<σ1<1<σ2, cmax⁡≥cmin⁡>0, c0∈[cmin⁡,cmax⁡], X0∈S+, Γ0:=(1+μ)∥Φ0(X0)∥, β0:=∥Φ0(X0)∥, κ:=2n, μ0:=((ε/(2Γ0κ))β02)2, and set k:=0.
(S1)
Find the solutionΔXk∈S of the subproblem
(10)min⁡ΔX∈Sn×n12ck〈ΔX,ΔX〉+Qμk(ΔX)s.t.Xk+ΔX⪰0.
If μk=0 and ΔXk=0, thenSTOP.
(S2)
Compute the ratio(11)rk≔Hμk(Xk)-Hμk(Xk+ΔXk)Hμk(Xk)-Qμk(ΔXk).
If rk≥ρ1, then the kth iteration is called successful, and we set Xk+1:=Xk+ΔXk; otherwise, the kth iteration is called unsuccessful, and we set Xk+1:=Xk.
(S3)
If(12)∥Φ0(Xk+1)∥≤max⁡{ηβk,ε-1∥Φ0(Xk+1)-Φμk(Xk+1)∥},
then set
(13)βk+1≔∥Φ0(Xk+1)∥
and choose μk+1 such that
(14)0<μk+1≤min⁡{(μ2Γ0κβk+12),μk4};
otherwise, let βk+1:=βk and μk+1:=μk.
(S4)
Updateck as follows.
Ifrk<ρ1, set ck+1≔σ2ck.
Ifrk∈[ρ1,ρ2), set ck+1≔mid(cmin⁡,ck,cmax⁡).
Ifrk≥ρ2, set ck+1≔mid(cmin⁡,σ1ck,cmax⁡).
(S5)
Setk←k+1, and go to (S1).
To verify that Algorithm1 is well-defined, we need the following properties of the smoothed Fischer-Burmeister function (5).Lemma 2 (see [1]).
Let(μ,X)∈R×Sn×n and Φμ(X) be defined by (5). Then(i)
ifμ>0, Φμ(X) is continuously differentiable at any X∈Sn×n;(ii)
for anyμ1,μ2>0 and X∈Sn×n, it follows that
(15)∥Φμ1(X)-Φμ2(X)∥≤2n|μ1-μ2|.Lemma 3.
LetXk be a given iterate and let ΔXk be the solution of the corresponding subproblem (10). Then
(16)Hμk(Xk)-Qμk(ΔXk)≥12ck〈ΔXk,ΔXk〉.Proof.
SinceXk≥0, the symmetric matrix ΔX:=0 is feasible for the subproblem (10). But ΔXk is a solution of this subproblem, so we obtain
(17)12ck〈ΔXk,ΔXk〉+Qμk(ΔXk)≤Qμk(0)=Hμk(Xk).
This proves our statement.The above lemma ensures that the denominator in the ratiork is always nonnegative. Note that this implies that the sequence {Hμk(Xk)} is monotonically nonincreasing. We next show that this denominator is equal to zero if and only if the termination criterion in step (S1) is satisfied. Hence, step (S2) is visited only if the denominator is positive, so that Algorithm 1 is well-defined.Lemma 4.
LetXk be a given iterate and ΔXk the solution of the corresponding subproblem. Then Hμk(Xk)-Qμk(ΔXk)=0 if and only if ΔXk=0.Proof.
First assume thatΔXk=0. Then Hμk(Xk)-Qμk(ΔXk)=0 since the definition of Qμk implies Qμk(0)=Hμk(Xk). Conversely, let Hμk(Xk)-Qμk(ΔXk)=0. Lemma 3 then implies 0=(1/2)ck〈ΔXk,ΔXk〉=(1/2)∥ΔXk∥2 and hence ΔXk=0.Next we have to justify our termination criterion in step (S1). To this end, we will show that this criterion is satisfied if and only if the current iterateXk is a stationary point of Hμk(X).Before we arrive at this result, we first take a closer look at subproblem. LetXk be a given iterate and let ΔXk be the unique solution of this subproblem. Since this subproblem is a convex program with a strictly feasible set, this problem is equivalent to its KKT conditions. In other words, ΔXk is a solution of subproblem if and only if there exist Lagrange multipliers Uk∈Sn×n such that the following KKT conditions hold:
(18)ckΔXk+DQμk(ΔXk)-Uk=0,Xk+ΔXk⪰0,Uk⪰0,〈Uk,Xk+ΔXk〉=0.
Now, if ΔXk=0 is the unique solution of this subproblem, then the system yields
(19)DQμk(0)-Uk=0,Xk⪰0,Uk⪰0,〈Uk,Xk〉=0.
However, these conditions are nothing, but the KKT conditions for the following problem are
(20)min⁡X∈S+∥Φμk(X)∥2.Summarizing these observations, we obtain the following result.Theorem 5.
Letμk=0. If ΔXk=0 is the (unique) solution of the subproblem for some ck>0, then Xk is a stationary point of the original problem. Conversely, if Xk is a stationary point of the original problem, then ΔXk=0 is the unique solution of subproblem for every ck>0.Proof.
The statements follow immediately from the preceding arguments.


## 3. Convergence Analysis

Throughout this section, we assume that Algorithm1 generates an infinite sequence {Xk}. Our aim is to establish a global convergence result for Algorithm 1. More precisely, we will show any accumulation point of {Xk} is a stationary point of the original problem.Lemma 6.
Let{Xk} be a sequence generated by Algorithm 1, and let μk and {Xk}k∈K be subsequences converging to 0 and some matrix X*, respectively, in such a way that {ck∥ΔXk∥}k∈K→0. Then X* is a stationary point of the original problem.Proof.
First note thatX* is symmetric positive semidefinite and hence feasible for original problem. Furthermore, since ck≥cmin⁡>0 for all k∈N, the assumption {ck∥ΔXk∥}k∈K→0 implies {∥ΔXk∥}k∈K→0. By continuity, we also have DHμk(Xk)→DH0(X*) as k→∞, k∈K. This together with system (18) implies that
(21)Uk=ckΔXk+DQμk(ΔXk),DQ0(0)=:U*
as k→∞, k∈K. Therefore, taking the limit k→∞ on the subsequence K in the KKT condition (18), we obtain
(22)DQ0(0)-U*=0,X*⪰0,U*⪰0,〈U*,X*〉=0.
Hence we conclude thatX* is a stationary point of the original problem.Another main step toward our global convergence result is contained in the following technical lemma.Lemma 7.
Let{Xk} be a sequence generated by Algorithm 1 and {Xk}k∈K a subsequence converging to some matrix X*. If X* is not a stationary point, then one has limsup⁡k→∞,k∈Kck<∞.Proof.
LetK-:={k-1∣k∈K}. Then we have {Xk+1}k∈K-→X*. We will show that limsup⁡k→∞,k∈K-ck+1<∞. Assume the contrary. Then, if necessary, we may suppose without loss of generality that
(23)lim⁡k→∞,k∈K-ck+1=∞.
The updating rule in step (S3) then implies that none of the iterations k∈K- with k sufficiently large is successful since otherwise we would have ck+1≤cmax⁡ for all these k∈K-. Hence, we have
(24)rk<ρ1
and Xk=Xk+1 for all k∈K- large enough. Since {Xk+1}k∈K-→X*, this implies {Xk}k∈K-→X*, too. Further, noticing that ck+1=σ2ck for all unsuccessful iterations, we also have
(25)lim⁡k→∞,k∈K-ck=∞
because of (23). We now want to show that
(26)rk⟶1ask⟶∞,k∈K-,
which would then lead to the desired contradiction to (24). To this end, we first note that
(27)liminf⁡k→∞,k∈K-ck∥ΔXk∥>0.
In fact, if ck∥ΔXk∥→0 on a subsequence, we would deduce from Lemma 6 that X* is a stationary point in contradiction to our assumption. Hence there is a constant γ>0 such that
(28)ck∥ΔXk∥≥γ,k∈K-.
By Lemma 2, this implies
(29)Hμk(Xk)-Qμk(ΔXk)≥12ck〈ΔXk,ΔXk〉≥12γ∥ΔXk∥
for all k∈K- sufficiently large.
We further note that{∥ΔXk∥}k∈K-→0. Otherwise, it would follow from (25) that ck∥Xk∥2→∞ on a suitable subsequence. This, in turn, would imply that the optimal value of the subproblem tends to infinity. However, this cannot be true since the feasible matrix ΔXk:=0 would give a smaller objective value. Hence we have {∥ΔXk∥}k∈K-→0.
Taking this into account and using{Xk}k∈K-→X* and the fact that F is continuously differentiable, we obtain through standard calculus arguments
(30)|Hμk(Xk)-Qμk(ΔXk)|=o(∥ΔXk∥)ask⟶∞,|Hμk(Xk)-Qμk(ΔXk)|=o(∥ΔXk∥)ffffflk∈K-.
Summarizing these observations, we get
(31)|rk-1|=|Hμk(0)-Hμk(Xk+ΔXk)Hμk(0)-Qμk(ΔXk)-1|=|Qμk(ΔXk)-Hμk(Xk+ΔXk)Hμk(0)-Qμk(ΔXk)|≤o(∥ΔXk∥)(1/2)γ∥ΔXk∥⟶0
as k→∞, k∈K-. This contradiction to (24) completes the proof.As a direct consequence of this lemma, we obtain the following result.Lemma 8.
Let{Xk} be sequence generated by Algorithm 1. Then there are infinitely many successful iterations.Proof.
If not, there would exist an indexk0∈N with rk<ρ1 and Xk=Xk0 for all k≥k0. This implies ck→∞ due to the updating rule in (S3). However, since Xk0 is not a stationary point and {Xk}→Xk0, we get a contradiction to Lemma 6.We are now in the position to prove the main convergence result for Algorithm1.Theorem 9.
Let{Xk} be a sequence generated by Algorithm 1. Then, any accumulation point of this sequence is a stationary point of the original problem.Proof.
LetX* be an accumulation point and {Xk}k∈K subsequence converging to X*. Since Xk=Xk+1 for all unsuccessful iterations k and since there are infinitely many successful iterations by Lemma 7, we may assume without loss of generality that all iterations k∈K are successful.
Assume thatX* is not a solution. Lemma 6 then implies
(32)limsup⁡k→∞,k∈Kck<∞.
Hence there is a constant γ>0 such that
(33)ck≥γ,k∈K.
Since each iteration k∈K is successful, we also have rl≥ρ1. Consequently, we obtain from Lemma 2(34)Hμk(Xk)-Hμk(Xk+ΔXk)≥ρ1(Hμk(Xk)-Qμk(ΔXk))≥12ρ1ck〈ΔXk,ΔXk〉≥12ρ1cmin⁡∥ΔXk∥2
for all k∈K. Since {Hμk(Xk)} is monotonically nonincreasing and bounded from below by, for example, H0(X*), we have Hμk(Xk)-Hmk(Xk+1)→0 as k→∞. Therefore, we obtain {ΔXk}k∈K→0 from (34). By (33), this also implies {ck∥Xk∥}k∈K→0. But then Lemma 6 shows that X* is a solution in contradiction to our assumption. This completes the proof.In the following, we will give one stronger global convergence result. Define the index set(35)ℵ≔{0}∪{{ηβk,μ-1∥Φ(Xk+1)-Φεk(Xk+1)∥}k∣∥Φ(Xk+1)∥{0}∪≤max⁡{ηβk,μ-1∥Φ(Xk+1)-Φεk(Xk+1)∥}}={k0=0<k1<k2<⋯}.Lemma 10.
IfF is a P0 function, then the sequence Xk generated by Algorithm 1 remains in the level set
(36)L0≔{X∈S∣H(X)≤(1+ε)2H(X0)}.Proof.
Letk be an arbitrary nonnegative integer, and let kj be the largest number in ℵ such that kj≤k, as ℵ is as defined from (35). It is easy to deduce from step 3 of Algorithm 1 that
(37)μk=μkj,βk=βkj,askj≤k<kj+1,∥Φμk(Xk)∥≤∥Φμkj(Xk)∥,askj≤k<kj+1.
Set
(38)Uj≔{X∈S∣∥Φμk(Xk)∥≤∥Φμkj(Xk)∥}.
As k is an arbitrary integer and xk∈Uj, it follows that Uj⊆L0.
Next, by induction, we will prove(39)Uj⊆L0,∀j≥0.
In view of Lemma2, we deduce that ∀X∈Uj,
(40)∥Φ(X)∥≤∥Φμkj(Xk)∥+κμkj≤∥Φμkj(Xkj)∥+κμkj≤∥Φ(Xkj)∥+2κμkj≤βkj+κμkj.
If j=0, then by (40) we have
(41)∥Φ(X)∥≤∥Φ(X0)∥+2κμ0≤(1+ε)∥Φ(X0)∥,∀X∈U0.
This proves U0∈L0.
SupposeUj-1⊆L0 for some j>0. Then, xkj-1∈L0 and hence βkj-1≤Γ0. Set
(42)ℵ1≔{k∈ℵ∣ηβk-1≥μ-1∥Φ(Xk)-Φμk-1(Xk)∥},ℵ2≔{k∈ℵ∣ηβk-1<μ-1∥Φ(Xk)-Φμk-1(Xk)∥}.
It follows from step 3 of Algorithm 1 and Lemma 2 that
(43)βkj≤ηβkj-1=ηβkj-1,ifkj∈ℵ1,
or
(44)βkj≤κμμkj-1=κμμkj-1≤12C0βkj-12≤12βkj-12,μkj-1=hlκμμkj-1≤12C0βkj-12ifkj∈ℵ2.
This implies that
(45)βkj≤δ3βkj-1,∀kj∈ℵ,
where δ3=max⁡{1/2,η}. Moreover, we have
(46)μkj≤14μkj-1=14μkj-1,∀kj∈ℵ.
From (45) and (46), we deduce that, for j>0,
(47)βkj≤δ3jβ0=δ3j∥Φ(X0)∥,μkj≤14jμ0≤ε24j+1(Γ0κ)2∥Φ(X0)∥4≤ε24j+1κ2∥Φ(X0)∥2.
Combining (47) with (40), we have
(48)∥Φ(X)∥≤δ3j∥Φ(X0)∥+ε2j∥Φ(X0)∥∥Φ(X)∥≤δ3j(1+ε)∥Φ(X0)∥,∀X∈Uj,
which implies Uj⊆L0. Hence (39) is proved and Lemma 10 is valid.It follows from Fischer [13] that if F is a P0 function or, more generally, an R0-function, then the level set L0 as defined in Lemma 10 is compact. Lemma 3 shows that the sequence {Hμk} is monotonically decreasing and converges.Theorem 11.
Assume thatF is a P0 function. Let {Xk} be a sequence generated by Algorithm 1. If there exists at least an accumulation point in the sequence {Xk}, then the index set ℵ defined by (35) is infinite:
(49)lim⁡k→∞μk=0,lim⁡k→∞Φ(Xk)=0,lim⁡k→∞Φμk(Xk)=0.Proof.
We first prove that setℵ is infinite. By contradiction, assume that ℵ is finite. Let k- be the largest number in ℵ. Then for all k≥k-, μk=μk- and βk=βk-. Denote
(50)μ-≔μk-,β-≔βk-,ϕ(X)≔Φ(X)-Φμ-(X).
Then ∀k>k-,
(51)∥Φ(Xk)∥>max⁡{ηβ-,μ-1∥ϕ(X)∥},(52)Φ(Xk)=ϕ(Xk)+Φε-(Xk).
From Theorem9, it follows that there exists at least an accumulation point X-∈L0 of {Xk} such that
(53)∇Hμ-(X-)=0.
Next, assume that subsequence{Xk}k∈ℵ^ converges to X-. In view of (53), we have {Φμ-(Xk)}k∈ℵ^→Φμ-(X-)=0 and hence there exists k^≥k- such that, for all k∈ℵ^ with k≥k^,
(54)∥Φμ-(Xk)∥≤(1-μ)ηβ-.
This together with (51) and (52) shows that, for all k∈ℵ^ with k≥k^,
(55)∥Φμ-(Xk)∥≤(1-ε)∥Φ(Xk)∥∥Φμ-(Xk)∥≤(1-ε)(∥Φμ-(Xk)∥+∥ϕ(Xk)∥);
that is,
(56)∥Φε-(Xk)∥<(ε-1-1)∥ϕ(Xk)∥,
which means
(57)∥Φ(Xk)∥≤∥Φε-(Xk)∥+∥ϕ(Xk)∥∥Φ(Xk)∥<ε-1∥ϕ(Xk)∥,ask∈ℵ^,k≥k^.
This contradicts (51). Hence the set ℵ is infinite.
Next,{μk}→0 follows immediately from the updating rule of μk and the fact that the set ℵ is infinite. Moreover, by the proof of Lemma 10, we deduce
(58)∥Φ(Xk)∥≤δ4j(1+ε)∥Φ(X0)∥,askj≤k<kj+1.
Because the set ℵ is infinite, it follows from Lemma 10 and (58) that
(59)lim⁡k→∞Φ(Xk)=0,lim⁡k→∞Φμk(Xk)=0.
This completes the proof.As a consequence of the above theorem, we get the following global convergence result.Corollary 12.
Assume thatF is a P0 function. Let {Xk} be a sequence generated by Algorithm 1. Then every accumulation point of the sequence {Xk} is a solution of NCP(F).


## 4. Numerical Experiments

### 4.1. The Reformulation for Subproblem

To test the numerical performance of Algorithm1, we implemented the method in MATLAB (Version 7.0) using the SDPT3-Solver (Version 3.0) for the corresponding subproblems. First, we will give the reformulation of the subproblem. In order to solve nonlinear semidefinite programs of the form (1) by Algorithm 1, we have to be able to deal with a subproblem given by
(60)min⁡ΔX∈Sn×n12ck〈ΔX,ΔX〉+Qμk(ΔX)s.t.Xk+ΔX≥0.
For this purpose, we would like to use the SDPT3-Solver (version 3.0) from [14]. This software is designed to solve linear semidefinite programs with cone constraints of the form
(61)min⁡∑j=1ns‍〈Cjs,Xjs〉+∑i=1nq‍(ciq)Txiq+(cl)xls.t.∑j=1ns‍(Ajs)Tsvec(Xjs)+∑i=1nq‍(Aiq)Txiq+(Al)Txl=b,Xjs∈S+sj×sj,∀j,xiq∈Kqqi∀i,xl∈R+nl,
where Xjs, Xjs are symmetric matrices of dimension sj; ciq, xiq are vectors in Rqi; S+sj×sj denotes the sj-dimension positive semidefinite cone defined by S+sj×sj:={X∈Ssj×sj:X⪰0}; Dqqi denotes the qi-dimensional second-order cone defined by Kqqi:={x=(x1,x2:qiT)T∈Rqi:x1≥∥x2:qi∥}; cl and xl are vectors in Rnl; Ajs are s-j×m matrices with s-j=sj(sj+1)/2; Aiq and Al are qi×m and l×m matrices, respectively; and svec is the operator defined by svec(X):=(X(1,1),2X(1,2),X(2,2),2X(1,3),2X(2,3),X(3,3),…)T∈Rn(n+1)/2 for any symmetric matrix X∈Sn×n.We now want to rewrite the problem (60) in the form of (61). To this end, we need to make some reformulations, which will be described step by step in the following.First, we drop the constant from the objective function without affecting the problem. Next, we introduce the auxiliary variableS∈Sn×n and set Xk+ΔX=S. Because ΔX needs only to be symmetric and not to be positive semidefinite, we set Δx=svec(ΔX) and write the problem in terms of Δx∈Rn- with n-:=n(n+1)/2. Then problem (60) is equivalent to
(62)min⁡(12ck+∥svec(DΦμk(Xk))∥2)∥Δx∥2+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.svec(Xk)+Δx=svec(S),Δx∈Rn-,s⪰0.
By introducing the second-order cone constraint ∥Δx∥t, the above problem can be further rewritten as
(63)min⁡(12ck+∥svec(DΦmuk(Xk))∥2)t2+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.svec(Xk)+Δx=svec(S),∥Δx∥≤t,Δx∈Rn-,S⪰0,t∈R.Unfortunately, the termt2 is not linear as required in (61). So we replace t2 by the new variable s≥0 and add the constraint t2≤s. But this constraint can be rewritten as the semidefinite constraint
(64)(stt1)⪰0.
Introducing once again an auxiliary variable, problem (62) and hence the original subproblem (60) are equivalent to
(65)min⁡(12ck+∥svec(DΦμk(Xk))∥2)s+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.-svec(Xk)=Δx-svec(S),∥Δx∥≤t,(stt1)-W=0,Δx∈Rn-,W⪰0,S⪰0,t∈R,s∈R+.
We write the equality constraint
(66)(stt1)-W=0,
in the svec-notation. Then we get
(67)min⁡(12ck+∥svec(DΦμk(Xk))∥2)t+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.-svec(Xk)=Δx-svec(S),∥Δx∥≤t,(st1)-svec(W)=0,Δx∈Rn-,W⪰0,S⪰0,t∈R,s∈R+.We are now in a position to give the explicit correspondence between the parameters, variables, and input data in our last problem formulation (67) and those in the SDPT3 standard form. The problem parameters are given by
(68)ns≔2,nq≔1,s1≔n,s2≔2,q1≔1+n-,l≔1+2m.
The variables are given by
(69)X1s≔S∈S+n×n,X2s≔W∈S2×2,x1q≔(t,ΔxT)T∈Kq1+n-,xl≔(s,ξT,ωT)T∈R1+2m.
The input datum in the objective function is given by
(70)C1s≔0∈Sn×n,C2s≔0∈S2×2,c1q≔(0,svec(Df(Xk))T)T∈R1+n-,cl≔(12ck,αke,0)∈R1+2m
with
(71)e=(1,…,1)T∈Rm.
Finally, the matrices A1s∈Rn-×(n-+3+m), A2s∈R3×(n-+3+m), A1q∈R(1+n-)×(n-+3+m), and Al∈R(1+m+m)×(n-+3+m) and the vector b∈Rn-+3+m are given by
(72)A1s=(-I|0|0),A2s=(0|-I|0),A1q=(0I|020000⋮⋮⋮000|0⋯0-svec(W)⋯-svec(w)),Al=(000|10000|0⋯0I-I),b=(-svec(Xk)T|00-1|g(Xk)T)T.
This is the desired reformulation.It may be worth mentioning that problem (60) can also be transformed as
(73)min⁡(12ck+∥svec(DΦμk(Xk))∥2)t+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.-svec(Xk)=Δx-svec(S),∥Δx∥≤t,(stt1)-W=0,Δx∈Rn-,W⪰0,S⪰0,t∈R,s∈R+.
Since the constraint ∥Δx∥2≤t is equivalent to
(74)(tΔxTΔxI),
problem (73) can further be reformulated as a linear semidefinite program that involves a semidefinite cone constraint instead of a second-order cone constraint. However, such a semidefinite representation is much more expensive in terms of memory requirement. Therefore, we adopted the reformulation (73) in our numerical experiments.


### 4.2. Numerical Results

We present some numerical tests using Algorithm1. All the codes are written in MATLAB 7.10. The tests are conducted on a DELL computer with Intel(R)Core(TM)i5-2400 processor (3.10 GHz) and 4.00 GB of memory on Windows 7.Consider the following nonlinear semidefinite complementarity problem withF(X):=XX-2X+I. It is obvious that the solution set of SDCP(F) is nonempty, since X=I is its one solution. The parameters in the algorithm can be presented as follows: ρ1 and ρ2 can be randomly generated from [0.1,0.5] and [0.6,1], respectively; σ1 and σ2 are randomly generated from [0.5,1] and [1.5,2], respectively; cmin⁡ and cmax⁡ are randomly generated from [0.01,1] and [500,1000], respectively; μ is randomly chosen from [10,20]; c0 is randomly generated from [cmin⁡,cmax⁡] and X0:=ATA, where A∈Rm×n with every entry being randomly generated from [0,1]. The stopping criterion is set as ∥Xk∥≤10-6 and μk≤10-6.For the purpose of comparison, we also solve this problem by the following descent algorithm based on the method proposed in [15].Algorithm 13 (decent direction method).
(S0)
(Initialization) Choose0<ρ<1, 0<α<1, X0∈S+, Γ0:=(1+μ)∥Φ0(X0)∥, β0:=∥Φ0(X0)∥, κ:=2n, μ0:=((ε/(2Γ0κ))β02)2, and set k:=0.
(S1)
Find the solutionΔXk∈S of the subproblem
(75)min⁡ΔX∈Sn×n12〈ΔX,ΔX〉+Qμk(ΔX)s.t.Xk+ΔX⪰0.
If ΔXk=0, thenSTOP.
(S2)
Computeαk=max⁡{1,α,α2,…} such that
(76)Hμk(Xk)≥Hμk(Xk+αkΔXk)+ραkQμk(ΔXk).
(S3)
Letxk+1:=Xk+αkΔXk. If
(77)∥Φ0(Xk+1)∥≤max⁡{ηβk,ε-1∥Φ0(Xk+1)-Φμk(Xk+1)∥},
then set
(78)βk+1:=∥Φ0(Xk+1)∥
and choose μk+1 such that
(79)0<μk+1≤min⁡{(μ2Γ0κβk+12),μk4};
otherwise, let βk+1:=βk and μk+1:=μk.
(S4)
Setk←k+1, and go to (S1).The above descent algorithm uses different NCP function from the one used in [15]. The parameters in this algorithm are set as follows: ρ can be randomly generated from [0.1,0.5]; α:=0.5; the starting point X0 and the stopping criteria are the same as Algorithm 1.We now solve this problem 40 times by Algorithms1 and 13, respectively, with the initial point X0 being randomly generated as above. Table 1 lists the numerical results for the applications of Algorithms 1 and 13. The average number of iterations and the average computational time (CPU time) are reported in Table 1. The results generally show that our method is efficient in solving this problem.Table 1
Numerical results.
n
m
Algorithm1
Algorithm13
Average iteration
Average CPU
Average iteration
Average CPU
5
2
5.4
0.03
15
0.04
10
5
13.54
2.22
20.57
3.88
10
10
12.88
2.25
19.27
4.00
15
5
15.97
1.02
55.92
10.23
15
10
20.74
1.05
34.27
15.17
20
5
30.11
10.22
100.68
25.55
20
10
37.35
11.05
157.39
30.11
25
5
41.99
24.54
217.22
44.67
25
10
47.02
30.22
348.04
61.22
25
15
45.11
27.66
477.93
110.84


## 4.1. The Reformulation for Subproblem

To test the numerical performance of Algorithm1, we implemented the method in MATLAB (Version 7.0) using the SDPT3-Solver (Version 3.0) for the corresponding subproblems. First, we will give the reformulation of the subproblem. In order to solve nonlinear semidefinite programs of the form (1) by Algorithm 1, we have to be able to deal with a subproblem given by
(60)min⁡ΔX∈Sn×n12ck〈ΔX,ΔX〉+Qμk(ΔX)s.t.Xk+ΔX≥0.
For this purpose, we would like to use the SDPT3-Solver (version 3.0) from [14]. This software is designed to solve linear semidefinite programs with cone constraints of the form
(61)min⁡∑j=1ns‍〈Cjs,Xjs〉+∑i=1nq‍(ciq)Txiq+(cl)xls.t.∑j=1ns‍(Ajs)Tsvec(Xjs)+∑i=1nq‍(Aiq)Txiq+(Al)Txl=b,Xjs∈S+sj×sj,∀j,xiq∈Kqqi∀i,xl∈R+nl,
where Xjs, Xjs are symmetric matrices of dimension sj; ciq, xiq are vectors in Rqi; S+sj×sj denotes the sj-dimension positive semidefinite cone defined by S+sj×sj:={X∈Ssj×sj:X⪰0}; Dqqi denotes the qi-dimensional second-order cone defined by Kqqi:={x=(x1,x2:qiT)T∈Rqi:x1≥∥x2:qi∥}; cl and xl are vectors in Rnl; Ajs are s-j×m matrices with s-j=sj(sj+1)/2; Aiq and Al are qi×m and l×m matrices, respectively; and svec is the operator defined by svec(X):=(X(1,1),2X(1,2),X(2,2),2X(1,3),2X(2,3),X(3,3),…)T∈Rn(n+1)/2 for any symmetric matrix X∈Sn×n.We now want to rewrite the problem (60) in the form of (61). To this end, we need to make some reformulations, which will be described step by step in the following.First, we drop the constant from the objective function without affecting the problem. Next, we introduce the auxiliary variableS∈Sn×n and set Xk+ΔX=S. Because ΔX needs only to be symmetric and not to be positive semidefinite, we set Δx=svec(ΔX) and write the problem in terms of Δx∈Rn- with n-:=n(n+1)/2. Then problem (60) is equivalent to
(62)min⁡(12ck+∥svec(DΦμk(Xk))∥2)∥Δx∥2+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.svec(Xk)+Δx=svec(S),Δx∈Rn-,s⪰0.
By introducing the second-order cone constraint ∥Δx∥t, the above problem can be further rewritten as
(63)min⁡(12ck+∥svec(DΦmuk(Xk))∥2)t2+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.svec(Xk)+Δx=svec(S),∥Δx∥≤t,Δx∈Rn-,S⪰0,t∈R.Unfortunately, the termt2 is not linear as required in (61). So we replace t2 by the new variable s≥0 and add the constraint t2≤s. But this constraint can be rewritten as the semidefinite constraint
(64)(stt1)⪰0.
Introducing once again an auxiliary variable, problem (62) and hence the original subproblem (60) are equivalent to
(65)min⁡(12ck+∥svec(DΦμk(Xk))∥2)s+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.-svec(Xk)=Δx-svec(S),∥Δx∥≤t,(stt1)-W=0,Δx∈Rn-,W⪰0,S⪰0,t∈R,s∈R+.
We write the equality constraint
(66)(stt1)-W=0,
in the svec-notation. Then we get
(67)min⁡(12ck+∥svec(DΦμk(Xk))∥2)t+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.-svec(Xk)=Δx-svec(S),∥Δx∥≤t,(st1)-svec(W)=0,Δx∈Rn-,W⪰0,S⪰0,t∈R,s∈R+.We are now in a position to give the explicit correspondence between the parameters, variables, and input data in our last problem formulation (67) and those in the SDPT3 standard form. The problem parameters are given by
(68)ns≔2,nq≔1,s1≔n,s2≔2,q1≔1+n-,l≔1+2m.
The variables are given by
(69)X1s≔S∈S+n×n,X2s≔W∈S2×2,x1q≔(t,ΔxT)T∈Kq1+n-,xl≔(s,ξT,ωT)T∈R1+2m.
The input datum in the objective function is given by
(70)C1s≔0∈Sn×n,C2s≔0∈S2×2,c1q≔(0,svec(Df(Xk))T)T∈R1+n-,cl≔(12ck,αke,0)∈R1+2m
with
(71)e=(1,…,1)T∈Rm.
Finally, the matrices A1s∈Rn-×(n-+3+m), A2s∈R3×(n-+3+m), A1q∈R(1+n-)×(n-+3+m), and Al∈R(1+m+m)×(n-+3+m) and the vector b∈Rn-+3+m are given by
(72)A1s=(-I|0|0),A2s=(0|-I|0),A1q=(0I|020000⋮⋮⋮000|0⋯0-svec(W)⋯-svec(w)),Al=(000|10000|0⋯0I-I),b=(-svec(Xk)T|00-1|g(Xk)T)T.
This is the desired reformulation.It may be worth mentioning that problem (60) can also be transformed as
(73)min⁡(12ck+∥svec(DΦμk(Xk))∥2)t+2svec(Φμk(Xk)DΦμk(Xk))TΔxs.t.-svec(Xk)=Δx-svec(S),∥Δx∥≤t,(stt1)-W=0,Δx∈Rn-,W⪰0,S⪰0,t∈R,s∈R+.
Since the constraint ∥Δx∥2≤t is equivalent to
(74)(tΔxTΔxI),
problem (73) can further be reformulated as a linear semidefinite program that involves a semidefinite cone constraint instead of a second-order cone constraint. However, such a semidefinite representation is much more expensive in terms of memory requirement. Therefore, we adopted the reformulation (73) in our numerical experiments.


## 4.2. Numerical Results

We present some numerical tests using Algorithm1. All the codes are written in MATLAB 7.10. The tests are conducted on a DELL computer with Intel(R)Core(TM)i5-2400 processor (3.10 GHz) and 4.00 GB of memory on Windows 7.Consider the following nonlinear semidefinite complementarity problem withF(X):=XX-2X+I. It is obvious that the solution set of SDCP(F) is nonempty, since X=I is its one solution. The parameters in the algorithm can be presented as follows: ρ1 and ρ2 can be randomly generated from [0.1,0.5] and [0.6,1], respectively; σ1 and σ2 are randomly generated from [0.5,1] and [1.5,2], respectively; cmin⁡ and cmax⁡ are randomly generated from [0.01,1] and [500,1000], respectively; μ is randomly chosen from [10,20]; c0 is randomly generated from [cmin⁡,cmax⁡] and X0:=ATA, where A∈Rm×n with every entry being randomly generated from [0,1]. The stopping criterion is set as ∥Xk∥≤10-6 and μk≤10-6.For the purpose of comparison, we also solve this problem by the following descent algorithm based on the method proposed in [15].Algorithm 13 (decent direction method).
(S0)
(Initialization) Choose0<ρ<1, 0<α<1, X0∈S+, Γ0:=(1+μ)∥Φ0(X0)∥, β0:=∥Φ0(X0)∥, κ:=2n, μ0:=((ε/(2Γ0κ))β02)2, and set k:=0.
(S1)
Find the solutionΔXk∈S of the subproblem
(75)min⁡ΔX∈Sn×n12〈ΔX,ΔX〉+Qμk(ΔX)s.t.Xk+ΔX⪰0.
If ΔXk=0, thenSTOP.
(S2)
Computeαk=max⁡{1,α,α2,…} such that
(76)Hμk(Xk)≥Hμk(Xk+αkΔXk)+ραkQμk(ΔXk).
(S3)
Letxk+1:=Xk+αkΔXk. If
(77)∥Φ0(Xk+1)∥≤max⁡{ηβk,ε-1∥Φ0(Xk+1)-Φμk(Xk+1)∥},
then set
(78)βk+1:=∥Φ0(Xk+1)∥
and choose μk+1 such that
(79)0<μk+1≤min⁡{(μ2Γ0κβk+12),μk4};
otherwise, let βk+1:=βk and μk+1:=μk.
(S4)
Setk←k+1, and go to (S1).The above descent algorithm uses different NCP function from the one used in [15]. The parameters in this algorithm are set as follows: ρ can be randomly generated from [0.1,0.5]; α:=0.5; the starting point X0 and the stopping criteria are the same as Algorithm 1.We now solve this problem 40 times by Algorithms1 and 13, respectively, with the initial point X0 being randomly generated as above. Table 1 lists the numerical results for the applications of Algorithms 1 and 13. The average number of iterations and the average computational time (CPU time) are reported in Table 1. The results generally show that our method is efficient in solving this problem.Table 1
Numerical results.
n
m
Algorithm1
Algorithm13
Average iteration
Average CPU
Average iteration
Average CPU
5
2
5.4
0.03
15
0.04
10
5
13.54
2.22
20.57
3.88
10
10
12.88
2.25
19.27
4.00
15
5
15.97
1.02
55.92
10.23
15
10
20.74
1.05
34.27
15.17
20
5
30.11
10.22
100.68
25.55
20
10
37.35
11.05
157.39
30.11
25
5
41.99
24.54
217.22
44.67
25
10
47.02
30.22
348.04
61.22
25
15
45.11
27.66
477.93
110.84


## 5. Conclusion

In this paper, we propose a trust-region method to solve nonlinear semidefinite complementarity problem. The well-posedness of the new method is proved and the global convergence is also presented. The numerical comparisons with the descent algorithm show the efficiency of the proposed method. For further study, we will discuss the convergent rate of the algorithm.



---

*Source: 589731-2014-05-07.xml*