The Application of Fuzzication for Solving Quadratic Functions

This paper describes an approximation algorithm for solving standard quadratic optimization problems(StQPs) over the standard simplex by using fuzziﬁcation technique. We show that the approximate solution of the algorithm is an ǫ -critical point and satisﬁes ǫ - δ condition. The algorithm is compared with IBM ILOG CPLEX (short for CPLEX). The computational results indicate that the new algorithm is faster than CPLEX. Especially for infeasible problems. Furthermore, we calculate 100 instances for diﬀerent size StQP problems. The numerical experiments show that the average computational time of the new algorithm for calculating the ﬁrst local minimizer is in (cid:13) ( n ) when the size of the problems is less or equal to 450.


The problem
A standard quadratic optimization problem (StQP) is described as follows. (1.1) where Q is an arbitrary symmetric n × n matrix and e is the n-vector of all ones. The research of approximate algorithms is important for solving NP-hard problems. A StQP is a NP-hard problem when the matrix is indefinite (Bomze , 1998;Nowak , 1999). Ahmadi and Zhang show that it is NP-hard to decide if a quadratic function has a local minimizer over a polyhedron (Ahmadi and Zhang , 2020). Furthermore, they indicate that if there is a polynomial-time algorithm that finds a point within Euclidean distance c n (for any constant c ≥ 0) of a local minimizer of an n-variate quadratic function over a polytope, then P = N P (Ahmadi and Zhang , 2020).
On the other hand, researchers have proved that there are polynomial time approximate scheme (PTAS) for solving StQPs over the standard simplex (Bomze and de Klerk , 2002;De Klerk , 2008;De Klerk et al. , 2008).
Calculating local optimal solution is a major topic in nonconvex optimization problems (Agarwal et al. , 2016;Ahmadi and Zhang , 2020). Agarwal et al. present that finding approximate local minima for nonconvex optimization problems can be done in linear time (Agarwal et al. , 2016) and prove that the total time of computing local minima over nonconvex domains is ( nd ǫ 7/4 ), where d is the degree and n is the number of variables. It is clear that the computational time will become larger when ǫ becomes smaller such that if ǫ → 0 then time → ∞.
Fuzzy sets theory and possibility theory were invented half century ago (Zadeh , 1978). Fuzzy sets theory has been applied for solving optimization problems (Gao-1 , 2020;Gao-2 , 2020;Gao and Kawarada , 1995;Kaufmann and Gupta , 1998;Liu , 2007;Momot and Momot , 2009). A fuzzy number is a fuzzy set whose membership function is defined in R. Triangular fuzzy number(TFN) is a special case of a fuzzy number whose membership function forms a triangle (Dijkman et al. , 1983). The concept of fuzzification is commonly used in fuzzy sets theory and fuzzy control theory (Gao , 1999;Zimmermann , 2001). A fuzzification is the process of transforming system variables into fuzzy sets which are usually fuzzy numbers or TFNs. For multiple variables, fuzzification can map x ∈ R n to a set of TFNs. A set of TFNs is usually called a possibility distribution. In practice, normalized possibility distributions are used (Gao-1 , 2020). This paper is the continuation research of the previous paper (Gao-1 , 2020). We take advantage of piecewise linearity of TFNs and apply the concept of normalized possibility distribution to find optimal values for StQP problems. This paper addresses the following issues.
(1) The approximation solution of our algorithm is an ǫ-critical point.
(3) Our algorithm is compared with CPLEX. This paper is organized as follows. Section 2 gives a preliminary to basic concepts that are used in this paper. Section 3 describes the algorithm in detail. Section 4 provides the analysis of the algorithm and shows that the approximate solution satisfies ǫ − δ condition. Section 5 gives numerical experiment results. Section 6 gives the conclusion.

Preliminary
Definition 1 (ǫ−approximate local minimum (Agarwal et al. , 2016)) Suppose that function f (x)(x ∈ R n ) is twice differentiable. A point x * ∈ R n is called an ǫ−approximate local minimum if the following is satisfied.
where · denotes the Euclidean norm of a vector and I is the identity matrix. We say that a point x * is an ǫ−critical point if inequality (2.1) is satisfied.
In this paper, we only use the concept of ǫ-critical point such that the inequality (2.1).
where ˙ denotes the Euclidean norm of a vector.
Note that Hessian matrix plays a important role in optimization problems. The property of Hessian matrix determines if a solution is an optimal solution.
Definition 3 For a function f : R n → R and x ∈ R n , if f (x) is twice differentiable, then the Hessian matrix of f (x) at x = x 0 is defined as follows.
In this paper we will use the following results. If f (x 0 ) is a local minimum value, then the Hessian matrix is positive semidefinite.
If f (x 0 ) is a local maximum value, then the Hessian matrix is negative semidefinite.

