A descent method for a reformulation of the second-order cone complementarity problem

Journal of Computational and Applied Mathematics, vol. 213, pp. 547-558, 2008.

A descent method for a reformulation of the second-order cone complementarity problem

Jein-Shan Chen ¹ Department of Mathematics National Taiwan Normal University

Taipei 11677, Taiwan

Shaohua Pan²

School of Mathematical Sciences South China University of Technology

Guangzhou 510640, China

July 28, 2006

(first revised September 19, 2006) (second revised January 31, 2007)

Abstract Analogous to the nonlinear complementarity problem (NCP) and the semi- definite complementarity problem (SDCP), a popular approach to solving the second- order cone complementarity problem (SOCCP) is to reformulate it as an unconstrained minimization of a certain merit function over IRⁿ. In this paper, we present a descent method for solving the unconstrained minimization reformulation of the SOCCP which is based on the Fischer-Burmeister merit function associated with second-order cone [4], and prove its global convergence. Particularly, we compare the numerical performance of the method for the symmetric affine SOCCP generated randomly with the Fischer- Burmeister merit function approach [4]. The comparison results indicate that, if a scaling strategy is imposed on the test problem, the descent method proposed is comparable with the merit function approach in the CPU time for solving test problems although the former may require more function evaluations.

Key words. Second-order cone, complementarity, merit function, descent method

1Member of Mathematics Division, National Center for Theoretical Sciences, Taipei Office.

The author’s work is partially supported by National Science Council of Taiwan. E-mail:

jschen@math.ntnu.edu.tw, FAX: 886-2-29332342.

2E-mail: shhpan@scut.edu.cn.

1 Introduction

In this paper, we consider the following SOCCP of finding ζ ∈ IRⁿ satisfying

hF (ζ), ζi = 0, F (ζ) ∈ K, ζ ∈ K, (1)

where h·, ·i is the Euclidean inner product, F : IRⁿ→ IRⁿ is a smooth (i.e., continuously differentiable) mapping, and K is the Cartesian product of second-order cones (SOC), also called Lorentz cones [9]. In other words,

K = Kⁿ¹ × · · · × K^nm, (2)

where m, n₁, . . . , n_m ≥ 1, n₁+ · · · + n_m = n, and

Kⁿⁱ := {(x₁, x₂) ∈ IR × IRⁿⁱ⁻¹ | kx₂k ≤ x₁}, (3) with k · k denoting the Euclidean norm and K¹ denoting the set of nonnegative reals IR+. A special case of (2) is K = IRⁿ₊, the nonnegative orthant in IRⁿ, which corresponds to m = n and n₁ = · · · = n_m = 1. If K = IRⁿ₊, then (1) reduces to the nonlinear complementarity problem (NCP). The NCP plays a fundamental role in optimization theory and has many applications in engineering and economics; see, e.g., [7, 10, 11, 12].

Unless otherwise stated, in the first three sections of this paper, we assume K = Kⁿ for simplicity, i.e., K is a single second-order cone (all the analysis can be carried over to the case where K is a product of second-order cones without difficulty).

There have been proposed various methods for solving the SOCCP. They include interior-point methods [1, 25, 28, 29, 31], reformulating SOC constraints as smooth con- vex constraints [32], and (non-interior) smoothing Newton methods [6, 16, 20]. These methods require solving a nontrivial system of linear equations at each iteration. In the recent paper [4], an alternative approach based on reformulating the SOCCP as an un- constrained smooth minimization problem was studied. In particular, they were finding a smooth function ψ : IRⁿ× IRⁿ → IR₊ such that

ψ(x, y) = 0 ⇐⇒ hx, yi = 0, x ∈ Kⁿ, y ∈ Kⁿ. (4) We call such a ψ a merit function. Then SOCCP can be expressed as an unconstrained smooth (global) minimization problem:

ζ∈IRminⁿ ψ(F (ζ), ζ). (5)

Various gradient methods such as conjugate gradient methods and quasi-Newton methods [2, 15] can be applied to (5). For this approach to be effective, the choice of ψ is crucial.

In the case of NCP, a popular choice is ψ_FB(a, b) = 1

2

Xn i=1

φ_FB(ai, bi)²

for all a = (a₁, ..., a_n)^T ∈ IRⁿ and b = (b₁, ..., b_n)^T ∈ IRⁿ, where φ_FB is the well-known Fischer-Burmeister NCP-function [13, 14] defined by

φ_FB(a_i, b_i) =

q

a²_i + b²_i − a_i− b_i.

It has been shown that ψ_FB is smooth (even though φ_FB is not differentiable) and is a merit function for NCP [8, 22, 23]. These two functions can be extended to the case of SOCCP via Jordan algebra shown as below. For any x = (x₁, x₂), y = (y₁, y₂) ∈ IR × IRⁿ⁻¹, we define their Jordan product associated with Kⁿ as

x ◦ y := (hx, yi, y₁x₂+ x₁y₂).

The identity element under this product is e := (1, 0, . . . , 0)^T ∈ IRⁿ. We write x² to mean x ◦ x and write x + y to mean the usual componentwise addition of vectors. It is known that x² ∈ Kⁿ for all x ∈ IRⁿ. Moreover, if x ∈ Kⁿ, then there exists a unique vector in Kⁿ, denoted by x^1/2, such that (x^1/2)² = x^1/2◦ x^1/2 = x. Then,

