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

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.

1. Introduction

This paper deals with the semidefinite complementarity problem (SDCP) with respect to a mapping 𝐹∶𝑆→𝑆, denoted by SDCP(𝐹), to find an 𝑋∈𝑆 such that (𝑋,𝐹(𝑋))∈𝑆₊×𝑆₊,⟨𝑋,𝐹(𝑋)⟩=0, (1) where 𝑆⊂𝑅^𝑛×𝑛 is a set comprising those 𝑋∈𝑅^𝑛×𝑛 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(𝐹) 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 function 𝐹∶𝑆→𝑆 is said to be monotone if ⟨𝑋−𝑌,𝐹(𝑋)−𝐹(𝑌)⟩≥0 (2) for any 𝑋,𝑌∈𝑆. An SDCP(𝐹) is called a monotone SDCP(𝐹) if the involved function is a monotone function. The Frobenius norm of a matrix 𝑋 is defined by ‖𝑋‖_𝐹∶=‖𝑋‖∶=⟨𝑋,𝑋⟩. (3) Let 𝐷𝐹(𝑋)∶𝑆→𝑆 be a linear operator satisfying lim_Δ𝑋→0‖𝐹(𝑋+Δ𝑋)−𝐹(𝑋)−𝐷𝐹(𝑋)Δ𝑋‖‖Δ𝑋‖=0; (4) then, 𝐹 is said to be Fréchet differentiable at 𝑋 and 𝐷𝐹(𝑋) is the Fréchet derivative of 𝐹 at 𝑋. The function 𝐹 is said to be differentiable if it is differentiable at each 𝑋∈𝑆 and to be continuously differentiable if also 𝐷𝐹(𝑋) is continuous at each 𝑋∈𝑆. In this paper, we suppose that 𝐹∶𝑆→𝑆 is a continuously differentiable monotone function.

Recently, there has been much interest in SDCP(𝐹). 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: Φ_𝜇(𝑋)=𝑋+𝐹(𝑋)−𝑋²+𝐹(𝑋)²+2𝜇²𝐼, (5) where (𝜇,𝑋,𝑌)∈𝑅×𝑆×𝑆 and 𝐼 is the 𝑛×𝑛 identity matrix. This smoothing function was introduced by Kanzow [9] in the case of the NCP based on the Fischer-Burmeister function. Let 𝐻_𝜇(𝑋)∶=Φ_𝜇(𝑋)². (6) From Lemma 1 of [3], we know that if 𝜇→0, then 𝐻_𝜇(𝑋)⟶𝐻₀𝑋^∗∶=𝑋^∗−𝑋^∗−𝐹𝑋^∗₊,0^𝑇, (7) where [𝑋^∗−𝑌^∗]₊ denotes the orthogonal projection of 𝑋^∗−𝑌^∗ at 𝑆₊, whereas by Lemma 2.1 of [4], 𝐻₀𝑋^∗=0⟺𝑋^∗slovesSDCP(𝐹). (8) Thus, we can solve SDCP(𝐹) by using the following approach: reformulate SDCP(𝐹) as a system of nonsmooth equation 𝐻₀(𝑋)=0 and then approximate nonsmooth equations by parameterized smooth equations 𝐻_𝜇(𝑋)=0; we solve the smooth equations at each iteration and make ‖𝐻_𝜇(𝑋)‖ decrease gradually by reducing the smoothing parameter 𝜇 to zero. In practice, however, it is usually impossible to solve the equation 𝐻_𝜇(𝑋)=0 exactly for 𝜇>0.

In this paper, we present a continuous and approximate method to solve SDCP(𝐹). 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. Section 2 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(𝐹) and prove that the proposed algorithm is well-defined.