The Algorithm
This section consists of three parts. The first part is to describe the fuzzification that transfers the space x ∈ R n into the space of a normalized possibility distribution. The second part is to prove that the approximate solution is an ǫ−critical point. The final part is to describe the transformed Hessian matrix at the ǫ−critical point after the fuzzification.

The Fuzzification
Fuzzification is the process of assigning input variables of a system to fuzzy sets by using membership functions. In this paper, the variable x ∈ R n is assigned to a set of normalized TFNs.
Suppose that λ i = (d l i , d m i , d r i ) is a TFN, and λ(t) = (λ 1 (t), · · · , λ n (t)) is a set of TFNs, that is considered as a possibility distribution. Since x ∈ △ in problem (1.1), it is required to normalize a given possibility distribution. A normalized TFN µ i (t) is defined as follows.
The fuzzification is to map x i to a normalized TFN µ i (t). As a result, variable x ∈ R n is mapped to a set of normalized TFNs. The mathematical form of the fuzzification is as follows.
A possibility distribution λ(t) is as follows.
The coefficients a i and b i are denoted as follows.
Since Q is a symmetric matrix, equation (3.6) becomes the following.
where K(t) and V are denoted in (3.8) and (3.9). (3.9) The necessary and sufficient condition of t * i being a solution of equation (3.7) is as follows (Gao-1 , 2020).
, respectively, Q is the symmetric matrix, and A and B are defined in (3.3).

µ(t * ) is an ǫ−critical points
In this section, we prove that if t * is a solution of (3.6), then µ(t * ) is an ǫcritical point of f (µ(t)). According to the definition 1 in section 2, for a given ǫ > 0, we need to prove the following. (3.11) Since we know that Let us calculate ∂f ∂µi (i = 1, 2, · · · , n). The first order partial derivative is as follows. ∂f Since x i and x j are independent, we assume that µ i (t) and µ j (t) are independent too. Therefore, ∂µ(t) ∂µi(t) is as follows: Since Q is a symmetric matrix, equation (3.12) becomes as follows: where Q T i = (q i1 , q i2 , · · · , q in ). If we use the notations K(t) and V in (3.8) and (3.9), then (3.13) becomes as follows.
We have the following theorem.
Theorem 1 Suppose that t * ∈ D i is a solution of (3.6). For a give ǫ > 0, If 14) then µ(t * ) is an ǫ−critical point.
Proof Based on the definition of ǫ−critical point, we need to prove the following.