φ_FB(x, y) := (x² + y²)^1/2− x − y (6) is well-defined for all (x, y) ∈ IRⁿ× IRⁿ and maps IRⁿ× IRⁿ to IRⁿ. It was shown in [16]

that φ_FB(x, y) = 0 if and only if hx, yi = 0, x ∈ Kⁿ, y ∈ Kⁿ. Hence, ψ_FB : IRⁿ× IRⁿ→ IR₊ given by

ψ_FB(x, y) := 1

2kφ_FB(x, y)k², (7)

is a merit function for SOCCP because ψ_FB satisfies (4) as well. Therefore, the SOCCP is equivalent to the global minimization problem:

ζ∈IRminⁿ f_FB(ζ) := ψ_FB(F (ζ), ζ). (8)

It was also shown in the paper [4] that, like the NCP case, ψ_FB is smooth and, when

∇F is positive semi-definite, every stationary point of (8) solves the SOCCP. For SDCP, which is a natural extension of NCP where IRⁿ₊ is replaced by the cone of positive semi- definite matrices S₊ⁿ and the partial order ≤ is also changed by ¹_Sn

+ (a partial order associated with S₊ⁿ where A ¹_Sn

+ B means B − A ∈ S₊ⁿ) accordingly, the aforementioned features hold for the following analog of the SDCP merit function studied by Yamashita and Fukushima in [33]:

ψ_YF(x, y) := ψ0(hx, yi) + ψ_FB(x, y), (9) where ψ₀ : IR → [0, ∞) is any smooth function satisfying

ψ0(t) = 0 ∀t ≤ 0 and ψ₀⁰(t) > 0 ∀t > 0. (10) In [33], ψ₀(t) = ¹₄(max{0, t})⁴ was considered. In fact, the function ψ_YF, which was recently studied in [4], is also a SOCCP version merit function that enjoys favorable

properties as what ψ_FB has. Moreover, ψ_YF possesses properties of bounded level sets and error bound.

In this paper, we focus on the following equivalent reformulation of SOCCP, which arises via the merit function ψ_YF defined as in (9)-(10):

ζ∈IRminⁿf_YF(ζ) := ψ_YF(F (ζ), ζ). (11) We are motivated by the work [33] showing a descent method for the SDCP. Thus, the main purpose of the paper is to explore the extension to SOCCP. In other words, we wish to adopt the algorithm therein to solve the equivalent reformulation (11) of the SOCCP and prove its global convergence (see Sec. 3). In particular, we also compare the numer- ical performance of the descent algorithm for the symmetric affine SOCCPs generated randomly with the Fischer-Burmeister merit function approach [4]. Here it is worth of pointing out that the proposed algorithm does not work for the other reformulation (8).

The reason is that f_FB(ζ) lacks property of bounded level sets and does not provide error bound due to the absence of the term ψ0.

Some words about our notation. Throughout this paper, IRⁿ denotes the space of n-dimensional real column vectors. For any differentiable function f : IRⁿ → IR, ∇f (x) denotes the gradient of f at x. For any differentiable mapping F : IRⁿ→ IR^m, ∇F (x) is a n × m matrix which denotes the transposed Jacobian of F at x.

2 Preliminaries

As mentioned in the introduction, ψ_YF satisfies (4), so the SOCCP can be recast as an equivalent global minimization (11). It was shown in [4] that the function f_YF is smooth, has bounded level sets, and provides error bound for the unconstrained minimization reformulation. Moreover, every stationary point of problem (11) is a solution of the SOCCP. In this section, we review some basic concepts and properties that will be used for proving the convergence results of the descent algorithm later. Since the work of [4]

already includes as special cases the following lemmas, we here omit the proofs.

Lemma 2.1 [4, Prop. 3.2] Let φ_FB, ψ_FB be given by (6) and (7), respectively, and ψ_YF be given by (9)-(10). Then ψ_FB and ψ_YF are both smooth on IRⁿ× IRⁿ.

Lemma 2.2 [4, Prop. 4.2] Let ψ_YF be given by (9)-(10) and f_YF(ζ) be defined as (11).

Then, for every ζ ∈ IRⁿ such that ∇F (ζ) is positive semi-definite, either f_YF(ζ) = 0 or

∇f_YF(ζ) 6= 0 with hd(ζ), ∇f_YF(ζ)i < 0, where d(ζ) := −

µ

ψ⁰₀(hF (ζ), ζi)ζ + ∇_xψ_FB(F (ζ), ζ)

¶

. (12)

In what follows, we say that F is monotone if

hF (ζ) − F (ξ), ζ − ξi ≥ 0 ∀ζ, ξ ∈ IRⁿ and F is strongly monotone if there exists ρ > 0 such that

hF (ζ) − F (ξ), ζ − ξi ≥ ρkζ − ξk² ∀ζ, ξ ∈ IRⁿ.

It is well known that, when F is differentiable, F is monotone if and only if ∇F (ζ) is positive semi-definite for all ζ ∈ IRⁿ while F is strongly monotone if and only if ∇F (ζ) is positive definite for all ζ ∈ IRⁿ.

Lemma 2.3 [4, Prop. 5.2] Suppose that F is a differentiable and monotone mapping from IRⁿ to IRⁿ. Suppose also that the SOCCP (1) is strictly feasible, i.e., there exists ζ ∈ IRb ⁿ such that F (ζ),^b ζ ∈ int(K^{b n}). Then the level set

L(γ) := {ζ ∈ IRⁿ | f_YF(ζ) ≤ γ}