Define 𝑄_𝜇𝑘(Δ𝑋)∶=Φ_𝜇𝑘𝑋_𝑘+𝐷Φ_𝜇𝑘𝑋_𝑘Δ𝑋². (9)

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) Choose 0<𝜌₁<𝜌₂<1, 0<𝜎₁<1<𝜎₂, 𝑐_max≥𝑐_min>0, 𝑐₀∈[𝑐_min,𝑐_max], 𝑋⁰∈𝑆₊, Γ₀∶=(1+𝜇)‖Φ₀(𝑋⁰)‖, 𝛽₀∶=‖Φ₀(𝑋⁰)‖, 𝜅∶=√2𝑛, 𝜇₀∶=((𝜀/(2Γ₀𝜅))𝛽²₀)², and set 𝑘∶=0. (S1)Find the solution Δ𝑋^𝑘∈𝑆 of the subproblem min_{Δ𝑋∈𝑆𝑛×𝑛}12𝑐_𝑘⟨Δ𝑋,Δ𝑋⟩+𝑄_𝜇𝑘(Δ𝑋)s.t.𝑋^𝑘+Δ𝑋⪰0. (10) If 𝜇_𝑘=0 and Δ𝑋^𝑘=0, then STOP. (S2)Compute the ratio 𝑟_𝑘𝐻∶=^𝜇𝑘𝑋^𝑘−𝐻_𝜇𝑘𝑋^𝑘+Δ𝑋^𝑘𝐻_𝜇𝑘𝑋^𝑘−𝑄_𝜇𝑘Δ𝑋^𝑘. (11) If 𝑟_𝑘≥𝜌₁, then the 𝑘th iteration is called successful, and we set 𝑋^𝑘+1∶=𝑋^𝑘+Δ𝑋^𝑘; otherwise, the 𝑘th iteration is called unsuccessful, and we set 𝑋^𝑘+1∶=𝑋^𝑘. (S3)If Φ₀𝑋^𝑘+1≤max𝜂𝛽_𝑘,𝜀⁻¹Φ₀𝑋^𝑘+1−Φ_𝜇𝑘𝑋^𝑘+1, (12) then set 𝛽_𝑘+1∶=Φ₀𝑋^𝑘+1 (13) and choose 𝜇_𝑘+1 such that 0<𝜇_𝑘+1𝜇≤min2Γ₀𝜅𝛽²_𝑘+1,𝜇_𝑘4; (14) otherwise, let 𝛽_𝑘+1∶=𝛽_𝑘 and 𝜇_𝑘+1∶=𝜇_𝑘. (S4)Update 𝑐_𝑘 as follows. If 𝑟_𝑘<𝜌₁, set 𝑐_𝑘+1∶=𝜎₂𝑐_𝑘. If 𝑟_𝑘∈[𝜌₁,𝜌₂), set 𝑐_𝑘+1∶=mid(𝑐_min,𝑐_𝑘,𝑐_max). If 𝑟_𝑘≥𝜌₂, set 𝑐_𝑘+1∶=mid(𝑐_min,𝜎₁𝑐_𝑘,𝑐_max). (S5)Set 𝑘←𝑘+1, and go to (S1).
To verify that Algorithm 1 is well-defined, we need the following properties of the smoothed Fischer-Burmeister function (5).

Lemma 2 (see [1]). Let (𝜇,𝑋)∈𝑅×𝑆^𝑛×𝑛 and Φ_𝜇(𝑋) be defined by (5). Then(i)if 𝜇>0, Φ_𝜇(𝑋) is continuously differentiable at any 𝑋∈𝑆^𝑛×𝑛;(ii)for any 𝜇₁,𝜇₂>0 and 𝑋∈𝑆^𝑛×𝑛, it follows that Φ_𝜇1(𝑋)−Φ_𝜇2(𝑋)≤√2𝑛𝜇₁−𝜇₂. (15)

Lemma 3. Let 𝑋^𝑘 be a given iterate and let Δ𝑋^𝑘 be the solution of the corresponding subproblem (10). Then 𝐻_𝜇𝑘𝑋^𝑘−𝑄_𝜇𝑘Δ𝑋^𝑘≥12𝑐_𝑘Δ𝑋^𝑘,Δ𝑋^𝑘. (16)

Proof. Since 𝑋^𝑘≥0, the symmetric matrix Δ𝑋∶=0 is feasible for the subproblem (10). But Δ𝑋^𝑘 is a solution of this subproblem, so we obtain 12𝑐_𝑘Δ𝑋^𝑘,Δ𝑋^𝑘+𝑄_𝜇𝑘Δ𝑋^𝑘≤𝑄_𝜇𝑘(0)=𝐻_𝜇𝑘𝑋^𝑘. (17) This proves our statement.

The above lemma ensures that the denominator in the ratio 𝑟_𝑘 is always nonnegative. Note that this implies that the sequence {𝐻_𝜇𝑘(𝑋^𝑘)} 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. Let 𝑋^𝑘 be a given iterate and Δ𝑋^𝑘 the solution of the corresponding subproblem. Then 𝐻_𝜇𝑘(𝑋^𝑘)−𝑄_𝜇𝑘(Δ𝑋^𝑘)=0 if and only if Δ𝑋^𝑘=0.

Proof. First assume that Δ𝑋^𝑘=0. Then 𝐻_𝜇𝑘(𝑋^𝑘)−𝑄_𝜇𝑘(Δ𝑋^𝑘)=0 since the definition of 𝑄_𝜇𝑘 implies 𝑄_𝜇𝑘(0)=𝐻_𝜇𝑘(𝑋^𝑘). Conversely, let 𝐻_𝜇𝑘(𝑋^𝑘)−𝑄_𝜇𝑘(Δ𝑋^𝑘)=0. Lemma 3 then implies 0=(1/2)𝑐_𝑘⟨Δ𝑋^𝑘,Δ𝑋^𝑘⟩=(1/2)‖Δ𝑋^𝑘‖² and hence Δ𝑋^𝑘=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 iterate 𝑋^𝑘 is a stationary point of 𝐻_𝜇𝑘(𝑋).

Before we arrive at this result, we first take a closer look at subproblem. Let 𝑋^𝑘 be a given iterate and let Δ𝑋^𝑘 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, Δ𝑋^𝑘 is a solution of subproblem if and only if there exist Lagrange multipliers 𝑈^𝑘∈𝑆^𝑛×𝑛 such that the following KKT conditions hold: 𝑐_𝑘Δ𝑋^𝑘+𝐷𝑄_𝜇𝑘Δ𝑋^𝑘−𝑈^𝑘𝑋=0,_𝑘+Δ𝑋^𝑘⪰0,𝑈^𝑘𝑈⪰0,_𝑘,𝑋^𝑘+Δ𝑋^𝑘=0. (18) Now, if Δ𝑋^𝑘=0 is the unique solution of this subproblem, then the system yields 𝐷𝑄_𝜇𝑘(0)−𝑈^𝑘𝑋=0,_𝑘⪰0,𝑈^𝑘𝑈⪰0,_𝑘,𝑋^𝑘=0. (19) However, these conditions are nothing, but the KKT conditions for the following problem are min_𝑋∈𝑆+Φ_𝜇𝑘(𝑋)². (20)