Hessian matrix after the fuzzification
In this section, we will discuss how Hessian matrix is transformed into a different form of Hessian matrix after the fuzzification.
Hessian matrix Transformed Hessian matrices Fuzzification (3.2) Fig. 1: The transformation of Hessian matrix a solution of equation (3.7), then the transformed Hessian matrix is as follows.
where q ij is the element of matrix Q, V i is an element of V that is defined in (3.9), and C is a constant which is denoted as follows.
Proof According to equation (3.13), the second order partial derivative is as follows.
If we denote that C = 2 ( n j=1 λj (t)) 4 ) | t=t * and use the elements of V in (3.9), then the second order derivative becomes as follows.
Therefore, the transformed Hessian matrix in domain D i is as follows.
We know that µ(t * ) is an ǫ−critical point if t * is a solution of (3.6) and the inequalities in (3.14) is satisfied. It is necessary to investigate the transformed Hessian matrix after the fuzzification. This leads to the following theorem.
Theorem 3 If t * ∈ D i is a solution of equation (3.6) and S = V T QA > 0, then the transformed Hessian matrix is positive semi-definite in D i .
Proof According to the Theorem 3.2 in the paper (Gao-1 , 2020), when S = V T QA > 0, then f (µ(t * )) is a local minimum value such that f (µ(t) > f (µ(t * ) for all t ∈ D i . That is, the transformed Hessian matrix H ′ in D i is positive semi-definite.
Let us see an example. The following Hessian matrix H is indefinite because it has one negative eigenvalue and two positive eigenvalues.
If we define λ 3 (t) = 0 such that a 3 = 0 and b 3 = 0 and we establish a possibility distribution such as (λ 1 (t), λ 2 (t), 0). After the fuzzification, the Hessian matrix H becomes H ′ that is positive semi-definite because all of its eigenvalues are non-negative.
Theorem 3 provides a method to verify if a ǫ-critical point is an approximate minimum solution or not. The algorithm is depicted in Algorithm 1.

Algorithm 1
Given integer number M ≥ 1 while M > 0 do Step 1. Using random number generator generates vectors A and B to fzzify variable x ∈ R n with a normalized possibility distribution µ(t) ∈ R n . Determine sub domain D i = [d m i−1 , d m i ](i = 1, · · · , M ).
Step 2. Calculate V in D i based on (3.9).
Step 4. Verify if t * satisfies the inequalities (3.10). If t * / ∈ D i , go to next sub domain.
Step 5. Calculate S = V T QA, and decide f (µ(t * )) is a local minimum value or a local maximum value in D i .

Algorithm Analysis
In this section, we show that the approximate solution µ(t * ) satisfies δ − ǫ proof conditions when µ(t * ) is an ǫ−approximate local minimizer.
where ∇f (µ(t * )) is the gradient of f at µ(t * )), H(µ(t * )) is the transformed Hessian matrix at t * ∈ D i . We have the following.
The inequality in (4.3) becomes as follows.
Then, we have the following.
Let us denote δ as follows.

Numerical Results
In this section, first, we will show the comparison between the new algorithm and CPLEX that can be downloaded from IBM website(IBM , 2020). Secondly, we calculate 100 instances for each size from n=50 to n=450 to find out the first local minimizer, and the averages of computational times are measured. The numerical experiment results will be discussed. CPLEX is a well-known tool for computing optimal values of non-linear functions such as quadratic polynomial functions(IBM , 2020). The C++ methods of CPLEX are called with Microsoft Visual Studio 2015.

Comparison with CPLEX
In this section, numerical experiments are performed with a computer with Intel i5-7200U processor and 8 GB RAM is used.
When we use CPLEX, the parameters Param::MIP::Tolerances::AbsMIPGap and Param::MIP::Tolerances::MIPGap are set with 0.001 and 0.01. The parameter Param::MIP::Tolerances::AbsMIPGap indicates the absolute MIP gap tolerance. It sets an absolute tolerance on the gap between the best integer objective and the objective of the best node remaining. The default value of this parameter is 1e-06(IBM , 2020). The parameter Param::MIP::Tolerances::MIPGap indicates a relative MIP gap tolerance. It Sets a relative tolerance on the gap between the best integer objective and the objective of the best node remaining. The default value of this parameter is 1e-04 (IBM , 2020).
For CPLEX, The smaller the two parameters, the longer the computational time. If we use the default values for the two parameters, the computation of CPLEX will take longer times.
According to Table 1, when n = 10, for the 20 instances, CPLEX solution is same as the new algorithm solution for each instance. However, the new algorithm is faster than CPLEX. Based on the average time for the 20 instances, the new algorithm is ten times faster than CPLEX. Especially, for some instances, for example, instance number 7, 14, 18, and 20, CPLEX takes quite long time to calculate the solutions.
According to Table 2, when n = 20, for the 20 instances, compare with CPLEX, some of the new algorithm solutions have discrepancy. Especially for the instance number 5, 12, and 19, the new algorithm solutions are not global optimal solutions but local optimal solutions. If we ignore these discrepancies, the new algorithm is faster than CPLEX. For example, the instance number 2, 3, 6, 8, and 13. For some instances, for example, instance number 2, 8, 13, and 18, CPLEX takes quite long time to calculate the solutions.
The reason why CPLEX takes long time is because these instances are infeasible. A optimization problem is said to be infeasible if no solution which satisfies all the constraints exists(J. W. Chinneck , 2008). CPLEX is an effective algorithm for solving infeasible problems. The technique that CPLEX uses is called infeasibility isolation (IBM , 2020).  Compare with CPLEX, the infeasibility does not affect the new algorithm. The reason why the new algorithm is not affected by the infeasibility is that the new algorithm uses random number generators for the coefficients of the possibility distributions. A further research is required for clarifying the reason why infeasibility issue does not affect the new algorithm.

Numerical experiments
The numerical experiments in this section are performed with a computer with 6 core Intel i7-8700 CPUs and 16 GB RAM. The algorithm is implemented in C++ with Microsoft Visual Studio 2015.
First, we generate 100 instances of symmetric matrices for each size of 50, 100, 150, 200, 250, 300, 350, 400, and 450 StQPs by using random number generator. Secondly, we measure the computational time when the first minimizer of each instance has been found by running the new algorithm. Finally, we calculate the average computational time for each size. The average computation times for each size of StQPs are shown in Table 3. The data is visualized in Figure 1. Figure 1 indicates that the average computational time for calculating the first minimizer of StQP in in (n).  We propose a new algorithm for finding approximate solutions for StQP problems. This paper proves that (1) the approximate solution of the new algorithm is a ǫ-critical point and (2) The approximate solution satisfies ǫ-δ condition. Furthermore, we show that when the approximate solution is local minimum, then the transformed Hessian matrix is positive semi-definite. Table 1 and Table 2 indicate that the new algorithm has better performance than CPLEX in most cases. Note that CPLEX uses the Branch and Bound technique that takes more time when it deals with the infeasibility issue (J. W. Chinneck , 2008). On the other hand, our algorithm uses normalized possibility distributions that are generated with randomized algorithm technique which is capable to circumvent the infeasibility and achieve better performances.
Future research will involve the application of the new algorithm to solve practical problems such as portfolio optimization problems and the analysis of the new algorithm from the perspective of computational complexity theory.