to appear in Journal of Applied Mathematics and Computing, 2016
Numerical comparisons based on four smoothing functions for absolute value equation
B. Saheya1
College of Mathematical Science Inner Mongolia Normal University Hohhot 010022, Inner Mongolia, P. R. China
E-mail: saheya@imnu.edu.cn
Cheng-He Yu
Department of Mathematics National Taiwan Normal University
Taipei 11677, Taiwan E-mail: 60240031S@ntnu.edu.tw
Jein-Shan Chen 2 Department of Mathematics National Taiwan Normal University
Taipei 11677, Taiwan E-mail: jschen@math.ntnu.edu.tw
March 25, 2016
(revised on August 18, 2016)
Abstract The system of absolute value equation, denoted by AVE, is a non-differentiable NP-hard problem. Many approaches have been proposed during the past decade and most of them focus on reformulating it as complementarity problem and then solve it accordingly. Another approach is to recast the AVE as a system of nonsmooth equations and then tackle with the nonsmooth equations. In this paper, we follow this path. In particular, we rewrite it as a system of smooth equations and propose four new smoothing functions along with a smoothing-type algorithm to solve the system of equations. The
1The author’s work is supported by Natural Science Foundation of Inner Mongolia (Award Number:
2014MS0119).
2Corresponding author. The author’s work is supported by Ministry of Science and Technology, Taiwan.
main contribution of this paper focuses on numerical comparisons which suggest a better choice of smoothing function along with the smoothing-type algorithm.
Key words. Smoothing function, smoothing algorithm, singular value, convergence.
1 Introduction
The absolute value equation (AVE) is in the form of
Ax + B|x| = b, (1)
where A ∈ Rn×n, B ∈ Rn×n, B 6= 0, and b ∈ Rn. Here |x| means the componentwise absolute value of vector x ∈ Rn. When B = −I, where I is the identity matrix, the AVE (1) reduces to the special form:
Ax − |x| = b. (2)
It is known that the AVE (1) was first introduced by Rohn in [18] and recently has been investigated by many researchers, for example, Hu and Huang [5], Jiang and Zhang [7], Ketabchi and Moosaei [8], Mangasarian [9–13], Mangasarian and Meyer [14], Prokopyev [15], and Rohn [20].
In particular, Mangasarian and Meyer [14] show that the AVE (1) is equivalent to the bilinear program, the generalized LCP (linear complementarity problem), and the stan- dard LCP provided 1 is not an eigenvalue of A. With these equivalent reformulations, they also show that the AVE (1) is NP-hard in its general form and provide existence results. Prokopyev [15] further improves the above equivalence which indicates that the AVE (1) can be equivalently recast as LCP without any assumption on A and B, and also provides a relationship with mixed integer programming. In general, if solvable, the AVE (1) can have either unique solution or multiple (e.g., exponentially many) solutions.
Indeed, various sufficiency conditions on solvability and non-solvability of the AVE (1) with unique and multiple solutions are discussed in [14, 15, 19]. Some variants of the AVE, like the absolute value equation associated with second-order cone and the absolute value programs, are investigated in [3] and [21], respectively.
As for its numerical solvers, many numerical methods for solving the AVEs (1)–(2) have been proposed. A parametric successive linearization algorithm for the AVE (1) that terminates at a point satisfying necessary optimality conditions is studied in [10].
The generalized Newton algorithm for the AVE (2) is investigated in [11], in which it was proved that this algorithm converges linearly from any starting point to the unique solution of the AVE (2) under the condition that kA−1k < 14. The generalized Newton algorithm with semismooth and smoothing Newton steps combined into the algorithm
is considered in [23]. The smoothing-type algorithms for solving the AVEs (1)-(2) are studied in [1, 5, 7]. A branch and bound method for the absolute value programs (AVP), which is an extension of the AVE, is studied in [21].
Among the aforementioned approaches, many of them focus on reformulating it as complementarity problem and then solve it accordingly. An alternative approach is to recast the AVE as a system of nonsmooth equations and then tackle with the nonsmooth equations by applying nonsmooth Newton algorithm [16] or smoothing Newton algorithm [17]. In this paper, we follow the latter pathway. More specifically, we rewrite it as a system of smooth equations and propose four new smoothing functions along with a smoothing-type algorithm to solve the system of equations. To see this, motivated by the approach in [1, 7], we define Hi : Rn+1 → Rn+1 as
Hi(µ, x) =
µ
Ax + BΦi(µ, x) − b
for µ ∈ R and x ∈ Rn (3) where Φi : Rn+1 → Rn is given by
Φi(µ, x) :=
φi(µ, x1) φi(µ, x2)
... φi(µ, xn)
for µ ∈ R and x ∈ Rn (4)
with four various smoothing functions φi : R2 → R that will be introduced later. The role of φi looks similar to the function φp used in [7]. However, they are substantially different. More specifically, the function φp employed in [7] is strongly semismooth on R2, whereas each φi proposed in this paper is continuously differentiable on R2. Now, we present the exact form for each function φi, which is defined as below:
φ1(µ, t) = µh
ln(1 + e−µt) + ln(1 + eµt)i
(5)
φ2(µ, t) =
t if t ≥ µ2, t2
µ + µ
4 if −µ2 < t < µ2,
−t if t ≤ −µ2.
(6)
φ3(µ, t) = p
4µ2+ t2 (7)
φ4(µ, t) =
t2
2µ if |t| ≤ µ,
|t| − µ2 if |t| > µ.
(8)
Some of the smoothing functions have appeared in other contexts for other optimization problems, but they are all novel ones for dealing with the AVE (1). The main idea in this paper is showing that the AVE (1) has a solution if and only if Hi(µ, x) = 0, φi
is continuously differentiable at any (µ, t) ∈ R++× R, and limµ↓0φi(µ, x) = |x|. Then, with these four new smoothing functions, we consider the smoothing-type algorithm studied in [6, 22] to solve the AVE (1). In other words, we reformulate the AVE (1) as parameterized smooth equations and then employ a smoothing-type algorithm to solve it. In addition, we show that the algorithm is well-defined under the assumption that the minimal singular value of the matrix A is strictly greater than the maximal singular value of the matrix B. We also show that the proposed algorithm is globally and locally quadratically convergent no matter which smoothing function φi is used. Numerical implementations and comparisons based on these four different φi are reported as well.
From the numerical results, we conclude that φ2 is the best choice of smoothing function when we apply the proposed smoothing-type algorithm. More detailed reports will be seen in Section 4.
2 Smoothing Reformulation
In this section, we depict the graphs of φifor i = 1, 2, 3, 4 and investigate their properties.
Then, we show the equivalent reformulation that Hi(µ, x) = 0 if and only if x solves the AVE (1), and talk about the condition to guarantee the unique solvability of the AVE (1). We begin with showing the pictures of φi for i = 1, 2, 3, 4, see Figure 1.
|t|
ϕ1(μ,t) ϕ2(μ,t) ϕ3(μ,t) ϕ4(μ,t)
-1.0 -0.5 0.0 0.5 1.0
0.0 0.2 0.4 0.6 0.8 1.0
t ϕi(μ,t)
Figure 1: Graphs of |t| and all four φi(µ, t) with µ = 0.1.
From Figure 1, we see that φ2 is the one which best approximates the function |t|
under the sense that it is closest to |t| among all φi for i = 1, 2, 3, 4. To see this, we adopt the max norm to measure the distance of two real-valued functions. In other words, for given two real-valued functions f and g, the distance between them is defined as
kf − gk∞ = max
t∈R{f (t) − g(t)}.
Now, for any fixed µ > 0, we know that lim
|t|→∞
φi(µ, t) − |t|
= 0, for i = 1, 2, 3.
This implies that
maxt∈R
φi(µ, t) − |t|
= |φi(µ, 0)| , for i = 1, 2, 3.
Since, φ1(µ, 0) = (2 ln 2)µ ≈ 1.4µ, φ2(µ, 0) = µ4, and φ3(µ, 0) = 2µ, we obtain φ1(µ, t) − |t|
∞ = (2 ln 2)µ ≈ 1.4µ φ2(µ, t) − |t|
∞ = µ 4 φ3(µ, t) − |t|
∞ = 2µ On the other hand, we see that
t→∞lim
φ4(µ, t) − |t|
= µ
2 and φ4(µ, 0) = 0, which says
maxt∈R
φ4(µ, t) − |t|
= µ 2. Hence, we obtain
φ4(µ, t) − |t|
∞ = µ 2. From all the above, we conclude that
φ3(µ, t) − |t|
∞ >
φ1(µ, t) − |t|
∞>
φ4(µ, t) − |t|
∞>
φ2(µ, t) − |t|
∞. (9) This shows that φ2 is the function among φi, i = 1, 2, 3, 4 which best approximates the function |t|. In fact, for fixed µ > 0, there has the local behavior that
φ3(µ, t) > φ1(µ, t) > φ2(µ, t) > |t| > φ4(µ, t). (10) A natural question arises here, does the smoothing algorithm based on φ2 perform best among all φ1, φ2, φ3, φ4? This will be answered in Section 4.
Proposition 2.1. Let φi : R2 → R for i = 1, 2, 3, 4 be defined as in (5), (6), (7) and (8), respectively. Then, we have
(a) φi is continuously differentiable at (µ, t) ∈ R++× R;
(b) lim
µ↓0φi(µ, t) = |t| for any t ∈ R.
Proof. (a) In order to prove the continuous differentiability of φi, we need to write out the expressions of ∂φi∂t(µ,t) and ∂φi∂µ(µ,t); and then show the continuity of ∂φi∂t(µ,t) and ∂φi∂µ(µ,t). (i) For i = 1, we compute that
∂φ1(µ, t)
∂t = 1
1 + e−µt
− 1
1 + eµt ,
∂φ1(µ, t)
∂µ = h
ln(1 + e−µt) + ln(1 + eµt)i + t
µ
−1 1 + e−µt
+ 1
1 + eµt
.
Then, it is clear to see that ∂φ1∂t(µ,t) and ∂φ1∂µ(µ,t) are continuous. Hence, φ1 is continuously differentiable.
(ii) For i = 2, we compute that
∂φ2(µ, t)
∂t =
1 if t ≥ µ2,
2t
µ if −µ2 < t < µ2,
−1 if t ≤ −µ2.
∂φ2(µ, t)
∂µ =
0 if t ≥ µ2,
−
t µ
2
+14 if −µ2 < t < µ2, 0 if t ≤ −µ2.
Then, it can be verified that ∂φ2∂t(µ,t) and ∂φ2∂µ(µ,t) are continuous because
lim
t→µ2
∂φ2(µ, t)
∂t = lim
t→µ2
2t µ = 1, lim
t→−µ2
∂φ2(µ, t)
∂t = lim
t→−µ2
2t
µ = −1.
and
lim
t→µ2
∂φ2(µ, t)
∂µ = lim
t→µ2
"
− t µ
2
+ 1 4
#
= 0,
lim
t→−µ2
∂φ2(µ, t)
∂µ = lim
t→−µ2
"
− t µ
2
+ 1 4
#
= 0.
Hence, φ2 is continuously differentiable.
(iii) For i = 3, we compute that
∂φ3(µ, t)
∂t = t
p4µ2+ t2,
∂φ3(µ, t)
∂µ = 4µ
p4µ2+ t2.
Again it is clear to see that ∂φ3∂t(µ,t) and ∂φ3∂µ(µ,t) are continuous. Hence, φ3 is continuously differentiable.
(iv) For i = 4, we compute that
∂φ4(µ, t)
∂t =
1 if t > µ,
t
µ if −µ ≤ t ≤ µ,
−1 if t < −µ.
∂φ4(µ, t)
∂µ =
−12 if t > µ,
−12
t µ
2
if −µ ≤ t ≤ µ,
−12 if t < −µ.
Then, we conclude that ∂φ4∂t(µ,t) and ∂φ4∂µ(µ,t) are continuous by checking
t→µlim
∂φ4(µ, t)
∂t = lim
t→µ
t µ = 1,
t→−µlim
∂φ4(µ, t)
∂t = lim
t→−µ
t
µ = −1.
and
limt→µ
∂φ4(µ, t)
∂µ = lim
t→µ
"
−1 2 × t
µ
2#
= −1 2,
t→−µlim
∂φ4(µ, t)
∂µ = lim
t→−µ
"
−1 2 × t
µ
2#
= −1 2. Hence, φ4 is continuously differentiable.
From all the above, we prove that φi is continuously differentiable at (µ, t) ∈ R++× R.
(b) For i = 1, 2, 3, 4, we always have the following:
µ→0lim
∂φi(µ, t)
∂t =
1 if t > 0,
−1 if t < 0, which verifies part (b). 2
For subsequent needs in convergence analysis and numerical implementations, we summarize the gradient of each φi as below.
∇φ1(µ, t) =
h
ln(1 + e−µt) + ln(1 + eµt)i +µt
−1 1+e− tµ
+ 1
1+eµt
1 1+e− tµ
− 1
1+e
t µ
.
∇φ2(µ, t) = ξ1 ξ2
, where ξ1 =
0 if t ≥ µ2,
−
t µ
2
+14 if −µ2 < t < µ2, 0 if t ≤ −µ2.
ξ2 =
1 if t ≥ µ2,
2t
µ if −µ2 < t < µ2,
−1 if t ≤ −µ2.
∇φ3(µ, t) =
√ 4µ 4µ2+t2
√ t 4µ2+t2
.
∇φ4(µ, t) = v1 v2
, where v1 =
−12 if t > µ,
−12
t µ
2
if −µ ≤ t ≤ µ,
−12 if t < −µ.
v2 =
1 if t > µ,
t
µ if −µ ≤ t ≤ µ,
−1 if t < −µ.
In fact, Proposition 2.1 can be also depicted by geometric views. In particular, from Figures 2-5, we see that when µ ↓ 0, φi is getting closer to |t|, which verifies Proposition 2.1(b).
μ=0.5 μ=0.3 μ=0.1 μ=0.01
-2 -1 0 1 2
0.0 0.5 1.0 1.5 2.0
t ϕ1(μ,t)
Figure 2: Graphs of φ1(µ, t) with µ = 0.01, 0.1, 0.3, 0.5.
μ=0.5 μ=0.3 μ=0.1 μ=0.01
-2 -1 0 1 2
0.0 0.5 1.0 1.5 2.0
t ϕ2(μ,t)
Figure 3: Graphs of φ2(µ, t) with µ = 0.01, 0.1, 0.3, 0.5.
μ=0.5 μ=0.3 μ=0.1 μ=0.01
-2 -1 0 1 2
0.0 0.5 1.0 1.5 2.0
t ϕ3(μ,t)
Figure 4: Graphs of φ3(µ, t) with µ = 0.01, 0.1, 0.3, 0.5.
μ=0.5 μ=0.3 μ=0.1 μ=0.01
-2 -1 0 1 2
0.0 0.5 1.0 1.5 2.0
t ϕ4(μ,t)
Figure 5: Graphs of φ4(µ, t) with µ = 0.01, 0.1, 0.3, 0.5.
Now, in light of Proposition 2.1, we obtain the equivalent reformulation Hi(µ, x) = 0 for the AVE (1).
Proposition 2.2. Let Φi(µ, x) for i = 1, 2, 3, 4 be defined as in (4). Then, we have
(a) Hi(µ, x) = 0 if and only if x solves the AVE (1);
(b) Hi is continuously differentiable on Rn+1\ {0} with the Jacobian matrix given by
∇Hi(µ, x) :=
1 0
B ∇1Φi(µ, x) A + B ∇2Φi(µ, x)
(11) where
∇1Φi(µ, x) :=
∂φi(µ,x1)
∂µ
∂φi(µ,x2)
∂µ...
∂φi(µ,xn)
∂µ
,
∇2Φi(µ, x) :=
∂φi(µ,x1)
∂x1 0 · · · 0
0 ∂φi∂x(µ,x2)
2 · · · 0
... ... . .. ... 0 · · · 0 ∂φi∂x(µ,xn)
n
.
Proof. This result follows from Proposition 2.1 immediately and the computation of the Jacobian matrix is straightforward. 2
For completeness, we also talk about the unique solvability of the AVE (1), which is presumed in our numerical implementations. The following assumption and proposition are both employed from [7]. The Assumption 2.3 will be also used to guarantee that
∇Hi(µ, x) is invertible at any (µ, x) ∈ R++× Rn, see Proposition 3.2 in Section 3.
Assumption 2.3. The minimal singular value of the matrix A is strictly greater than the maximal singular value of the matrix B.
Proposition 2.4. ([7, Prop. 2.3]) The AVE (1) is uniquely solvable for any b ∈ Rn if Assumption 2.3 is satisfied.
3 A smoothing-type algorithm
From Proposition 2.2, we know that the AVE (1) is equivalent to Hi(µ, x) = 0. Ac- cordingly, in this section, we consider the smoothing-type algorithm as in [1, 7] to solve Hi(µ, x) = 0. In fact, this type of algorithm has been also proposed for solving other kinds of problems, see [2, 6, 22] and references therein.
Algorithm 3.1. (A Smoothing-Type Algorithm)
Step 0. Choose δ, σ ∈ (0, 1), µ0 > 0, x0 ∈ Rn. Set z0 := (µ, x0). Denote e0 := (1, 0) ∈ R × Rn. Choose β > 1 such that (min {1, kHi(z0)k})2 ≤ βµ0. Set k := 0.
Step 1. If kHi(zk)k = 0, stop.
Step 2. Set τk := min1, kHi(zk)k , and compute 4zk := (4µk, 4xk) ∈ R × Rn by using
∇Hi(zk)4zk = −Hi(zk) + (1/β)τk2e0, (12) where ∇Hi(·) is defined as in (11).
Step 3. Let αk be the maximum of the values 1, δ, δ2, · · · such that
kHi(zk+ αk4zk)k ≤ [1 − σ(1 − 1/β)αk] kHi(zk)k (13) Step 4. Set zk+1 := zk+ αk4zk and k := k + 1. Back to Step 1.
Following the same arguments as in [4, 6], the line search (13) in the above scheme is well-defined. In other words, the Algorithm 3.1 is well-defined and possesses some nice properties.
Proposition 3.2. (a) Suppose that Assumption 2.3 holds. Then, the Algorithm 3.1 is well-defined.
(b) Let the sequence zk be generated by Algorithm 3.1. Then, (i) both kHi(zk)k and {τk} are monotonically decreasing;
(ii) τk2 ≤ βµk holds for all k;
(iii) the sequence {µk} is monotonically decreasing, and µk > 0 for all k.
Proof. Please refer to [6, Remark 2.1] or [7, Prop. 3.1]. 2
The key point in the above scheme is the solvability of Newton equations (12) in Step 2. The following result is regarding this issue. Since the Φi function plays almost the same role as the function Φp used in [7], the below Proposition 3.3 can be obtained by mimicking the same arguments as in [7, Theorem 3.2]. We omit its proof and only state it.
Proposition 3.3. Let Hi and ∇Hi be given as in (3) and (11), respectively. Suppose that Assumption 2.3 holds. Then, ∇Hi(µ, x) is invertible at any (µ, x) ∈ R++× Rn.
Next, we discuss the global and local convergence. Again, although the function Φi
here is continuously differentiable and the function Φp used in [7] is only semismooth, their roles in the proof are almost the same. Consequently, the arguments for convergence analysis are almost the same. Hence, we also omit the detailed proof and only present the convergence result.
Proposition 3.4. Suppose that Assumption 2.3 holds and that the sequence zk is generated by Algorithm 3.1. Then,
(a) zk is bounded;
(b) any accumulation point of zk is a solution of the AVE (1).
(c) The whole sequence zk convergence to z∗ with kzk+1 − zkk = o kzk− z∗k and µk+1 = µ2k.
4 Numerical implementations
In this section, we report the numerical results of Algorithm 3.1 for solving the AVE (1) and (2). All numerical experiments are carried out in Mathematica 10.0 running on a PC with Intel i5 of 3.00GHz CPU processor, 4.00GB Memory and 32-bit Windows 7 operating system.
In our numerical experiments, the stoping criteria for Algorithm 3.1 is kHi(zk)k ≤ 1.0e − 6. We also stop programs when the total iteration is more than 100. Throughout the computational experiments, the following parameters are used:
δ = 0.5, σ = 0.0001, µ0 = 0.1, β = max1, 1.01 ∗ τ02/µ .
4.1 Experiments on the AVE Ax − |x| = b
In this subsection we consider the simplified form of AVE (2). Consider the ordinary differential equation [24, Example 4.2]:
d2x
dt2 − |x| = (1 − t2), x(0) = −1, x(1) = 0, t ∈ [0, 1]. (14) As explained in [24, Example 4.2], after descretization (by using finite difference method), the above ODE can be recast an AVE in form of
Ax − |x| = b, (15)
where the matrix A is given by
ai,j =
−242, i = j, 121, |i − j| = 1, 0, otherwise.
(16)
We implement the above problems by using φi, i = 1, 2, 3, 4 and n = 2, 5, 10, 20, . . . , 100, respectively. Every starting point x is randomly generated 10 times from a uniform distribution on x ∈ [−2, 2]. The results are put together in Table 1, where
Table 1: The numerical results of ordinary differential equation (14)
Dim N φ1 T φ1 Ar φ1 N φ2 T φ2 Ar φ2 N φ3 T φ3 Ar φ3 N φ4 T φ4 Ar φ4
2 5.1 0.0967 3.30E-07 3.9 0.0015 6.92E-08 5.1 0.0016 5.93E-08 4 0.0062 5.99E-08 5 5.9 0.3697 2.23E-07 4.1 0.0031 7.47E-08 5.6 0.0062 2.21E-08 4.2 0.0016 6.54E-08 10 6.4 0.4851 2.98E-07 4.3 0.0094 2.10E-07 5.9 0.0031 1.05E-07 4.5 0.0031 4.67E-08 20 5.2 0.4290 2.41E-07 4.9 0.0078 1.10E-08 6.3 0.0078 2.13E-09 5 0.0094 2.46E-09 40 8.8 4.4117 4.66E-07 6.1 0.5210 5.28E-08 7.3 0.0172 6.59E-08 6.3 0.0156 1.88E-07 60 9.1 2.4289 2.31E-07 6.8 0.0281 4.49E-08 9 0.0312 1.20E-08 7.7 0.0312 1.31E-07 80 9.8 2.0514 3.61E-07 7.4 0.0374 3.21E-10 9.3 0.0452 3.21E-08 9.2 0.0593 3.15E-08 100 9.8 8.2306 4.44E-07 7.8 0.0577 8.78E-08 10 0.0671 2.26E-07 9.5 0.0827 2.83E-08
Dim denotes the size of problem,
N φi denotes the average number of iterations,
T φi denotes the average value of the CPU time in seconds,
Ar φi denotes the average value ofkH(zk)k when Algorithm 3.1 stop.
From Table 1, in terms of the average number of iterations, the efficiency of φ2(µ, t) is best, followed by φ4(µ, t), φ3(µ, t) and φ1(µ, t). This is especially true for the problem of high dimension ordinary differential equation (14). In terms of time efficiency, φ1(µ, t) is still better than other functions too. In other words, for the AVE (2) arising from the ODE (15), we have
φ2(µ, t) > φ4(µ, t) > φ3(µ, t) > φ1(µ, t) where “>” means “better performance”.
To compare the performance of smoothing function φi(µ, t), i = 1, 2, 3, 4, we adopt the performance profile which is introduced in [25] as a means. In other words, we regard Algorithm 3.1 corresponding to a smoothing function φi(µ, t), i = 1, 2, 3, 4 as a solver, and assume that there are ns solvers and np test problems from the test set P which is generated randomly. We are interested in using the iteration number as performance measure for Algorithm 3.1 with different φi(µ, t). For each problem p and solver s, let
fp,s = iteration number required to solve problem p by solver s.
We employ the performance ratio
rp,s:= fp,s
min{fp,s: s ∈ S},
where S is the four solvers set. We assume that a parameter rp,s ≤ rM for all p, s are chosen, and rp,s = rM if and only if solver s does not solve problem p. In order to obtain an overall assessment for each solver, we define
ρs(τ ) := 1
npsize{p ∈ P : rp,s≤ τ },
which is called the performance profile of the number of iteration for solver s. Then, ρs(τ ) is the probability for solver s ∈ S that a performance ratio fp,s is within a factor τ ∈ R of the best possible ratio.
We then need to test the four functions for ODE (14) at random starting points.
In particular, starting points for each dimension are randomly chosen 20 times from a uniform distribution on x ∈ [−2, 2]. In order to obtain an overall assessment for the four functions, we are interested in using the number of iterations as a performance measure for Algorithm 3.1 with φ1(µ, t), φ2(µ, t), φ3(µ, t), and φ4(µ, t), respectively. The performance plot based on iteration number is presented in Figure 6. From this figure, we see that φ2(µ, t) working with Algorithm 3.1 has the best numerical performance, followed by φ4(µ, t). In other words, in view of “iteration numbers”, there has
φ2(µ, t) > φ4(µ, t) > φ3(µ, t) > φ1(µ, t) where “>” means “better performance”.
1.0 1.5 2.0 2.5 3.0
0.0 0.2 0.4 0.6 0.8 1.0
τ
ρs(τ) ϕ1(μ,t)
ϕ2(μ,t) ϕ3(μ,t) ϕ4(μ,t)
Figure 6: Performance profile of iteration numbers of Algorithm 3.1 for the ODE (14).
We are also interested in using the computing time as performance measure for Algo- rithm 3.1 with different φi(µ, t), i = 1, 2, 3, 4. The performance plot based on computing time is presented in Figure 7. From this figure, we can also see the function φ2(µ, t) has best performance, then followed by φ3(µ, t). Note that the time efficiency of φ1(µ, t) is very bad. In other words, in view of “computing time”, there has
φ2(µ, t) > φ3(µ, t) > φ4(µ, t) > φ1(µ, t) where “>” means “better performance”.
In summary, for the special AVE (2) arising from the ODE (14), no matter the number of iterations or the computing time is taken into account, the function φ2(µ, t) is the best choice for the Algorithm 3.1.
2 4 6 8 10 12 14 0.0
0.2 0.4 0.6 0.8 1.0
τ
ρs(τ) ϕ1(μ,t)
ϕ2(μ,t) ϕ3(μ,t) ϕ4(μ,t)
Figure 7: Performance profile of computing time of Algorithm 3.1 for the ODE (14).
4.2 Experiments on the general AVE Ax + B|x| = b
In this subsection we consider the general AVE (1): Ax + B|x| = b. Here matrix A (or B) is equal to a normal distribution random matrix minus another one so that we can randomly generate the testing problems.
In order to ensure that Assumption 2.3 holds, we further modify the matrix A in light of the below conditions.
• If min{wii : i = 1, . . . , n} = 0 with {u, w, v} = SingularValueDecomposition[A], then we set A = u(w + 0.01 ∗ IdentityMatrix[n])v.
• Set A = λmaxλ (BTB)+0.01
min(ATA) A.
Then, it is clear to verify that Assumption 2.3 is satisfied for such A. Moreover, we set p =2RandomVariate[NormalDistribution[ ],{n, 1}] and b = Ap + B|p| so that the testing problems are solvable.
We implement the above problems for φi, i = 1, 2, 3, 4 and n = 2, 5, 10, 20, . . . , 100, re- spectively. Every case is randomly generated 10 times for testing. The numerical results are listed in Table 2. From Table 2, in terms of the number of iterations and computa- tion time, the efficiency of φ2(µ, t) is best, followed by φ4(µ, t). The iteration number of φ1(µ, t) is less than φ3(µ, t), but the computing time of φ1(µ, t) is more than φ3(µ, t).
Figure 8 shows the performance profile of iteration number in Algorithm 3.1 in the range of τ ∈ [1, 15] for four solvers on 100 test problem which are generated randomly.
The four solvers correspond to Algorithm 3.1 with φ1(µ, t), φ2(µ, t), φ3(µ, t), and φ4(µ, t), respectively. From this figure, we see that φ2(µ, t) working with Algorithm 3.1 has the best numerical performance, followed by φ4(µ, t). In summary, from the viewpoint of
Table 2: The numerical results of experiments
Dim N φ1 T φ1 Ar φ1 N φ2 T φ2 Ar φ2 N φ3 T φ3 Ar φ3 N φ4 T φ4 Ar φ4
2 6.2 0.4596 5.00E-7 3.6 0.0031 8.56E-8 7.1 0.0016 1.79E-7 3.9 0 8.04E-8 5 7.4 0.2246 6.05E-7 4.1 0.0031 8.39E-8 9.6 0.0094 4.73E-7 4.3 0.0016 7.53E-8 10 10.2 1.0733 2.23E-7 4.3 0.0062 8.26E-8 17.2 0.0187 4.79E-7 4.7 0.0031 7.53E-8 20 19.8 3.7830 5.00E-7 4.8 0.0062 9.95E-8 26.3 0.0499 1.86E-7 5.9 0.0094 1.06E-7 30 28.7 5.0575 4.46E-7 5.6 0.0140 1.00E-7 43.2 0.1295 5.22E-8 9.3 0.0265 1.82E-7 40 38.6 3.0935 6.52E-7 7.1 0.0234 5.60E-8 54.1 0.2137 1.65E-7 11.9 0.0374 9.14E-8 50 42.7 1.9016 5.37E-7 5.3 0.0218 7.73E-8 61.5 0.3120 1.93E-8 10.4 0.0437 5.88E-8 60 52.1 2.5272 5.61E-7 6.6 0.0359 5.90E-8 78.7 0.4976 1.05E-8 13.9 0.0718 1.15E-7 70 60.2 3.7050 6.10E-7 9.9 0.0624 1.12E-7 94.4 0.7332 1.80E-7 18.7 0.1264 1.26E-7 80 58.0 4.1246 4.31E-7 8.9 0.0640 6.03E-8 98.5 0.8845 3.88E-8 17.5 0.1420 5.35E-8 90 78.2 11.170 6.28E-7 10.0 0.0905 2.23E-7 114.3 1.2745 1.46E-7 20.9 0.2028 1.46E-7 100 72.2 12.211 4.77E-7 7.5 0.0709 1.62E-7 110.8 1.6477 1.31E-7 16.9 0.1881 1.34E-7
0 2 4 6 8 10 12 14
0.0 0.2 0.4 0.6 0.8 1.0
τ
ρs(τ) ϕ1(μ,t)
ϕ2(μ,t) ϕ3(μ,t) ϕ4(μ,t)
Figure 8: Performance profile of iteration numbers of Algorithm 3.1 for general AVE.
“iteration numbers”, we conclude that
φ2(µ, t) > φ4(µ, t) > φ1(µ, t) > φ3(µ, t) where “>” means “better performance”.
Finally, we are also interested in using the computing time as performance measure for Algorithm 3.1 with different φi(µ, t), i = 1, 2, 3, 4. The performance plot based on computing time is presented in Figure 9. From this figure, we can also see the function φ2(µ, t) has best performance, then followed by φ4(µ, t). Note that the time efficiency of φ1(µ, t) is very bad. Again, from the viewpoint of “computing time”, we conclude that
φ2(µ, t) > φ4(µ, t) > φ3(µ, t) > φ1(µ, t) where “>” means “better performance”.
0 5 10 15 20 25 30 0.0
0.2 0.4 0.6 0.8 1.0
τ
ρs(τ) ϕ1(μ,t)
ϕ2(μ,t) ϕ3(μ,t) ϕ4(μ,t)
Figure 9: Performance profile of computing time of Algorithm 3.1 for general AVE.
5 Conclusion
In this paper, we recast the AVE (1) as a system of smooth equations. Accordingly, we have proposed four smoothing functions along with a smoothing-type algorithm studied in [1, 7] to solve it. As mentioned in Section 2, there holds the local behavior shown as in (10):
φ3(µ, t) > φ1(µ, t) > φ2(µ, t) > |t| > φ4(µ, t).
and φ2(µ, t) is the one which best approximates the function |t| shown as in (9), i.e., φ3(µ, t) − |t|
∞>
φ1(µ, t) − |t|
∞>
φ4(µ, t) − |t|
∞>
φ2(µ, t) − |t|
∞. Surprisingly, φ2(µ, t) is also the best choice of smoothing function no matter when the iteration number or the computing time is taken into account. For the “iteration” aspect, the order of numerical performance from good to bad is
φ2(µ, t) > φ4(µ, t) > φ1(µ, t) > φ3(µ, t), for th AVE (1).
φ2(µ, t) > φ4(µ, t) > φ3(µ, t) > φ1(µ, t), for th AVE (2).
whereas for the “time” aspect, the order of numerical performance from good to bad is
φ2(µ, t) > φ4(µ, t) > φ3(µ, t) > φ1(µ, t), for th AVE (1).
φ2(µ, t) > φ3(µ, t) > φ4(µ, t) > φ1(µ, t), for th AVE (2).
In other words, φ2(µ, t) is the best choice of smoothing function to work with the proposed smoothing-type algorithm, meanwhile it also best approximate the function |t|. This is a very interesting discovery which may be helpful in other contexts. One of future directions is to check whether such phenomenon occurs in other types of algorithms.
References
[1] L. Caccetta, B. Qu, and G.-L. Zhou, A globally and quadratically convergent method for absolute value equations, Computational Optimization and Applications, vol. 48, 45–58, 2011.
[2] J.-S. Chen, C.-H. Ko, Y.-D. Liu, and S.-P. Wang, New smoothing functions for solving a system of equalities and inequalities, Pacific Journal of Optimization, vol. 12, pp. 185–206, 2016.
[3] S.-L. Hu, Z.-H. Huang, and Q. Zhang, A generalized Newton method for absolute value equations associated with second order cones, Journal of Computational and Applied Mathematics, vol. 235, pp. 1490–1501, 2011.
[4] Z.-H. Huang, Locating a maximally complementary solution of the monotone NCP by using non-interior-point smoothing algorithms, Mathematical Methods of Opera- tion Research, vol 61, pp. 41–45, 2005.
[5] S.-L. Hu and Z-H. Huang, A note on absolute value equations, Optimization Letters, vol. 4, pp. 417–424, 2010.
[6] Z-H. Huang, Y. Zhang, and W. Wu, A smoothing-type algorithm for solving system of inequalities, Journal of Computational and Applied Mathematics, vol. 220, pp. 355–363, 2008.
[7] X. Jiang and Y. Zhang, A smoothing-type algorithm for absolute value equations, Journal of Industrial and Management Optimization, vol. 9, pp. 789–798, 2013.
[8] S. Ketabchi and H. Moosaei, Minimum norm solution to the absolute value equation in the convex case, Journal of Optimization Theory and Applications, vol.
154, pp. 1080–1087, 2012.
[9] O. L. Mangasarian, Absolute value programming, Computational Optimization and Applications, vol. 36, pp. 43–53, 2007.
[10] O. L. Mangasarian, Absolute value equation solution via concave minimization, Optimization Letters, vol. 1, pp. 3–5, 2007.
[11] O. L. Mangasarian, A generalized Newton method for absolute value equations, Optimization Letters, vol. 3, pp. 101–108, 2009.
[12] O. L. Mangasarian, Primal-dual bilinear programming solution of the absolute value equation, Optimization Letters, vol. 6, pp. 1527–1533, 2012.
[13] O. L. Mangasarian, Absolute value equation solution via dual complementarity, Optimization Letters, vol. 7, pp. 625–630, 2013.
[14] O. L. Mangasarian and R. R. Meyer, Absolute value equation, Linear Algebra and Its Applications, vol. 419, pp. 359–367, 2006.
[15] O. A. Prokopyev, On equivalent reformulations for absolute value equations, Computational Optimization and Applications, vol. 44, pp. 363–372, 2009.
[16] L. Qi, Convergence analysis of some algorithms for solving nonsmooth equations, Mathematics of Operations Research, vol. 18, pp. 227–244, 1993.
[17] L. Qi, D. Sun, and G.-L. Zhou, A new look at smoothing Newton methods for nonlinear complementarity problems and box constrained variational inequality problems, Mathematical Programming, vol. 87, pp. 1–35, 2000.
[18] J. Rohn, A theorem of the alternatives for the equation Ax + B|x| = b, Linear and Multilinear Algebra, vol. 52, 421–426, 2004.
[19] J. Rohn, Solvability of systems of interval linear equations and inequalities, in Lin- ear Optimization Problems with Inexact Data edited by. M. Fiedler, J. Nedoma, J.
Ramik, J. Rohn and K. Zimmermann, Springer, pp. 35–77, 2006.
[20] J. Rohn, An algorithm for solving the absolute value equation, Eletronic Journal of Linear Algebra, vol. 18, 589–599, 2009.
[21] S. Yamanaka and M. Fukushima, A brancd and bound method for the absolute value programs, Optimization, vol. 63, pp. 305–319, 2014.
[22] Y. Zhang and Z-H. Huang, A nonmonotone smoothing-type algorithm for solv- ing a system of equalities and inequalities, Journal of Computational and Applied Mathematics, vol. 233, pp. 2312–2321, 2010.
[23] C. Zhang and Q.-J. Wei, Global and finite convergence of a generalized Newton method for absolute value equations, Journal of Optimization Theory and Applica- tions, vol. 143, pp. 391–403, 2009.
[24] F.K. Haghani, On Generalized Traubs Method for Absolute Value Equations, Jour- nal of Optimization Theory and Applications, vol. 166, pp. 619–625, 2015.
[25] E. D. Dolan and J. J. More, Benchmarking optimization software with perfor- mance profiles, Mathematical Programming, vol. 91, pp. 201–213, 2002.