Summarizing these observations, we obtain the following result.

Theorem 5. Let 𝜇_𝑘=0. If Δ𝑋^𝑘=0 is the (unique) solution of the subproblem for some 𝑐_𝑘>0, then 𝑋^𝑘 is a stationary point of the original problem. Conversely, if 𝑋^𝑘 is a stationary point of the original problem, then Δ𝑋^𝑘=0 is the unique solution of subproblem for every 𝑐_𝑘>0.

Proof. The statements follow immediately from the preceding arguments.

3. Convergence Analysis

Throughout this section, we assume that Algorithm 1 generates an infinite sequence {𝑋^𝑘}. Our aim is to establish a global convergence result for Algorithm 1. More precisely, we will show any accumulation point of {𝑋^𝑘} is a stationary point of the original problem.

Lemma 6. Let {𝑋^𝑘} be a sequence generated by Algorithm 1, and let 𝜇_𝑘 and {𝑋^𝑘}_𝑘∈𝐾 be subsequences converging to 0 and some matrix 𝑋^∗, respectively, in such a way that {𝑐_𝑘‖Δ𝑋^𝑘‖}_𝑘∈𝐾→0. Then 𝑋^∗ is a stationary point of the original problem.

Proof. First note that 𝑋^∗ is symmetric positive semidefinite and hence feasible for original problem. Furthermore, since 𝑐_𝑘≥𝑐_min>0 for all 𝑘∈𝑁, the assumption {𝑐_𝑘‖Δ𝑋^𝑘‖}_𝑘∈𝐾→0 implies {‖Δ𝑋^𝑘‖}_𝑘∈𝐾→0. By continuity, we also have 𝐷𝐻_𝜇𝑘(𝑋^𝑘)→𝐷𝐻₀(𝑋^∗) as 𝑘→∞, 𝑘∈𝐾. This together with system (18) implies that 𝑈^𝑘=𝑐_𝑘Δ𝑋^𝑘+𝐷𝑄_𝜇𝑘Δ𝑋^𝑘,𝐷𝑄₀(0)=∶𝑈^∗ (21) as 𝑘→∞, 𝑘∈𝐾. Therefore, taking the limit 𝑘→∞ on the subsequence 𝐾 in the KKT condition (18), we obtain 𝐷𝑄₀(0)−𝑈^∗𝑋=0,_∗⪰0,𝑈^∗𝑈⪰0,_∗,𝑋^∗=0. (22)
Hence we conclude that 𝑋^∗ 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 {𝑋^𝑘} be a sequence generated by Algorithm 1 and {𝑋^𝑘}_𝑘∈𝐾 a subsequence converging to some matrix 𝑋^∗. If 𝑋^∗ is not a stationary point, then one has limsup_{𝑘→∞,𝑘∈𝐾}𝑐_𝑘<∞.

Proof. Let 𝐾∶={𝑘−1∣𝑘∈𝐾}. Then we have {𝑋^𝑘+1}_𝑘∈𝐾→𝑋^∗. We will show that limsup_{𝑘→∞,𝑘∈𝐾}𝑐_𝑘+1<∞. Assume the contrary. Then, if necessary, we may suppose without loss of generality that lim_{𝑘→∞,𝑘∈𝐾}𝑐_𝑘+1=∞. (23) The updating rule in step (S3) then implies that none of the iterations 𝑘∈𝐾 with 𝑘 sufficiently large is successful since otherwise we would have 𝑐_𝑘+1≤𝑐_max for all these 𝑘∈𝐾. Hence, we have 𝑟_𝑘<𝜌₁ (24) and 𝑋^𝑘=𝑋^𝑘+1 for all 𝑘∈𝐾 large enough. Since {𝑋^𝑘+1}_𝑘∈𝐾→𝑋^∗, this implies {𝑋^𝑘}_𝑘∈𝐾→𝑋^∗, too. Further, noticing that 𝑐_𝑘+1=𝜎₂𝑐_𝑘 for all unsuccessful iterations, we also have lim_{𝑘→∞,𝑘∈𝐾}𝑐_𝑘=∞ (25) because of (23). We now want to show that 𝑟_𝑘⟶1as𝑘⟶∞,𝑘∈𝐾, (26) which would then lead to the desired contradiction to (24). To this end, we first note that liminf_{𝑘→∞,𝑘∈𝐾}𝑐_𝑘Δ𝑋^𝑘>0. (27) In fact, if 𝑐_𝑘‖Δ𝑋^𝑘‖→0 on a subsequence, we would deduce from Lemma 6 that 𝑋^∗ is a stationary point in contradiction to our assumption. Hence there is a constant 𝛾>0 such that 𝑐_𝑘Δ𝑋^𝑘≥𝛾,𝑘∈𝐾. (28) By Lemma 2, this implies 𝐻_𝜇𝑘𝑋^𝑘−𝑄_𝜇𝑘Δ𝑋^𝑘≥12𝑐_𝑘Δ𝑋^𝑘,Δ𝑋^𝑘≥12𝛾Δ𝑋^𝑘 (29) for all 𝑘∈𝐾 sufficiently large.
We further note that {‖Δ𝑋^𝑘‖}_𝑘∈𝐾→0. Otherwise, it would follow from (25) that 𝑐_𝑘‖𝑋^𝑘‖²→∞ 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 Δ𝑋^𝑘∶=0 would give a smaller objective value. Hence we have {‖Δ𝑋^𝑘‖}_𝑘∈𝐾→0.
Taking this into account and using {𝑋^𝑘}_𝑘∈𝐾→𝑋^∗ and the fact that 𝐹 is continuously differentiable, we obtain through standard calculus arguments 𝐻_𝜇𝑘𝑋^𝑘−𝑄_𝜇𝑘Δ𝑋^𝑘=𝑜Δ𝑋^𝑘as𝑘⟶∞,𝑘∈𝐾. (30) Summarizing these observations, we get 𝑟_𝑘=𝐻−1^𝜇𝑘(0)−𝐻_𝜇𝑘𝑋^𝑘+Δ𝑋^𝑘𝐻_𝜇𝑘(0)−𝑄_𝜇𝑘Δ𝑋^𝑘=𝑄−1_𝜇𝑘Δ𝑋^𝑘−𝐻_𝜇𝑘𝑋^𝑘+Δ𝑋^𝑘𝐻_𝜇𝑘(0)−𝑄_𝜇𝑘Δ𝑋^𝑘≤𝑜Δ𝑋^𝑘(1/2)𝛾Δ𝑋^𝑘⟶0 (31) as 𝑘→∞, 𝑘∈𝐾. This contradiction to (24) completes the proof.