is nonempty and bounded for all γ ≥ 0, where f_YF is given by (11).

Remark 2.1 It is known that Lemma 2.3 is also true if the conditions of monotonicity and strict feasibility is replaced by strong monotonicity.

We next recall some basic results about the spectral factorization associated with Kⁿ. For any x = (x1, x2) ∈ IR × IRⁿ⁻¹, it admits a spectral factorization of the form

x = λ₁(x) · u⁽¹⁾_x + λ₂(x) · u⁽²⁾_x , (13) where λ_i(x) and u⁽ⁱ⁾_x for i = 1, 2 are the spectral values and the associated spectral vectors of x given by

λ_i(x) = x₁+ (−1)ⁱkx₂k, u⁽ⁱ⁾_x =









1 2

µ

1, (−1)ⁱ x₂ kx2k

¶

if x₂ 6= 0

1 2

µ

1, (−1)ⁱw2

¶

if x2 = 0

(14)

with w₂ being any vector in IRⁿ⁻¹ satisfying kw₂k = 1. If x₂ 6= 0, the factorization is unique. The set {u⁽¹⁾_x , u⁽²⁾_x } is called a Jordan frame and has the following properties.

Property 2.1 For any x = (x₁, x₂) ∈ IR × IRⁿ⁻¹ with the spectral values λ₁(x), λ₂(x) and spectral vectors u⁽¹⁾_x , u⁽²⁾_x given as in (14), we have

(a) u⁽¹⁾_x and u⁽²⁾_x are orthogonal under Jordan product and have length 1/√

2 , i.e., u⁽¹⁾_x ◦ u⁽²⁾_x = 0 , ku⁽¹⁾_x k = ku⁽²⁾_x k = 1

√2 .

(b) u⁽¹⁾_x and u⁽²⁾_x are idempotent under Jordan product, i.e. u⁽ⁱ⁾_x ◦ u⁽ⁱ⁾_x = u⁽ⁱ⁾_x for i = 1, 2.

The spectral factorization (13)-(14) of x, as well as x² and x^1/2 have various interesting properties; see [16]. For instances, for any x = (x1, x2) ∈ IR × IRⁿ⁻¹, with spectral values λ₁(x), λ₂(x) and spectral vectors u⁽¹⁾_x , u⁽²⁾_x , the following results hold:

(1) x² = λ²₁(x) u⁽¹⁾_x + λ²₂(x) u⁽²⁾_x ∈ Kⁿ.

(2) If x ∈ Kⁿ, then 0 ≤ λ₁(x) ≤ λ₂(x) and x^1/2=^qλ₁(x) u⁽¹⁾_x +^qλ₂(x) u⁽²⁾_x .

To close this section, we present a property of ψ_FB associated with the spectral value.

Lemma 2.4 [4, Lemma 9(a)] For any {(x^k, y^k)}^∞_k=1 ⊆ IRⁿ × IRⁿ, let λ1(x^k) ≤ λ2(x^k) and µ₁(y^k) ≤ µ₂(y^k) denote the spectral values of x^k and y^k, respectively. Then, if λ₁(x^k) → −∞ or µ₁(y^k) → −∞, we have ψ_FB(x^k, y^k) → ∞.

3 Main Results

In this section, we propose a descent method for solving the unconstrained minimiza- tion reformulation (11) of the SOCCP and prove its global convergence. The proposed method uses d(ζ) defined as (12) as its direction. Now let us describe the algorithm.

Algorithm 3.1:

(Step 0) Choose ζ⁰ ∈ IRⁿ, ε ≥ 0, σ ∈ (0, 1/2), β ∈ (0, 1) and set k := 0.

(Step 1) If f_YF(ζ^k) ≤ ε, then stop.

(Step 2) Compute d(ζ^k) := −

µ

ψ₀⁰(hF (ζ^k), ζ^ki)ζ^k+ ∇_xψ_FB(F (ζ^k), ζ^k)

¶

.

(Step 3) Find a step-size tk := β^mk, where mk is the smallest nonnegative integer m satisfying the Armijo’s rule:

f_YF(ζ^k+ β^md(ζ^k)) ≤ (1 − σβ^2m) f_YF(ζ^k). (15) (Step 4) Set ζ^k+1 := ζ^k+ t_k d(ζ^k), k := k + 1 and go to Step 1.

Note that the above algorithm is ∇F -free, i.e., there is no need to compute the Ja- cobian matrix of F , and moreover, the computation work in each iteration is very small, i.e., only several vector multiplications. In fact, this type of algorithm was also studied for the NCP (see [17]) and the SDCP (see [33]) and the most remarkable feature of this

type of algorithm is that not only the step-size but also the search direction itself is adjusted via the Armijo’s rule. In practical experience, σ is usually chosen close to zero, and β is usually chosen in (₁₀¹ ,¹₂) depending on the confidence we have on the quality of the initial step-size (see [2]).

Next, we prove the global convergence of Algorithm 3.1. Without any loss of gener- ality, we suppose ε = 0 so that the algorithm generates an infinite sequence {ζ^k}.

Proposition 3.1 Suppose that F is monotone and the SOCCP (1) is strictly feasible.

Then the sequence {ζ^k} generated by Algorithm 3.1 has at least one accumulation point, and any accumulation point is a solution of the SOCCP (1).

Proof. The proof is standard and can be found in [2]. For completeness, we here present its proof by the following three steps.

(i) First, we show that, whenever ζ^k is not a solution, there exists a nonnegative integer mk in Step 3 of Algorithm 3.1. Suppose not, then for any positive integer m, we have

