In this paper, the study of fractional order partial differential equations is made by using the reliable algorithm of the new iterative method (NIM). The fractional derivatives are considered in the Caputo sense whose order belongs to the closed interval [0, 1]. The proposed method is directly extended to study the fractional-order Roseau-Hyman and fractional order inhomogeneous partial differential equations without any transformation to convert the given problem into integer order. The obtained results are compared with those obtained by Variational Iteration Method (VIM), Homotopy Perturbation Method (HPM), Laplace Variational Iteration Method (LVIM) and the Laplace Adominan Decomposition Method (LADM). The results obtained by NIM, show higher accuracy than HPM, LVIM and LADM. The accuracy of the proposed method improves by taking more iterations.
1.
Introduction
In this paper, we consider the following bound constrained nonlinear systems of equations:
where F(x)=(F1(x),F2(x),⋯,Fn(x))T, and Fi:Rn→R is a nonlinear continuously differentiable function whose gradient is available. We denote by F′(x)=(∇F1(x),∇F2(x),⋯,∇Fn(x))T the Jacobian matrix of F at a given point x. The set Ω⊆Rn is defined as
for some given lower and upper bounds satisfying −∞≤li<ui≤+∞ for all i=1,2,⋯,n.
The bound constrained nonlinear equation of the type (1.1) is an important problem in the practical problems. There are a couple of different mathematical programming problems like Karush-Kuhn-Tucker systems and complementarity problems can be reformulated as the problem (1.1), see [1,2,3,4,5,6,7,8,9]. On the other hand, in many cases, the function Fi(x) is not always defined on the whole space Rn, and one usually puts some suitable bounds on some or all of the variables.
The Newton type method is one of the most important numerical methods for problem (1.1) and many researchers are interested in this method [7,10,11,12,13,14]. Given a current iterate xk∈Ω, the Newton method considers the least-squares solutions dk of the following nonlinear constrained equation:
We set the next iterate to be xk+1=xk+dk and call (1.2) the constrained Gauss-Newton method.
Another natural possibility is to consider solving the basic unconstrained Newton equation:
Denote the solution of (1.3) by dkN if it exists and then define xk+1 as the projection of xk+dkN onto Ω. This scheme can be called the projected Newton method.
On the other hand, there are many versions of the Newton type method, such as the constrained Levenberg-Marquardt method [6,15,16] usually used to solve the following subproblems:
where σ is a positive constant.
Along with constrained versions of the methods in question, one can also consider their projected variants. The projected Levenberg-Marquardt method has been proposed in [15] and its iteration consists of finding the solution dkLM of the unconstrained subproblem
and then defines the next iterate xk+1 as the projection of xk+dkLM onto Ω.
The Newton iteration can be costly, since partial derivatives must be computed and the linear system (1.3) must be solved at every iteration. This fact motivates the development of quasi-Newton methods [10,14,17] which are defined as the generalizations of (1.3) given by
In quasi-Newton methods, the matrices Bk are intended to be approximations of F′(xk) and be updated by some quasi-Newton formulas. Another well known algorithm is the trust region type algorithm, for example, [3,4,9,18,19,20,21].
Whether Newton method or quasi-Newton, one has to solve a linear system with full dimension, which will be expensive for large scale problems. To overcome this drawback, the active set methods are developed by many authors [7,12,13,22]. Since only a reduced dimension linear system to be dealt with at each iteration, the active set Newton methods are more efficient than the full Newton method especially for large scale problems.
To prove global convergence of the method outlined above, one often assumes that the iteration sequence is contained in a bounded set. If li and ui are bounded and the algorithm generates a feasible sequence, the assumption holds naturally. Otherwise, one often makes an assumption that the level set is bounded. For unconstrained nonlinear equations system, M.Solodov designed a Newton method with projection technique, the method can generate a bounded iteration sequence without additional assumption and the global convergence is obtained. Motivated by the idea of M.Solodov [23], in this paper, we extend the method to constrained equations (1.1). By using this active set strategy, we only need to solve a linear system with reduced dimension at each iteration. The algorithm generates a bounded sequence automatically even if li and ui are infinite. We obtain the global convergence and give numerical tests to show the efficiency of the proposed algorithm.
The paper is organized as follows: In section 2, we describe our algorithm in detail. In section 3, we prove the global convergence of the proposed algorithm. Some numerical tests are shown in Section 4 and a conclusion is given in section 5. Throughout this paper, we use ‖⋅‖ to denote the 2−norm and E denotes the identity matrix.
2.
Algorithm
We now describe our active set quasi-Newton method with projection technique in detail. To describe our algorithm, we introduce the definition of projection operator which is defined as a mapping form Rn to a nonempty closed convex subset Ω:
A well-known property of the operator is that it is nonexpensive, namely,
Given a current iterate xk, let
where δ and c are positive constants such that
and define the index sets
The precise statement of our algorithm is as follows:
Algorithm 2.1: (Active Set-type Quasi-Newton Method)
(S.0) Choose a positive definite matrix Bk, x0∈[l,u], choose parameters β∈(0,1), λ∈(0,1), δ>0, c>0, ε>0, μk>0, and ρk∈[0,1), and set k:=0.
(S.1) If ‖F(xk)‖≤ε, stop.
(S.2) Try to compute a vector dk∈Rn in the following way:
For i∈Ak, set
For i∈Ik, set solve the linear system
where
(S.3) Find zk=xk+αkdk, where αk=βmk with mk being the smallest nonnegative integer m such that
(S.4) Compute
(S.5) Update Bk+1, set k:=k+1, go to (S.1).
Just as mentioned in [13], throughout this paper, we assume that the parameter δ>0 is chosen sufficiently small such that
This implies that we cannot have xki−li≤δk and ui−xki≤δk for the same index i∈Ak.
Our algorithm is somewhat different from the traditional active set Newton method as described in [13], where the search step dk in (S.2) is computed in the following formulas:
For i∈Ak, set
For i∈Ik, solve the linear system
As described in [13], in order to understand the formula for the computation of the components dki for i∈Ik, note that, after a possible permutation of the rows and columns, [13] rewrite the standard (unconstrained) Newton equation F′(xk)d=−F(xk) as
Here we replace (2.7) by (2.3) and (2.8) by (2.4), the main proposal is to guarantee that the inequality (2.5) holds. On the other hand, we compute dIk by (2.4) instead of (2.8) which can be seen as an inexact Newton method.
The matrix Bk is updated by the well known rank two secant type formula updated by the well known BFGS formula
where yk=F(xk+1)−F(xk) and sk=xk+1−xk.
3.
Convergence
In the section, we prove the global convergence of Algorithm 2.1, we make the following assumption.
Assumption:
(A1) The function F(x) is Lipschitz continuous and monotone, i.e., there exists a positive constant L such that
and
(A2) The sequence of matrices {Bk} is positive definite and bounded, i.e., there exists a positive constant κ such that ‖Bk‖≤κ for all k.
We first show that the algorithm is feasible, i.e., there exists a positive m such that (2.5) holds.
Lemma 3.1. The Algorithm 2.1 is well defined.
Proof. We prove that the inequality (2.5) will hold with a nonnegative integer m. Suppose that for some index k this is not the case, which means, for all integer m, we have
We further get
Now we take the limit of both sides of (3.4) as m→∞, when (3.4) holds which implies that λ≥1, which contradicts the choice of λ∈(0,1). Hence we have that the inequality holds for some integer m, and the whole algorithm is well defined.
In what follows, we assume that the algorithm generates an infinite iteration sequence. The following result shows that the algorithm generates a bounded sequence automatically and the proof is similar to Lemma 3.2 in [24] and we omit it here.
Theorem 3.2. Suppose assumptions (A1) and (A2) hold, sequences {xk} and {zk} are generated by Algorithm 2.1, then {xk} and {zk} are both bounded. Furthermore, for any ¯x such that F(¯x)=0, it holds that
and
Now we give the global convergence result of the Algorithm 2.1.
Lemma 3.3. Let {xk} be generated by Algorithm 2.1, assume Assumption (A1) and (A2) hold, and there exists constants 0<ρ_<¯ρ<1, and μ_<¯μ such that ρ_≤ρk≤¯ρ, and μ_≤μk≤¯μ. Then {xk} converges to some x∗ such that F(x∗)=0.
Proof. By the inequality (2.5), we have
By the definition of dk, we have that
and
Combining (3.9) and (3.10), we can assume that there exists a positive constant c1 such that
On the other hand, the definition of dk also gives that
From (2.4) and Assumption (A2), we have
Combining (3.12) and (3.13), we can assume that there exists a positive constant c2 such that
Now by (3.8), we obtain
By the continuity of F(x), the bound of sequence {zk} and (3.6), we have
We consider the two possible cases:
In the first case, the continuity of F and the boundness of {xk} imply that the sequence {xk} has some accumulation point x∗ such that F(x∗)=0. Since ¯x was an arbitrary solution, we can choose ¯x=x∗ in (3.5). The sequence {‖xk−x∗‖} converges and since x∗ is an accumulation point of {xk}, it must be the case that {xk} converges to x∗.
Now consider the second case. From (3.14), we have
Hence by (3.16), we have
(The following proof is very similar to the last part in Theorem 2.1 [23], for complement, we list it here.) By the step rule, we have the inequality (2.5) is not valid for the value βmk−1, i.e.,
Let k→∞, we get
Here x∗, d∗, ρ∗, μ∗ denote the limits of the corresponding sequence respectively. On the other hand, by (3.4), we get
that contradicts the choice for λ∈(0,1). Hence the case lim infk→∞‖F(xk)‖ is impossible.
This completes the proof.
4.
Numerical experiments
In this section, we demonstrate the numerical performance of Algorithm 2.1 (AQN) and its computational advantage by comparing with the modified Kanzow [13] ACTN method (denoted as AKP) and the classical Quasi-Newton method with project (denoted as CQN). All presented codes are written in MATLAB2019 and run on a PC with 3.30GHz CPU processor, 4.0GB memory and Windows 8 operation system.
We consider ten problems with dimension n = 1000, 5000, 10000. We use six different starting points, that is:
After several parameter selection experiments, we select the initial parameters that can make the three algorithms have better performance :
Set the terminating criterion for the iteration process as ||F(xk)||≤10−6. The problems are listed as follows.
Problem 1. [25]
where Ω=Rn+.
Problem 2. [25]
where Ω=Rn+.
Problem 3. [25]
where Ω=Rn+.
Problem 4. [25]
where Ω=Rn+.
Problem 5. [25]
where Ω=Rn+.
Problem 6. [25]
where h=1n+1 and Ω=Rn+.
Problem 7. [25]
where Ω=Rn+.
Problem 8. [26]
where Ω=Rn+.
Problem 9. [26]
where Ω=Rn+.
Problem 10. [24]
where Ω=Rn+.
Comprehensive results of our numerical experiment are presented in Tables 1–10. The columns of the presented tables have the following definitions:
IP: the initial points.
DIM: the dimension of the problem.
NI: the iterative number.
NF: the iterative number of function evaluation.
CPU: the CPU time in seconds when the algorithm terminate.
NORM: the final norm equation.
We denote result by '−' whenever the number of iterations exceeds 500 or the terminating criterion has not been satisfied. Among these results, none of the three methods were able to solve Problem 9 when initial point is x3=(2,2,...,2)T. Therefore, Table 9 does not include the case when the initial point is x3. Meanwhile, in the drawing process, when the result was denoted by '−', its NI, NF, CPU and NORM are counted as ∞.
The performance of the three methods was evaluated using the performance profile which is presented by Dolan and Moré [27]. We comparing three methods with the same problem, dimension and initial point in an experiment, and recoding information of interest such as NI, NF, CPU and NORM.
We denote the set of problems as P and the set of methods as M. For example, for each problem p and method m, we define
Compare the performance on problem p by method m with the best performance by any method on this problem, that is, we use the performance ratio
We assume that a parameter R≥rp,m for all p,m is chosen, and rp,m=R if and only if method m does not solve problem p. If method m can solve problem p successfully, we obtain an overall assessment of the performance between these methods. It can be described as follows:
where np represents the number of elements in set P, then ρm(τ) is the probability for method m∈M that a performance ratio rp,m is within a factor τ∈R of the best possible ratio. The function ρm is the (cumulative) distribution function for the performance ratio.
The performance profile ρm:R⟼[0,1] for a method is a nondecreasing, piecewise constant function, continuous from the right at each breakpoint. We are interested in methods with a high probability of solve success, then we need only to compare the values of ρm(τ) for all of the methods and choose the method with the largest, there means that we need to find which method's function ρm first rearch the line ρm(τ)=1. In the same way, we can obtain the performance profile with respect to NI, NF and NORM.
As can be seen from the information in the Table 1–10, AQN has more stable solving performance and can solve more problems, such as what AKP cannot solve: x3 and x5 of Problem 3 when n=10000; x2 of Problem 3 when n=5000,10000; x2 of Problem 5 when n=1000,5000,10000; x3 of Problem 5 when n=5000,10000; x2 of Problem 8 when n=5000,10000; x2 of Problem 10 when n=5000,10000. Compared with CQN, from Figure 1 and Figure 2, we can see that AQN reaches the line that ρm(1)=1 before CQN, which demonstrates AQN has a faster solution time(CPU) and the solution results of the final norm equation(NORM) are more accurate. Although from Figure 3 and Figure 4, there shows the iterative number of CQN is less than AQN, in practice, we pay more attention to the advantage of solution time. To sum up, AQN has a more stable and faster solving performance.
5.
Conclusions
In this paper, we propose an active set quasi-Newton method for the solution of optimization problem with bound constraints. The implementation of the method uses the quasi-Newton step as a trial step and the project step as the correction step. By using active set technique, we only need to solve a reduced dimension linear equation at each iteration to generate the search direction. We prove that the generated sequence is bounded automatically and obtain the global convergence of the proposed algorithm. Meanwhile, compared with other algorithms, our method has the most stable performance. There are some questions that need studying in the near future. Firstly, it is possible to get the global convergence of the proposed algorithm without the assumption of the positive definite of the matrix Bk. Secondly, how to get the local convergence of the proposed algorithm especially under some weak condition such as the local error bound condition needs further studying.
Conflict of interest
All authors declare no conflicts of interest in this paper.