As a direct consequence of this lemma, we obtain the following result.

Lemma 8. Let {𝑋^𝑘} be sequence generated by Algorithm 1. Then there are infinitely many successful iterations.

Proof. If not, there would exist an index 𝑘₀∈𝑁 with 𝑟_𝑘<𝜌₁ and 𝑋^𝑘=𝑋^𝑘0 for all 𝑘≥𝑘₀. This implies 𝑐_𝑘→∞ due to the updating rule in (S3). However, since 𝑋^𝑘0 is not a stationary point and {𝑋^𝑘}→𝑋^𝑘0, we get a contradiction to Lemma 6.

We are now in the position to prove the main convergence result for Algorithm 1.

Theorem 9. Let {𝑋^𝑘} be a sequence generated by Algorithm 1. Then, any accumulation point of this sequence is a stationary point of the original problem.

Proof. Let 𝑋^∗ be an accumulation point and {𝑋^𝑘}_𝑘∈𝐾 subsequence converging to 𝑋^∗. Since 𝑋^𝑘=𝑋^𝑘+1 for all unsuccessful iterations 𝑘 and since there are infinitely many successful iterations by Lemma 7, we may assume without loss of generality that all iterations 𝑘∈𝐾 are successful.
Assume that 𝑋^∗ is not a solution. Lemma 6 then implies limsup_{𝑘→∞,𝑘∈𝐾}𝑐_𝑘<∞. (32) Hence there is a constant 𝛾>0 such that 𝑐_𝑘≥𝛾,𝑘∈𝐾. (33) Since each iteration 𝑘∈𝐾 is successful, we also have 𝑟_𝑙≥𝜌₁. Consequently, we obtain from Lemma 2𝐻_𝜇𝑘𝑋^𝑘−𝐻_𝜇𝑘𝑋^𝑘+Δ𝑋^𝑘≥𝜌₁𝐻_𝜇𝑘𝑋^𝑘−𝑄_𝜇𝑘Δ𝑋^𝑘≥12𝜌₁𝑐_𝑘Δ𝑋^𝑘,Δ𝑋^𝑘≥12𝜌₁𝑐_minΔ𝑋^𝑘2 (34) for all 𝑘∈𝐾. Since {𝐻_𝜇𝑘(𝑋^𝑘)} is monotonically nonincreasing and bounded from below by, for example, 𝐻₀(𝑋^∗), we have 𝐻_𝜇𝑘(𝑋^𝑘)−𝐻_𝑚𝑘(𝑋^𝑘+1)→0 as 𝑘→∞. Therefore, we obtain {Δ𝑋^𝑘}_𝑘∈𝐾→0 from (34). By (33), this also implies {𝑐_𝑘‖𝑋^𝑘‖}_𝑘∈𝐾→0. But then Lemma 6 shows that 𝑋^∗ 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 ℵ∶={0}∪𝑘∣Φ𝑋^𝑘+1≤max𝜂𝛽_𝑘,𝜇⁻¹Φ𝑋^𝑘+1−Φ_𝜀𝑘𝑋^𝑘+1=𝑘₀=0<𝑘₁<𝑘₂<⋯. (35)

Lemma 10. If 𝐹 is a 𝑃₀ function, then the sequence 𝑋^𝑘 generated by Algorithm 1 remains in the level set 𝐿₀∶=𝑋∈𝑆∣𝐻(𝑋)≤(1+𝜀)²𝐻𝑋⁰. (36)