f_YF(ζ^k+ β^m d(ζ^k)) − f_YF(ζ^k) > −σβ^2m f_YF(ζ^k).

Dividing by β^m on both sides and letting m → ∞ yields

h∇f_YF(ζ^k), d(ζ^k)i ≥ 0. (16)

Since F is monotone which is equivalent to ∇F (ζ) is positive semi-definite, the inequality (16) contradicts Lemma 2.2. Hence, we can find an integer m_k in Step 3.

(ii) Secondly, we show that the sequence {ζ^k} generated by the algorithm has at least one accumulation point. By the descent property of Algorithm 3.1, the sequence {f_YF(ζ^k)}_k∈N is decreasing. Hence by Lemma 2.3, we have that {ζ^k} is bounded, and consequently has at least one accumulation point.

(iii) Finally, we prove that any accumulation point of {ζ^k} is a solution of the SOCCP (1). Let ζ^∗ be an arbitrary accumulation point of {ζ^k}_k∈N. In other words, there is a subsequence {ζ^k}_k∈K converging to ζ^∗, where K is a subset of N. We know d(·) is continuous (since ψ0 and ψ_FB are smooth) which implies {d(ζ^k)}k∈K converges to d(ζ^∗).

Next, we need to discuss two cases. First, we consider the case where there exists a constant ¯β such that β^mk ≥ ¯β > 0 for all k ∈ K. Then, from (15), we have

f_YF(ζ^k+1) ≤ (1 − σ ¯β²) f_YF(ζ^k)

for all k ∈ K and the entire sequence {f_YF(ζ^k)}_k∈K is decreasing. Thus, we obtain f_YF(ζ^∗) = 0 (by taking the limit) which says ζ^∗ is a solution of the SOCCP (1). Now, we consider the other case where there exists a further subsequence such that β^mk → 0.

Note that by Armijo’s rule (15) in Step 3, we have

f_YF(ζ^k+ β^mk−1 d(ζ^k)) − f_YF(ζ^k) > −σβ^2(mk−1)f_YF(ζ^k).

Dividing by β^mk−1 both sides and passing the limit on the further subsequence, we obtain h∇f_YF(ζ^∗), d(ζ^∗)i ≥ 0,

which yields that ζ^∗ is a solution of the SOCCP (1) by Lemma 2.2. 2

Proposition 3.2 Let F be a continuously differentiable and strongly monotone mapping.

Then the sequence {ζ^k} generated by Algorithm 3.1 converges to the unique solution of the SOCCP (1).

Proof. The proof is routine (see [7]), however, we present it for completeness. We know that the property of bounded level sets is also held when F is strongly monotone, so following the same arguments as in the proof of Prop. 3.1, we again obtain that {ζ^k} has at least one accumulation point and any accumulation point is a solution of the SOCCP (1).

On the other hand, the strong monotonicity of F implies that the SOCCP (1) has at most one solution. To see this, say there are two solutions ζ^∗, ξ^∗ ∈ IRⁿ such that