Proof. Let 𝑘 be an arbitrary nonnegative integer, and let 𝑘_𝑗 be the largest number in ℵ such that 𝑘_𝑗≤𝑘, as ℵ is as defined from (35). It is easy to deduce from step 3 of Algorithm 1 that 𝜇_𝑘=𝜇_𝑘𝑗,𝛽_𝑘=𝛽_𝑘𝑗,as𝑘_𝑗≤𝑘<𝑘_𝑗+1,Φ_𝜇𝑘𝑋^𝑘≤Φ_𝜇𝑘𝑗𝑋^𝑘,as𝑘_𝑗≤𝑘<𝑘_𝑗+1. (37) Set 𝑈_𝑗∶=𝑋∈𝑆∣Φ_𝜇𝑘𝑋^𝑘≤Φ_𝜇𝑘𝑗𝑋^𝑘. (38) As 𝑘 is an arbitrary integer and 𝑥_𝑘∈𝑈_𝑗, it follows that 𝑈_𝑗⊆𝐿₀.
Next, by induction, we will prove 𝑈_𝑗⊆𝐿₀,∀𝑗≥0. (39)
In view of Lemma 2, we deduce that ∀𝑋∈𝑈_𝑗, ‖Φ(𝑋)‖≤Φ_𝜇𝑘𝑗𝑋^𝑘+𝜅𝜇_𝑘𝑗≤Φ_𝜇𝑘𝑗𝑋^𝑘𝑗+𝜅𝜇_𝑘𝑗≤Φ𝑋^𝑘𝑗+2𝜅𝜇_𝑘𝑗≤𝛽_𝑘𝑗+𝜅𝜇_𝑘𝑗. (40) If 𝑗=0, then by (40) we have ‖Φ(𝑋)‖≤Φ𝑋⁰+2𝜅√𝜇₀≤(1+𝜀)Φ𝑋⁰,∀𝑋∈𝑈₀. (41) This proves 𝑈₀∈𝐿₀.
Suppose 𝑈_𝑗−1⊆𝐿₀ for some 𝑗>0. Then, 𝑥_𝑘𝑗−1∈𝐿₀ and hence 𝛽_𝑘𝑗−1≤Γ₀. Set ℵ₁∶=𝑘∈ℵ∣𝜂𝛽_𝑘−1≥𝜇⁻¹Φ𝑋^𝑘−Φ_𝜇𝑘−1𝑋^𝑘,ℵ₂∶=𝑘∈ℵ∣𝜂𝛽_𝑘−1<𝜇⁻¹Φ𝑋^𝑘−Φ_𝜇𝑘−1𝑋^𝑘. (42) It follows from step 3 of Algorithm 1 and Lemma 2 that 𝛽_𝑘𝑗≤𝜂𝛽_𝑘𝑗−1=𝜂𝛽_𝑘𝑗−1,if𝑘_𝑗∈ℵ₁, (43) or 𝛽_𝑘𝑗≤𝜅𝜇𝜇_𝑘𝑗−1=𝜅𝜇𝜇_𝑘𝑗−1≤12𝐶₀𝛽²_𝑘𝑗−1≤12𝛽²_𝑘𝑗−1,if𝑘_𝑗∈ℵ₂. (44) This implies that 𝛽_𝑘𝑗≤𝛿₃𝛽_𝑘𝑗−1,∀𝑘_𝑗∈ℵ, (45) where 𝛿₃=max{1/2,𝜂}. Moreover, we have 𝜇_𝑘𝑗≤14𝜇_𝑘𝑗−1=14𝜇_𝑘𝑗−1,∀𝑘_𝑗∈ℵ. (46)
From (45) and (46), we deduce that, for 𝑗>0, 𝛽_𝑘𝑗≤𝛿^𝑗₃𝛽₀=𝛿^𝑗₃Φ𝑋⁰,𝜇_𝑘𝑗≤14^𝑗𝜇₀≤𝜀²4^𝑗+1Γ₀𝜅²Φ𝑋⁰⁴≤𝜀²4^𝑗+1𝜅²Φ𝑋⁰². (47) Combining (47) with (40), we have ‖Φ(𝑋)‖≤𝛿^𝑗₃Φ𝑋⁰+𝜀2^𝑗Φ𝑋⁰≤𝛿^𝑗₃(1+𝜀)Φ𝑋⁰,∀𝑋∈𝑈_𝑗, (48) which implies 𝑈_𝑗⊆𝐿₀. Hence (39) is proved and Lemma 10 is valid.

It follows from Fischer [13] that if 𝐹 is a 𝑃₀ function or, more generally, an 𝑅₀-function, then the level set 𝐿₀ as defined in Lemma 10 is compact. Lemma 3 shows that the sequence {𝐻_𝜇𝑘} is monotonically decreasing and converges.

Theorem 11. Assume that 𝐹 is a 𝑃₀ function. Let {𝑋^𝑘} be a sequence generated by Algorithm 1. If there exists at least an accumulation point in the sequence {𝑋^𝑘}, then the index set ℵ defined by (35) is infinite: lim_𝑘→∞𝜇_𝑘=0,lim_𝑘→∞Φ𝑋^𝑘=0,lim_𝑘→∞Φ_𝜇𝑘𝑋^𝑘=0. (49)

Proof. We first prove that set ℵ is infinite. By contradiction, assume that ℵ is finite. Let 𝑘 be the largest number in ℵ. Then for all 𝑘≥𝑘, 𝜇_𝑘=𝜇_𝑘 and 𝛽_𝑘=𝛽_𝑘. Denote 𝜇∶=𝜇_𝑘,𝛽∶=𝛽_𝑘,𝜙(𝑋)∶=Φ(𝑋)−Φ_𝜇(𝑋). (50) Then ∀𝑘>𝑘, Φ𝑋^𝑘>max𝜂𝛽,𝜇⁻¹𝜙(𝑋),Φ𝑋^𝑘=𝜙𝑋^𝑘+Φ_𝜀𝑋^𝑘. (51)(52)
From Theorem 9, it follows that there exists at least an accumulation point 𝑋∈𝐿₀ of {𝑋^𝑘} such that ∇𝐻_𝜇𝑋=0. (53)
Next, assume that subsequence {𝑋_𝑘}_ℵ𝑘∈ converges to 𝑋. In view of (53), we have {Φ_𝜇(𝑋^𝑘)}_ℵ𝑘∈→Φ_𝜇(𝑋)=0 and hence there exists 𝑘≥𝑘 such that, for all 𝑘∈ℵ with 𝑘≥𝑘, Φ_𝜇𝑋^𝑘≤1−𝜇𝜂𝛽. (54) This together with (51) and (52) shows that, for all 𝑘∈ℵ with 𝑘≥𝑘, Φ_𝜇𝑋^𝑘≤(1−𝜀)Φ𝑋^𝑘≤(1−𝜀)Φ_𝜇𝑋^𝑘+𝜙𝑋^𝑘; (55) that is, Φ_𝜀𝑋^𝑘<𝜀⁻¹−1𝜙𝑋^𝑘, (56) which means Φ𝑋^𝑘≤Φ_𝜀𝑋^𝑘+𝜙𝑋^𝑘<𝜀⁻¹𝜙𝑋^𝑘,as𝑘∈ℵ,𝑘≥𝑘. (57) This contradicts (51). Hence the set ℵ is infinite.
Next, {𝜇_𝑘}→0 follows immediately from the updating rule of 𝜇_𝑘 and the fact that the set ℵ is infinite. Moreover, by the proof of Lemma 10, we deduce Φ𝑋^𝑘≤𝛿^𝑗₄(1+𝜀)Φ𝑋⁰,as𝑘_𝑗≤𝑘<𝑘_𝑗+1. (58) Because the set ℵ is infinite, it follows from Lemma 10 and (58) that lim_𝑘→∞Φ𝑋^𝑘=0,lim_𝑘→∞Φ_𝜇𝑘𝑋^𝑘=0. (59) This completes the proof.

As a consequence of the above theorem, we get the following global convergence result.

Corollary 12. Assume that 𝐹 is a 𝑃₀ function. Let {𝑋^𝑘} be a sequence generated by Algorithm 1. Then every accumulation point of the sequence {𝑋^𝑘} is a solution of NCP(𝐹).

4. Numerical Experiments

4.1. The Reformulation for Subproblem