( hF (ζ^∗), ζ^∗i = 0,

F (ζ^∗) ∈ Kⁿ, ζ^∗ ∈ Kⁿ and

( hF (ξ^∗), ξ^∗i = 0, F (ξ^∗) ∈ Kⁿ, ξ^∗ ∈ Kⁿ. Since F is strongly monotone, we have hF (ζ^∗) − F (ξ^∗), ζ^∗− ξ^∗i > 0. However,

hF (ζ^∗) − F (ξ^∗), ζ^∗− ξ^∗i

= hF (ζ^∗), ζ^∗i + hF (ξ^∗), ξ^∗i − hF (ζ^∗), ξ^∗i − hF (ξ^∗), ζ^∗i

= −hF (ζ^∗), ξ^∗i − hF (ξ^∗), ζ^∗i

≤ 0,

where the inequality is due to F (ζ^∗), ζ^∗, F (ξ^∗), ξ^∗ are all in K. Hence, it is a contradiction and therefore there is at most one solution for the SOCCP (1).

From all the above, it says there is a unique solution ζ^∗, so the entire sequence {x^k} must converge to ζ^∗. 2

Prop. 3.1-3.2 may not be so surprising since they seems as expected. Nonetheless, we do not take them for granted before we prove them even though we think they should be true. Now, the results of Prop. 3.1-3.2 do fill up the gap in the literature. We notice that Lemma 2.3 plays an important role in the proofs for them. In fact, the assumption of strict feasibility is necessary for Lemma 2.3 to be held. For example, when F (ζ) ≡ 0, every ζ ∈ Kⁿ is a solution of SOCCP (1) and hence the solution set is unbounded. In the following, we continue a further study of considering another (weaker) condition to replace this kind of “strict” condition by F being a R₀₁-function (will be defined in Def.

3.1) that is a new concept recently developed for linear and nonlinear transformations on Euclidean Algebra [19, 26, 30].

Definition 3.1 For a mapping F : IRⁿ → IRⁿ, it is called a (a) R01-function if for any sequence {ζ^k} such that

kζ^kk → ∞, (−ζ^k)₊

kζ^kk → 0, (−F (ζ^k))₊

kζ^kk → 0, (17)

we have

lim inf

k→∞

hζ^k, F (ζ^k)i kζ^kk² > 0;

(b) R₀₂-function if for any sequence {ζ^k} such that (17), we have lim inf

k→∞

ω(ζ^k◦ F (ζ^k)) kζ^kk² > 0.

The above concepts are extensions of the ones defined for NCP and for SDCP. It is also known that every R₀₁-function is R₀₂-function [24, Lemma 4]; and if F has the uniform Jordan P -property (see [19, 26, 30]), then F is R₀₂-function [24, Lemma 5]. However, it is not clear whether uniform P -property (see [19, 26, 30]) implies R₀₂-function or not.

With this new concept, Lemma 2.3 and Prop. 3.1 can be improved as Lemma 3.1 and Prop. 3.3, respectively. These results are significant not only they are brand-new but also there is no needs the assumption of strict feasibility therein.

Lemma 3.1 Let f_YF be given as in (11). Suppose that F is a R₀₁-function. Then the level set

L(γ) := {ζ ∈ IRⁿ | f_YF(ζ) ≤ γ}

is bounded for all γ ≥ 0.

Proof. We will prove this result by contradiction. Suppose there exists an unbounded sequence {ζ^k} ⊂ L(γ) for some γ ≥ 0. It can be seen that the sequence of the smaller spectral values of {ζ^k} and {F (ζ^k)} are bounded below. In fact, if not, it follows form Lemma 2.4 that f_YF(ζ^k) → ∞, which contradicts {ζ^k} ⊂ L(γ). Therefore, {(−ζ^k)₊} and {(−F (ζ^k))+} are bounded above, which says the conditions in (17) are satisfied. Then, by the assumption of R₀₁-function, we have

lim inf

k→∞

hζ^k, F (ζ^k)i kζ^kk² > 0.

This yields hζ^k, F (ζ^k)i → ∞, and hence f_YF(ζ^k) → ∞ by definition of f_YF given as in (11). Thus, it is a contradiction to {ζ^k} ⊂ L(γ). 2

Proposition 3.3 Let F be a continuously differentiable mapping. Suppose that F is R01-function. Then the sequence {ζ^k} generated by Algorithm 3.1 has at least one accu- mulation point, and any accumulation point is a solution of the SOCCP (1).

Proof. By applying Lemma 3.1 and follow the same arguments as in Prop. 3.1, the desired results hold. We omit it. 2

From [24, 30], the condition of R₀₁-function is weaker than strong monotonicity, and it is also weaker than monotonicity plus strict feasibility in certain sense. However, it is not clear yet whether R₀₁-function can be replaced by R₀₂-function in our brand-new results.

4 Numerical results

In this section, we report our computational experience with solving the symmetric affine SOCCPs generated randomly by the proposed algorithm, and compare the numerical performance with the Fischer-Burmeister merit function approach [4]. Unless otherwise stated, the function f_YF in Algorithm 3.1 is always defined as in (11), where ψ_YF is defined by (9)-(10) with ψ0(t) = ¹₂(max{0, t})².

The symmetric affine SOCCP is stated as follows: finding ζ ∈ IRⁿ such that

hF (ζ), ζi = 0, ζ ∈ K, F (ζ) = Mζ + q ∈ K (18) where M ∈ IR^n×n and q ∈ IRⁿ are a given symmetric positive semidefinite matrix and a vector, respectively. In our experiments, the matrix M and the vector q are generated by the following procedure. Elements of q were chosen randomly from the interval [−1, 1]

and the matrix M was obtained by setting M = NN^T, where N is a square matrix whose nonzero elements are chosen randomly from the interval [−1, 1]. In this procedure, the number of nonzero elements of N is determined so that the nonzero density of M can be approximately estimated.

All experiments were done at a PC with 2.8GHz CPU and 512MB memory. The computer codes were all written in Matlab 6.1. To improve the numerical behavior of Algorithm 3.1, we replaced the standard Armijo-rule by the nonmonotone line search as described in [18], i.e. we computed the smallest nonnegative integer m such that

f_YF(ζ^k+ β^md(ζ^k)) ≤ W_k− σβ^2mf_YF(ζ^k) (19) where W_k is given by

W_k = max

j=k−mk,...,kf_YF(ζ^j)

and where, for given nonnegative integers ˆm and s, we set m_k =

( 0 if k ≤ s

min {m_k−1+ 1, ˆm} otherwise . (20) Throughout the experiments, unless otherwise stated, we used the following parameters:

ˆ

m = 5, s = 5, β = 0.3, and σ = 1.0e − 4. (21) For the Fischer-Burmeister merit function (FBMF, for short) approach [4], we chose a limited-memory BFGS algorithm with 5 limited-memory vector-updates [3] to solve the unconstrained minimization reformulation (8) for the SOCCP (1). For the scaling matrix H⁰ = γI in the BFGS algorithm, we adopted the choice of γ = p^Tq/p^Tq recommended by [27, P. 226], where p := ζ − ζ^oldand q := ∇f_FB(ζ) − ∇f_FB(ζ^old). To ensure convergence, we revert to the steepest descent direction −∇f_FB(ζ) whenever the current direction 4 fails to satisfy the sufficient descent condition

∇f_FB(ζ)^T4 ≤ −10⁻⁴k∇f_FB(ζ)kk4k.

In addition, we also employed the same nonmonotone line search as above to seek a suit- able step-length, except that the parameter β is chosen as 0.2.

During the experiments, we started Algorithm 3.1 and the FBMF approach with the starting point ζ⁰ = 0.001(1, 1, · · · , 1)^T and terminated the iterate once one of the following conditions is satisfied:

(1) max{Ψ(ζ), |F (ζ)^Tζ|} ≤ 10⁻⁴, where Ψ represents f_YF or f_FB. (2) The number of iteration is over 50000.

(3) The step-length is lower than 10⁻¹⁶.

We have done the following three groups of experiments.

Experiment A. Testing the influence of the scaling strategy on Algorithm 3.1 and the FBMF method. Note that, when the matrix M and the vector q in (18) are replaced by

M = M/w and¯ q = q/w¯ (22)

where w ≥ 1 is a constant, the optimal solution of problem (18) does not change. Hence, in this experiment, we generated 10 test problems with sparsity 0.5% and 10% and m = 10, n₁ = n₂ = · · · = n_m = 100, and then solved each problem and their different scaled formulations with Algorithm 3.1 and the FBMF approach. Numerical results are summarized in Tables 1 and 2, where N0. represents the number of problem, Den denotes the approximate sparsity of M, Nf and Time respectively denote the total number of function evaluations and the CPU time for solving each problem.

Table 1: Numerical results of Algorithm 3.1 for the scaled problems

w = 1 w = 10 w = 50 w = 100

N0. Den Nf Time Nf Time Nf Time Nf Time

1 0.5% 21693 69.68 6716 21.71 4201 15.31 5608 21.89 2 0.5% 55916 175.9 26234 85.81 17836 68.15 24073 98.15

3 0.5% 11897 44.57 989 3.82 803 3.40 1168 5.34

4 0.5% 14860 53.68 998 4.04 776 3.28 1047 5.07

5 0.5% 13260 48.53 553 2.01 553 2.32 733 3.20

6 10% – – 2238 89.67 237 10.09 99 4.46

7 10% – – 2518 95.54 264 10.64 114 5.21

8 10% – – 8592 344.4 228 10.26 162 7.23

9 10% – – – – 273 12.78 81 6.18

10 10% – – 1982 82.60 239 10.98 125 5.56

where “−” means that the iteration was stopped since the step-length was less than 10⁻¹⁶. From Tables 1 and 2, we see that, when w > 1, i.e. imposing the scaling strategy on the original problems, Algorithm 3.1 and the FBMF approach require much less function evaluations. Therefore, the scaling strategy in (22) can greatly improve the numerical performance of Algorithm 3.1 and the merit function approach. In particular, for those problems to which Algorithm 3.1 fails due to too small step-length, using the scaling strategy can yield satisfying solutions. This implies that Algorithm 3.1 has more depen- dence on the scaling strategy than the MF approach.

Experiment B. Testing Algorithm 3.1 and the FBMF approach on the affine SOCCP (18) with various degree of sparsity. In this experiment, we generated 10 test problems with m = 1 and n = 1000 for each nonzero density 0.1%, 0.5%, 1%, 10%, 50% and 80%, and then solved each problem with Algorithm 3.1 and the FBMF approach. Numerical results were summarized in Tables 3-4, where Nf and Time are same as Experiment A, Gap denotes the value of |F (ζ)^Tζ| at the final iteration, and Scale in Table 4 denotes the value of w in (22). In particular, the values of Gap, Nf and Time in Table 4 are the averages of 10 trials for each sparsity.

From Table 3, it appears that Algorithm 3.1 and the FBMF approach have similar nu-

Table 2: Numerical results of the FBMF method for the scaled problems

w = 1 w = 10 w = 50 w = 100

N0. Den Nf Time Nf Time Nf Time Nf Time

1 0.5% 8135 56.96 4346 30.01 3076 22.14 4649 30.56 2 0.5% 9086 57.56 14020 91.06 16619 117.5 21284 149.9

3 0.5% 611 3.95 531 3.70 812 5.89 976 6.82

4 0.5% 1030 7.65 677 5.17 493 4.09 970 6.96

5 0.5% 769 5.56 403 2.98 583 4.09 591 4.32

6 10 % 6682 488.0 807 64.15 132 9.65 100 7.40

7 10 % 4668 337.7 737 56.85 247 19.21 185 16.37

8 10 % 5639 431.1 812 63.82 131 10.12 114 9.20

9 10 % 4616 347.4 723 57.21 112 9.20 81 6.18

10 10 % 5818 452.6 702 59.12 220 17.04 96 7.59

merical performance on those problems with sparsity 0.1%. However, from Table 4, we see that, under the scaling strategy shown, Algorithm 3.1 always needed less CPU time than the FBMF approach although the former may require more function evaluations. In addition, we also observe that the number of function evaluations required by Algorithm 3.1 will become less when the sparsity of M becomes higher.

Experiment C. Testing Algorithm 3.1 and the FBMF approach on the affine SOCCP (18) with various Cartesian structures of K. To construct SOCs of various types, we chose n_i and m such that n₁ = n₂ = · · · = n_m and n₁+ · · · + n_m = 2000. For each type of K, we solved 10 test problems with nonzero density 1% by Algorithm 3.1 and the FBMF approach, respectively. Numerical results were reported in Table 5, where Scale, Gap, Nf and Time are same as Experiment A, and particularly the values of Gap, Nf and Time are the averages of 10 trials for each type of K.

From Table 5, we see that, under the scaling strategy shown, Algorithm 3.1 is comparable with the FBMF method for the first five groups of test problems whether in the CPU time or in the number of function evaluations. For the last group of test problems, Algorithm 3.1 obviously required more CPU time and function evaluations than the FBMF approach. However, from Table 6, we see that if Scale is still chosen as 100 but

Table 3: Numerical results for the affine SOCCP with sparsity 0.1%

Algorithm 3.1 MF method Algorithm 3.1 MF method

N0. Nf Time Nf Time N0. Nf Time Nf Time

1 597 0.76 369 1.10 2 * * * *

3 539 0.85 325 0.98 4 * * * *

5 * * * * 6 * * * *

7 254 0.34 127 0.33 8 * * * *

9 * * * * 10 799 0.95 143 0.28

where “∗” means that the iteration was stopped since the number of iteration was over 50000.

Table 4: Numerical results for the affine SOCCP with different sparsity Algorithm 3.1 MF approach

Den Scale Gap Nf time Gap Nf time

0.5% 1 9.16e-5 1597.1 3.28 8.77e-5 914.9 4.09 1% 1 8.10e-5 7401.4 28.76 6.51e-5 5016.8 38.31 10% 10 6.89e-5 402.6 15.37 6.85e-5 312.0 21.92 50% 100 5.90e-5 472.2 17.07 6.59e-5 568.7 39.00 80% 100 4.81e-5 468.4 17.54 7.62e-5 668.7 45.23

the parameter β in the line search is chosen as 0.1 instead of 0.3, the numerical per- formance of Algorithm 3.1 will have a great improvment, and moreover, the CPU time and the number of function evaluations needed are comparable with those of the FBMF method.

To sum up, for the symmetric affine SOCCPs in (18), if a suitable scaling strategy and the parameter β are used, Algorithm 3.1 will be comparable with, even superior to, the FBMF method in the CPU time for solving test problems although the former may

require more function evaluations. Otherwise, the FBMF approach will be superior to Algorithm 3.1 whether in the CPU time or in the number of function evaluations.

Table 5: Numerical results for the affine SOCCP with different K Algorithm 3.1 MF approach

m Scale Gap Nf Time Gap Nf Time

1 10 7.95e-5 201.3 7.64 8.759e-5 167.5 10.58 10 10 8.24e-5 497.6 17.68 8.74e-5 217.1 14.03 50 10 9.34e-5 1193 49.12 9.71e-5 266.5 19.42 100 100 5.75e-5 116.4 6.62 7.86e-5 138.7 11.44 200 100 4.54e-5 129.2 9.28 7.93e-5 149.2 15.69 500 100 7.09e-5 9719.4 1115.4 7.93e-5 149.2 15.69

Table 6: Numerical results of Algorithm 3.1 for different β

β = 0.3 β = 0.1 β = 0.3 β = 0.1

N0. Nf Time Nf Time N0. Nf Time Nf Time

1 75 8.42 198 21.90 2 3081 355.5 142 15.23

3 19001 2231.4 210 22.25 4 6644 769.0 179 19.09

5 69 7.93 178 19.00 6 75 8.37 804 88.67

7 62169 7068.4 132 14.56 8 5939 689.6 295 32.23

9 77 8.54 208 23.03 10 64 7.12 144 16.03

5 Final Remarks

In this paper, we investigated a descent method for the equivalent reformulation (11) of the SOCCP which was also used for the NCP and the SDCP in literature, and proved

its global convergence under some mild assumptions. Numerical comparision with the Fischer-Burmeister merit function approach [4] for symmetric affine SOCCPs generated randomly indicate that the descent method is comparable with, even to superior to, the FBMF approach in the CPU time if a suitable scaling strategy and the parameter β in line search are adopted. We also expect that the method can be used to deal with large SOCCPs due to very small computational work per iteration. In addition, we notice that the proposed algorithm does not work for another reformulation (8) of the SOCCP since f_FB lacks property of bounded level sets (Lem. 2.3) where ψ₀ plays an important role therein.

Prop. 3.1-3.2 are more or less an afterthought of [4], nonetheless, it does parallel the extension to the SOCCP from the NCP and SDCP cases. On the other hand, this work does a further study based on replacing the conditions of monotonicity and strict feasibil- ity by a new (and weaker under certain sense) so-called R01-function. More specifically, under the new so-called R₀₁-function condition, the level sets of f_YF are still bounded and the proposed descent algorithm still has global convergence. These results are significant not only they are brand-new but also there is no needs the assumption of strict feasibility therein.

One future topic is to analyze the convergence rate theoretically which is more in- tractable. Other direction like weakening conditions which guarantees the property of bounded level sets is also interesting and worthwhile. There may have the direction as one referee pointed out which is to apply this optimization method to real-life studies, for example [5] and references therein.

Acknowledgement. The authors thank the referees for their careful reading of the paper and helpful suggestions.

References

[1] F. Alizadeh and S. Schmieta (2000), Symmetric cones, potential reduction meth- ods, and word-by-word extensions, in Handbook of Semidefinite Programming, edited by H. Wolkowicz, R. Saigal, and L. Vandenberghe, Kluwer Academic Publishers, Boston, pp. 195–233.

[2] D. P. Bertsekas (1999), Nonlinear Programming, 2nd edition, Athena Scientific, Belmont.

[3] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu (1995), A limited memory algorithm for bound constrained optimization, SIAM Journal of Scientific Computing, vol. 16, pp. 1190-1208.

[4] J.-S. Chen and P. Tseng (2005), An unconstrained smooth minimization refor- mulation of the second-order cone complementarity problem, Mathematical Program- ming, vol. 104, pp. 293-327.

[5] C-T. Cheng and K-W. Chau (2001), Fuzzy iteration methodology for reservoir flood control operation, Journal of the American Water Resources Association, vol.

37, pp. 1381-1388.

[6] X.-D. Chen, D. Sun, and J. Sun(2003), Complementarity functions and numerical experiments for second-order cone complementarity problems, Computational Opti- mization and Applications, vol. 25, pp. 39-56.

[7] F. Facchinei and J.-S. Pang(2003), Finite-Dimensional Variational Inequalities and Complementarity Problems, Volumes I and II, Springer-Verlag, New York.

[8] F. Facchinei and J. Soares(1997), A new merit function for nonlinear comple- mentarity problems and a related algorithm, SIAM Journal on Optimization, vol. 7, pp. 225-247.

[9] U. Faraut and A. Kor´anyi(1994), Analysis on Symmetric Cones, Oxford Math- ematical Monographs, Oxford University Press, New York.

[10] M. C. Ferris and J.-S. Pang(1997), Engineering and economic applications of complementarity problems, SIAM Review, vol. 39, pp. 669-713.

[11] M. C. Ferris and J.-S. Pang(1996), editors, Complementarity and Variational Problems: State of the Art, SIAM Publications, Philadelphia.

[12] M. C. Ferris, O. L. Mangasarian, and J.-S. Pang(2001), editors, Comple- mentarity: Applications, Algorithms and Extensions, Kluwer Academic Publishers, Dordrecht.

[13] A. Fischer(1992), A special Newton-type optimization methods, Optimization, vol.

24, pp. 269-284.

[14] A. Fischer(1997), Solution of the monotone complementarity problem with locally Lipschitzian functions, Mathematical Programming, vol. 76, pp. 513-532.

[15] R. Fletcher(1987), Practical Methods of Optimization, 2nd edition, Wiley- Interscience, Chichester.

[16] M. Fukushima, Z.-Q. Luo, and P. Tseng(2002), Smoothing functions for second-order cone complementarity problems, SIAM Journal on Optimization, vol.12, pp. 436-460.

[17] C. Geiger and C. Kanzow(1996), On the resolution of monotone complemen- tarity problems, Computational Optimization and Applications, vol.5, pp. 155-173.

[18] L. Grippo, F. Lampariello and S. Lucidi(1986), A nonmonotone line search technique for Newton’s method, SIAM Journal on Numerical Analysis, 1986, vol.23, pp. 707-716.

[19] M. Seetharama Gowda, Roman Sznajder, and J. Tao(2004), Some P- properties for the linear transformations on Euclidean Jordan algebras, Linear Algebra and its Applications, vol.393, pp. 203-232.

[20] S. Hayashi, N. Yamashita, and M. Fukushima(2001), On the coerciveness of merit functions for the second-order cone complementarity problem, Report, Depart- ment of Applied Mathematics and Physics, Kyoto University, Kyoto, Japan.

[21] S. Hayashi, T. Yamaguchi, N. Yamashita, and M. Fukushima(2005), A ma- trix splitting method for symmetric affine second-order cone complementarity problem, Journal of Computational and Applied Mathematics, vol.175, pp. 335-353.

[22] C. Kanzow(1994), An unconstrained optimization technique for large scale linearly constrained convex minimization problems, Computing, vol. 53, pp. 101-117.

[23] C. Kanzow(1996), Nonlinear complementarity as unconstrained optimization, Journal of Optimization Theory and Applications, vol. 88, pp. 139-155.

[24] Y.-J. Liu, Z.-W. Zhang, and Y.-H. Wang, Some properties of a class of merit functions for symmetric cone complementarity problems, to appear in Asia-Pacific Journal of Operational Research, 2005.

[25] M. S. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret(1998), Application of second-order cone programming, Linear Algebra and its Applications, vol. 284, pp.

193-228.

[26] M. Malik and S.R. Mohan, On Q and R₀ properties of a quadratic representation in the linear complementarity problems over the second-order cone, Linear Algebra and its Applications, vol. 397, pp. 85-97, 2005.

[27] J. Nocedal and S. J. Wright(1999), Numerical Optimization, Springer-Verlag, New York.

[28] R. D. C. Monteiro and T. Tsuchiya(2000), Polynomial convergence of primal- dual algorithms for the second-order cone programs based on the MZ-family of direc- tions, Mathematical Programming, vol. 88, pp. 61-83.

[29] S. Schmieta and F. Alizadeh(2001), Associative and Jordan algebras, and poly- nomial time interior-point algorithms for symmetric cones, Mathematics of Opera- tions Research, vol. 26, pp. 543-564.

[30] J. Tao and M. S. Gowda, Some P -properties for the nonlinear transformations on Euclidean Jordan Algebra, Technical Report, Department of Mathematics and Statistics, University of Maryland, 2004.

[31] T. Tsuchiya(1999), A convergence analysis of the scaling-invariant primal-dual path-following algorithms for second-order cone programming, Optimization Methods and Software, vol. 11, pp. 141–182.

[32] R. J. Vanderbei and H. Y. Benson(1999), On formulating semidefinite pro- gramming problems as smooth convex nonlinear optimization problems, ORFE 99-01, Department of Operations Research and Financial Engineering, Princeton University, Princeton.

[33] N. Yamashita and M. Fukushima(1999), A new merit function and a descent method for semidefinite complementarity problems, in Reformulation – Nonsmooth, Piecewise Smooth, Semismooth and Smoothing Methods edited by M. Fukushima and L. Qi, Kluwer Academic Publishers, Boston, 405–420.