To test the numerical performance of Algorithm 1, 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 min_{Δ𝑋∈𝑆𝑛×𝑛}12𝑐_𝑘⟨Δ𝑋,Δ𝑋⟩+𝑄_𝜇𝑘(Δ𝑋)s.t.𝑋^𝑘+Δ𝑋≥0. (60) 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 min^𝑛𝑠_𝑗=1𝐶^𝑠_𝑗,𝑋^𝑠_𝑗+^𝑛𝑞_𝑖=1𝑐^𝑞_𝑖^𝑇𝑥^𝑞_𝑖+𝑐^𝑙𝑥^𝑙s.t.^𝑛𝑠_𝑗=1𝐴^𝑠_𝑗^𝑇svec𝑋^𝑠_𝑗+^𝑛𝑞_𝑖=1𝐴^𝑞_𝑖^𝑇𝑥^𝑞_𝑖+𝐴^𝑙𝑇𝑥^𝑙𝑋=𝑏,_𝑠𝑗∈𝑆^{𝑠𝑗×𝑠𝑗}₊,∀𝑗,𝑥^𝑞_𝑖∈𝐾^𝑞𝑖_𝑞∀𝑖,𝑥^𝑙∈𝑅^𝑛𝑙₊, (61) where 𝑋^𝑠_𝑗, 𝑋^𝑠_𝑗 are symmetric matrices of dimension 𝑠_𝑗; 𝑐^𝑞_𝑖, 𝑥^𝑞_𝑖 are vectors in 𝑅^𝑞𝑖; 𝑆^{𝑠𝑗×𝑠𝑗}₊ denotes the 𝑠_𝑗-dimension positive semidefinite cone defined by 𝑆^{𝑠𝑗×𝑠𝑗}₊∶={𝑋∈𝑆^{𝑠𝑗×𝑠𝑗}∶𝑋⪰0}; 𝐷^𝑞𝑖_𝑞 denotes the 𝑞_𝑖-dimensional second-order cone defined by 𝐾^𝑞𝑖_𝑞∶={𝑥=(𝑥₁,𝑥^𝑇_2∶𝑞𝑖)^𝑇∈𝑅^𝑞𝑖∶𝑥₁≥‖𝑥_2∶𝑞𝑖‖}; 𝑐^𝑙 and 𝑥^𝑙 are vectors in 𝑅^𝑛𝑙; 𝐴^𝑠_𝑗 are 𝑠_𝑗×𝑚 matrices with 𝑠_𝑗=𝑠_𝑗(𝑠_𝑗+1)/2; 𝐴^𝑞_𝑖 and 𝐴^𝑙 are 𝑞_𝑖×𝑚 and 𝑙×𝑚 matrices, respectively; and svec is the operator defined by svec(𝑋)∶=(𝑋(1,1),√2𝑋(1,2),𝑋(2,2),√2𝑋(1,3),√2𝑋(2,3),𝑋(3,3),…)^𝑇∈𝑅^{𝑛(𝑛+1)/2} for any symmetric matrix 𝑋∈𝑆^𝑛×𝑛.

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 variable 𝑆∈𝑆^𝑛×𝑛 and set 𝑋^𝑘+Δ𝑋=𝑆. Because Δ𝑋 needs only to be symmetric and not to be positive semidefinite, we set Δ𝑥=svec(Δ𝑋) and write the problem in terms of Δ𝑥∈𝑅^𝑛 with 𝑛∶=𝑛(𝑛+1)/2. Then problem (60) is equivalent to 1min2𝑐_𝑘+svec𝐷Φ_𝜇𝑘𝑋^𝑘2‖Δ𝑥‖²+2svecΦ_𝜇𝑘𝑋^𝑘𝐷Φ_𝜇𝑘𝑋^𝑘𝑇𝑋Δ𝑥s.t.svec^𝑘+Δ𝑥=svec(𝑆),Δ𝑥∈𝑅^𝑛,𝑠⪰0. (62) By introducing the second-order cone constraint ‖Δ𝑥‖𝑡, the above problem can be further rewritten as 1min2𝑐_𝑘+svec𝐷Φ_𝑚𝑢𝑘𝑋^𝑘2𝑡²+2svecΦ_𝜇𝑘𝑋^𝑘𝐷Φ_𝜇𝑘𝑋^𝑘𝑇𝑋Δ𝑥s.t.svec^𝑘(𝑆),‖‖+Δ𝑥=svecΔ𝑥≤𝑡,Δ𝑥∈𝑅^𝑛,𝑆⪰0,𝑡∈𝑅. (63)

Unfortunately, the term 𝑡² is not linear as required in (61). So we replace 𝑡² by the new variable 𝑠≥0 and add the constraint 𝑡²≤𝑠. But this constraint can be rewritten as the semidefinite constraint 𝑠𝑡𝑡1⪰0. (64) Introducing once again an auxiliary variable, problem (62) and hence the original subproblem (60) are equivalent to 1min2𝑐_𝑘+svec𝐷Φ_𝜇𝑘𝑋^𝑘2𝑠+2svecΦ_𝜇𝑘𝑋^𝑘𝐷Φ_𝜇𝑘𝑋^𝑘𝑇𝑋Δ𝑥s.t.−svec^𝑘=Δ𝑥−svec(𝑆),‖Δ𝑥‖≤𝑡,𝑠𝑡𝑡1−𝑊=0,Δ𝑥∈𝑅^𝑛,𝑊⪰0,𝑆⪰0,𝑡∈𝑅,𝑠∈𝑅₊. (65) We write the equality constraint 𝑠𝑡𝑡1−𝑊=0, (66) in the svec-notation. Then we get 1min2𝑐_𝑘+svec𝐷Φ_𝜇𝑘𝑋^𝑘2𝑡+2svecΦ_𝜇𝑘𝑋^𝑘𝐷Φ_𝜇𝑘𝑋^𝑘𝑇𝑋Δ𝑥s.t.−svec^𝑘(𝑆),‖‖𝑠√=Δ𝑥−svecΔ𝑥≤𝑡,𝑡1−svec(𝑊)=0,Δ𝑥∈𝑅^𝑛,𝑊⪰0,𝑆⪰0,𝑡∈𝑅,𝑠∈𝑅₊. (67)

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 𝑛_𝑠∶=2,𝑛_𝑞∶=1,𝑠₁𝑠∶=𝑛,₂∶=2,𝑞₁∶=1+𝑛,𝑙∶=1+2𝑚. (68) The variables are given by 𝑋^𝑠₁∶=𝑆∈𝑆₊^𝑛×𝑛,𝑋^𝑠₂∶=𝑊∈𝑆^2×2,𝑥^𝑞₁∶=𝑡,Δ𝑥^𝑇𝑇∈𝐾^1+𝑛_𝑞,𝑥^𝑙∶=𝑠,𝜉^𝑇,𝜔^𝑇𝑇∈𝑅^1+2𝑚. (69) The input datum in the objective function is given by 𝐶^𝑠₁∶=0∈𝑆^𝑛×𝑛,𝐶^𝑠₂∶=0∈𝑆^2×2,𝑐^𝑞₁∶=0,svec𝐷𝑓𝑋^𝑘𝑇𝑇∈𝑅^1+𝑛,𝑐^𝑙1∶=2𝑐_𝑘,𝛼_𝑘𝑒,0∈𝑅^1+2𝑚 (70) with 𝑒=(1,…,1)^𝑇∈𝑅^𝑚. (71) Finally, the matrices 𝐴^𝑠₁∈𝑅^{𝑛×(𝑛+3+𝑚)}, 𝐴^𝑠₂∈𝑅^{3×(𝑛+3+𝑚)}, 𝐴^𝑞₁∈𝑅^{(1+𝑛)×(𝑛+3+𝑚)}, and 𝐴^𝑙∈𝑅^{(1+𝑚+𝑚)×(𝑛+3+𝑚)} and the vector 𝑏∈𝑅^𝑛+3+𝑚 are given by 𝐴^𝑠₁=−𝐼|0|0,𝐴^𝑠₂=0||0,𝐴−𝐼_𝑞1=0𝐼0√(𝑊)(𝑤),𝐴20000⋮⋮⋮0000⋯0−svec⋯−svec_𝑙=00000𝐼,𝑋1000⋯0−𝐼𝑏=−svec^𝑘𝑇00−1𝑔𝑋^𝑘𝑇𝑇. (72) This is the desired reformulation.

It may be worth mentioning that problem (60) can also be transformed as 1min2𝑐_𝑘+svec𝐷Φ_𝜇𝑘𝑋^𝑘2𝑡+2svecΦ_𝜇𝑘𝑋^𝑘𝐷Φ_𝜇𝑘𝑋^𝑘𝑇𝑋Δ𝑥s.t.−svec^𝑘=Δ𝑥−svec(𝑆),‖Δ𝑥‖≤𝑡,𝑠𝑡𝑡1−𝑊=0,Δ𝑥∈𝑅^𝑛,𝑊⪰0,𝑆⪰0,𝑡∈𝑅,𝑠∈𝑅₊. (73) Since the constraint ‖Δ𝑥‖²≤𝑡 is equivalent to 𝑡Δ𝑥^𝑇,Δ𝑥𝐼 (74) 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 Algorithm 1. 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 with 𝐹(𝑋)∶=𝑋𝑋−2𝑋+𝐼. It is obvious that the solution set of SDCP(𝐹) is nonempty, since 𝑋=𝐼 is its one solution. The parameters in the algorithm can be presented as follows: 𝜌₁ and 𝜌₂ can be randomly generated from [0.1,0.5] and [0.6,1], respectively; 𝜎₁ and 𝜎₂ are randomly generated from [0.5,1] and [1.5,2], respectively; 𝑐_min and 𝑐_max are randomly generated from [0.01,1] and [500,1000], respectively; 𝜇 is randomly chosen from [10,20]; 𝑐₀ is randomly generated from [𝑐_min,𝑐_max] and 𝑋⁰∶=𝐴^𝑇𝐴, where 𝐴∈𝑅^𝑚×𝑛 with every entry being randomly generated from [0,1]. The stopping criterion is set as ‖𝑋^𝑘‖≤10⁻⁶ and 𝜇_𝑘≤10⁻⁶.

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) Choose 0<𝜌<1, 0<𝛼<1, 𝑋⁰∈𝑆₊, Γ₀∶=(1+𝜇)‖Φ₀(𝑋⁰)‖, 𝛽₀∶=‖Φ₀(𝑋⁰)‖, 𝜅∶=√2𝑛, 𝜇₀∶=((𝜀/(2Γ₀𝜅))𝛽²₀)², and set 𝑘∶=0. (S1)Find the solution Δ𝑋^𝑘∈𝑆 of the subproblem min_{Δ𝑋∈𝑆𝑛×𝑛}12⟨Δ𝑋,Δ𝑋⟩+𝑄_𝜇𝑘(Δ𝑋)s.t.𝑋^𝑘+Δ𝑋⪰0. (75) If Δ𝑋^𝑘=0, then STOP. (S2)Compute 𝛼_𝑘=max{1,𝛼,𝛼²,…} such that 𝐻_𝜇𝑘𝑋^𝑘≥𝐻_𝜇𝑘𝑋^𝑘+𝛼_𝑘Δ𝑋^𝑘+𝜌𝛼_𝑘𝑄_𝜇𝑘Δ𝑋^𝑘. (76) (S3)Let 𝑥^𝑘+1∶=𝑋^𝑘+𝛼_𝑘Δ𝑋^𝑘. If Φ₀𝑋^𝑘+1≤max𝜂𝛽_𝑘,𝜀⁻¹Φ₀𝑋^𝑘+1−Φ_𝜇𝑘𝑋^𝑘+1, (77) then set 𝛽_𝑘+1∶=Φ₀𝑋^𝑘+1 (78) and choose 𝜇_𝑘+1 such that 0<𝜇_𝑘+1𝜇≤min2Γ₀𝜅𝛽²_𝑘+1,𝜇_𝑘4; (79) otherwise, let 𝛽_𝑘+1∶=𝛽_𝑘 and 𝜇_𝑘+1∶=𝜇_𝑘. (S4)Set 𝑘←𝑘+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 𝑋⁰ and the stopping criteria are the same as Algorithm 1.

We now solve this problem 40 times by Algorithms 1 and 13, respectively, with the initial point 𝑋⁰ 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.

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.

Conflict of Interests

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

Acknowledgments

This work was supported by the Natural Scientific Foundation of China (no. 71201040 and no. 11201099) and by the Major Program of the National Natural Science Foundation of China (no. 71031003).