# A Derivative-Free Affine Scaling LP Algorithm Based on Probabilistic Models for Linear Inequality Constrained Optimization

Authors: Jingbo Zhao; Peng Wang; Detong Zhu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, a derivative-free affine scaling linear programming algorithm based on probabilistic models is considered for solving linear inequality constrainted optimization problems. The proposed algorithm is designed to build probabilistic linear polynomial interpolation models using onlyn + 1 interpolation points for the objective function and reduce the computation cost for building interpolation function. We build the affine scaling linear programming methods which use probabilistic or random models and affine matrix within a classical linear programming framework. The backtracking line search technique can guarantee monotone descent of the objective function, and by using this technique, the new iterative points are located within the feasible region. Under some reasonable conditions, the global and local fast convergence of the algorithm is shown, and the results of numerical experiments are reported to show the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

### 1.1. Problem Description Motivation

In this paper, we consider the following nonlinear optimization problem with linear inequality constraints:(1)Pminfx,s.t.Ax≥b,where fx is the general nonlinear twice continuously differentiable functions, but none of their first-order or second-order derivatives is explicitly available. Ω=x|Ax>b is the feasible region where A=a1T,a2T,…,alTT∈Rl×n. Moreover, it is assumed that the interior of the feasible region intΩ=x|Ax>b is not empty.The focus of this paper is the analysis of a numerical scheme that utilizes affine scaling randomized models to minimize nonlinear functions with linear inequality constraints. The motivation comes from the algorithm for minimization of so-call black-box functions where values are computed. For such problems, function evaluations are costly and derivatives are typically unavailable and cannot be approximated. Hence, the derivative-free optimization is developed. Moreover, the list of applications including molecular geometry optimization, circuit design, groundwater community problems, medical image registration, dynamic pricing, and aircraft design [1]is diverse and growing. In recent years, some people developed derivative-free algorithms for solving unconstrained and constrained minimization problems. Kolda et al. [2] presented a new generating set search approach for minimizing functions subject to linear constraints. One of their main contributions is a new condition to define the set of conforming search directions that admits several computational advantages. Zhang and Conn [3] developed a framework for a class of derivative-free trust region algorithms for the least-squares minimization. These algorithms are designed to take advantage of the problem structure by building polynomial interpolation models for each function in the least-squares minimization. Gratton et al. [4] considered an implementation of a recursive model-based active-set trust-region method for solving bound-constrained nonlinear non-convex optimization problems without derivatives using the technique of self-correcting geometry proposed in [5]. Kolda et al. [2] analyzed the oracle complexity of first-order and derivative-free algorithms for smooth non-convex minimization. Billups et al. [6] proposed a derivative free algorithm for optimizing computationally expensive functions with computational error. The algorithm used the trust region regression method which used weighted regression to obtain more accurate model functions at each trust region iteration. Bueno et al. [7] introduced a new method for solving a constrained optimization problem in which the derivatives of the constraints were available but the derivatives of the objective function were not. The method is based on the inexact restoration framework, and the iteration was divided in two phases. The first phase considered improving the feasibility of the constraints, and the second phase minimized a suitable objective function subjective to a linear approximation of the constraints which must be solved using derivative-free methods. Wild and Shoemarker [8] introduced the global convergence of radial basis function trust region algorithms for derivative-free optimization. Their results extend the recent work of Conn et al. [9] to fully linear models that have a nonlinear term. Xue and Sun [10] proposed a derivative-free trust region algorithm for constrained minimization problems with separable structure, where derivatives of the objective function are not available and cannot be directly approximated. They constructed a quadratic interpolation model of the objective function around the current iterate and used the filter technique to ensure the feasibility and optimality of the iterative sequence.There is a variety of evidence supporting the claim that randomized models can yield both practical and theoretical benefits for deterministic optimization. As an example, the randomized coordinate descent method for large-scale convex deterministic optimization proposed in [11] yields better complexity results than cyclical coordinate descent. Most contemporary randomized methods generate random directions along which all that may be required is some minor level of descent in the objective f. Within direct search, the use of random positive spanning sets has also been recently investigated [12, 13] with gains in performance and convergence theory for nonsmooth problems. Bandeira et al. [14] considered the use of a probabilistic or random model within a classical trust region framework for the optimization of deterministic smooth general nonlinear functions.


### 1.2. Contribution and Structure of the Paper

In order to solve the linear inequality constrained optimization problem, the derivative-free affine scaling LP algorithm based on probabilistic models is introduced in this paper. Compared with the traditional derivative-free algorithm, the new derivative-free algorithm proposed in this paper has the following contributions:(1)
The randomized numerical algorithm is designed in this paper providing accurate enough approximations such that in each iteration, a sufficiently improving step is produced with higher probability than traditional random schemes(2)
The polynomial interpolation models are built only usingn + 1 interpolation points, which are far less than (n + 1) (n + 2)/2 interpolation points of building the quadratic interpolation function. The computational cost of building the interpolation function using n + 1 interpolation point is less than using (n + 1) (n + 2)/2 interpolation point.(3)
The backtracking line search technique is used such that the linear programming subproblem is solved only once for obtaining the search direction, and the total computational effort for completing one iteration is less than the traditional derivative-free algorithm. Meanwhile, this technique also ensures that each iteration point is a feasible interior point.This paper is organized as follows. In Section2, we introduce the first-order affine scaling linear programming method based on probabilistic models. The algorithm for solving linear inequality constrained optimization problems is given in Section 3. In Section 4, we give the global convergence of the algorithm. The fast local convergence of the algorithm is given in Section 5. Some numerical results are given in Section 6. Finally, we give some conclusions in Section 7.


## 1.1. Problem Description Motivation

In this paper, we consider the following nonlinear optimization problem with linear inequality constraints:(1)Pminfx,s.t.Ax≥b,where fx is the general nonlinear twice continuously differentiable functions, but none of their first-order or second-order derivatives is explicitly available. Ω=x|Ax>b is the feasible region where A=a1T,a2T,…,alTT∈Rl×n. Moreover, it is assumed that the interior of the feasible region intΩ=x|Ax>b is not empty.The focus of this paper is the analysis of a numerical scheme that utilizes affine scaling randomized models to minimize nonlinear functions with linear inequality constraints. The motivation comes from the algorithm for minimization of so-call black-box functions where values are computed. For such problems, function evaluations are costly and derivatives are typically unavailable and cannot be approximated. Hence, the derivative-free optimization is developed. Moreover, the list of applications including molecular geometry optimization, circuit design, groundwater community problems, medical image registration, dynamic pricing, and aircraft design [1]is diverse and growing. In recent years, some people developed derivative-free algorithms for solving unconstrained and constrained minimization problems. Kolda et al. [2] presented a new generating set search approach for minimizing functions subject to linear constraints. One of their main contributions is a new condition to define the set of conforming search directions that admits several computational advantages. Zhang and Conn [3] developed a framework for a class of derivative-free trust region algorithms for the least-squares minimization. These algorithms are designed to take advantage of the problem structure by building polynomial interpolation models for each function in the least-squares minimization. Gratton et al. [4] considered an implementation of a recursive model-based active-set trust-region method for solving bound-constrained nonlinear non-convex optimization problems without derivatives using the technique of self-correcting geometry proposed in [5]. Kolda et al. [2] analyzed the oracle complexity of first-order and derivative-free algorithms for smooth non-convex minimization. Billups et al. [6] proposed a derivative free algorithm for optimizing computationally expensive functions with computational error. The algorithm used the trust region regression method which used weighted regression to obtain more accurate model functions at each trust region iteration. Bueno et al. [7] introduced a new method for solving a constrained optimization problem in which the derivatives of the constraints were available but the derivatives of the objective function were not. The method is based on the inexact restoration framework, and the iteration was divided in two phases. The first phase considered improving the feasibility of the constraints, and the second phase minimized a suitable objective function subjective to a linear approximation of the constraints which must be solved using derivative-free methods. Wild and Shoemarker [8] introduced the global convergence of radial basis function trust region algorithms for derivative-free optimization. Their results extend the recent work of Conn et al. [9] to fully linear models that have a nonlinear term. Xue and Sun [10] proposed a derivative-free trust region algorithm for constrained minimization problems with separable structure, where derivatives of the objective function are not available and cannot be directly approximated. They constructed a quadratic interpolation model of the objective function around the current iterate and used the filter technique to ensure the feasibility and optimality of the iterative sequence.There is a variety of evidence supporting the claim that randomized models can yield both practical and theoretical benefits for deterministic optimization. As an example, the randomized coordinate descent method for large-scale convex deterministic optimization proposed in [11] yields better complexity results than cyclical coordinate descent. Most contemporary randomized methods generate random directions along which all that may be required is some minor level of descent in the objective f. Within direct search, the use of random positive spanning sets has also been recently investigated [12, 13] with gains in performance and convergence theory for nonsmooth problems. Bandeira et al. [14] considered the use of a probabilistic or random model within a classical trust region framework for the optimization of deterministic smooth general nonlinear functions.


## 1.2. Contribution and Structure of the Paper

In order to solve the linear inequality constrained optimization problem, the derivative-free affine scaling LP algorithm based on probabilistic models is introduced in this paper. Compared with the traditional derivative-free algorithm, the new derivative-free algorithm proposed in this paper has the following contributions:(1)
The randomized numerical algorithm is designed in this paper providing accurate enough approximations such that in each iteration, a sufficiently improving step is produced with higher probability than traditional random schemes(2)
The polynomial interpolation models are built only usingn + 1 interpolation points, which are far less than (n + 1) (n + 2)/2 interpolation points of building the quadratic interpolation function. The computational cost of building the interpolation function using n + 1 interpolation point is less than using (n + 1) (n + 2)/2 interpolation point.(3)
The backtracking line search technique is used such that the linear programming subproblem is solved only once for obtaining the search direction, and the total computational effort for completing one iteration is less than the traditional derivative-free algorithm. Meanwhile, this technique also ensures that each iteration point is a feasible interior point.This paper is organized as follows. In Section2, we introduce the first-order affine scaling linear programming method based on probabilistic models. The algorithm for solving linear inequality constrained optimization problems is given in Section 3. In Section 4, we give the global convergence of the algorithm. The fast local convergence of the algorithm is given in Section 5. Some numerical results are given in Section 6. Finally, we give some conclusions in Section 7.


## 2. First-Order Affine Scaling LP Method Based on Probabilistic Models

In this section, we first introduce the construction method of a stochastic fully linear model and then give the affine scaling linear programming method based on probabilistic models. In the classical trust-region method, there need (n + 1) (n + 2)/2 interpolation points to build a quadratic interpolation model mk of function f with in the ball Bxk,Δk where xk is the center point and Δk is the radius of the ball. It is obvious that the number of interpolation points is much larger than the dimension of problem (P) when the dimension increases. This brings great difficult to computer in storage and computation. In this section, we will introduce an affine linear programming model based on probabilistic models which need only n + 1 interpolation points. It means that we will build random fashion polynomial interpolation models. Firstly, we will give the models and some assumptions of the models.


### 2.1. Probabilistically Fully Linear Models

We propose the polynomialmx∈Pnd as proposed in [15], which is κef,κeg fully linear, and we consider the following linear models, written in the form:(2)mxk+s=mxk+gkTs,where gk=∇mxk. The following definition shows the accuracy of the model mxk.Definition 1 (see [14]).
We say that a functionm is a κeg,κef fully linear model of f on Bxk,Δk if, for every s∈B0,Δk,(3)∇fxk+s−mxk+s≤κegΔk,fxk+s−mxk+s≤κefΔk2.Definition 2.
(see [14]). We say that a sequence of random models Mk is probabilistically κeg,κef fully linear for a corresponding sequence Bxk,δk if the events
S^k = {Mk} is κeg,κef fully linear for a corresponding sequence Bxk,δk}satisfying the following submartingale-like condition:(4)PS^k|Fk−1M≥p,where Fk−1M=σM0,M1,…,Mk−1 is σ-algebra generated by M0,M1,…,Mk−1. Furthermore, if p≥1/2, then we say that the random models are probabilistically κeg,κef fully linear.
In Definition2, Xk and the interpolation radii δk are the random variables defined over the s-algebra generated by M0,M1,…,Mk−1. Mk depends on Xk and δk, and hence, Mk depends on M0,M1,…,Mk−1.


### 2.2. The Affine Scaling Linear Programming Method Based on Probabilistic Models

Using the scaling matrixDk in [15], the form is as follows:(5)Dk=diagAxk−b.According to the same convergence theory, for anyxk∈intΩ, the augmented affine scaling linear programming subproblem at the k-th iteration is defined as(6)Lkminϕ^kd=gkTd,s.t.−Δk≤Dk−1/2Adi≤Δk,i=1,2,…,m,−Δk≤di≤Δk,i=1,2,…,n,where gk=∇mxk.Same as reference [15], we have that a “good” decrease of the quadratic objective function in Lk along the gradient gk can lead to dual feasibility:(7)limk⟶∞gk−ATλk=0,andlimk⟶∞Dk1/2λk=0.


## 2.1. Probabilistically Fully Linear Models

We propose the polynomialmx∈Pnd as proposed in [15], which is κef,κeg fully linear, and we consider the following linear models, written in the form:(2)mxk+s=mxk+gkTs,where gk=∇mxk. The following definition shows the accuracy of the model mxk.Definition 1 (see [14]).
We say that a functionm is a κeg,κef fully linear model of f on Bxk,Δk if, for every s∈B0,Δk,(3)∇fxk+s−mxk+s≤κegΔk,fxk+s−mxk+s≤κefΔk2.Definition 2.
(see [14]). We say that a sequence of random models Mk is probabilistically κeg,κef fully linear for a corresponding sequence Bxk,δk if the events
S^k = {Mk} is κeg,κef fully linear for a corresponding sequence Bxk,δk}satisfying the following submartingale-like condition:(4)PS^k|Fk−1M≥p,where Fk−1M=σM0,M1,…,Mk−1 is σ-algebra generated by M0,M1,…,Mk−1. Furthermore, if p≥1/2, then we say that the random models are probabilistically κeg,κef fully linear.
In Definition2, Xk and the interpolation radii δk are the random variables defined over the s-algebra generated by M0,M1,…,Mk−1. Mk depends on Xk and δk, and hence, Mk depends on M0,M1,…,Mk−1.


## 2.2. The Affine Scaling Linear Programming Method Based on Probabilistic Models

Using the scaling matrixDk in [15], the form is as follows:(5)Dk=diagAxk−b.According to the same convergence theory, for anyxk∈intΩ, the augmented affine scaling linear programming subproblem at the k-th iteration is defined as(6)Lkminϕ^kd=gkTd,s.t.−Δk≤Dk−1/2Adi≤Δk,i=1,2,…,m,−Δk≤di≤Δk,i=1,2,…,n,where gk=∇mxk.Same as reference [15], we have that a “good” decrease of the quadratic objective function in Lk along the gradient gk can lead to dual feasibility:(7)limk⟶∞gk−ATλk=0,andlimk⟶∞Dk1/2λk=0.


## 3. The Algorithm

In Section2, we introduced the affine scaling LP method based on probabilistic models. In this section, we define the specific steps of the derivative-free affine scaling linear programming algorithm based on probabilistic models for solving problem (P).


### 3.1. Algorithm

#### 3.1.1. Initialization Step

Fix the positive parameterη1,η2,γ,Δmax with γ>1>η2>η1,β∈0,1/2, and ω∈0,1. Select initial k=0,Δ0≤Δmax, and x0∈intΩ.


#### 3.1.2. Main Step

(1)
Approximate the functionf in Bxk,Δk by the linear function m.(2)
Compute the affine scaling linear programming subproblem (Lk) and obtain the stepddk.(3)
Chooseαk=1,ω,ω2,…, until the following inequality is satisfied:(8)fxk+αkdk≤fxk+αkβgkTdkwithxk+αkdk∈Ω.(4)
Let(9)hk=αkdk,ifxk+αkdk∈intΩ,θkαkdk,otherwise,whereθk∈θl,1, for some 0<θl<1 and θk−1=Odk.(5)
Ifαk≥η1, then set xk+1=xk+hk and(10)Δk+1=γ−1Δk,ifgkTΓk1/2<η2Δk,minγΔk,Δmax,iffgkTΓk1/2≥η2Δk.(6)
Ifαk<η1, then set xk+1=xk and Δk+1=γ−1Δk. Increase k by one and go to step 1.Remark 1.
Each realization of the algorithm defines a sequence of realizations for the corresponding random variables, in particular,mk=Mkωk,xk=Xkωk, and Δk=δkωk.Remark 2.
Similarly, the definition ofαk in [15], for the search direction dk, we have(11)αk≤min−aiTxk−biaiTdk|−aiTxk−biaiTdk>0,i=1,2,…,l,
withαk=+∞ if −aiTxk−bi/aiTdk≤0 for all i. A key property of the scalar αk is that an arbitrary step αkdk to the point xk+αkdk does not violate any linear inequality constraints. To see this, we first observe that if −aiTxk−bi/aiTdk≤0 for some i=1,2,…,m, and hence aiTxk−bi>0 implies aiTdk>0. Therefore, for all αk∈0,+∞,(12)aiTxk+αkdk−bi=aiTxk−bi+αkaiTdk>0,which means that the i-th linear strict inequality constraint holds.
Ifαk=aiTxk−bi/aiTdk>0 for some i, hence aiTxk−bi>0 implies aiTdk<0. We have that from αk=aiTxk−bi/aiTdk,(13)aiTxk−bi≤−αkaiTdk.
Hence, (12) and (13) mean that no matter what cases the inequality aiTxk+αkdk>bi for all any i=1,2,…,l holds.


## 3.1. Algorithm

### 3.1.1. Initialization Step

Fix the positive parameterη1,η2,γ,Δmax with γ>1>η2>η1,β∈0,1/2, and ω∈0,1. Select initial k=0,Δ0≤Δmax, and x0∈intΩ.


### 3.1.2. Main Step

(1)
Approximate the functionf in Bxk,Δk by the linear function m.(2)
Compute the affine scaling linear programming subproblem (Lk) and obtain the stepddk.(3)
Chooseαk=1,ω,ω2,…, until the following inequality is satisfied:(8)fxk+αkdk≤fxk+αkβgkTdkwithxk+αkdk∈Ω.(4)
Let(9)hk=αkdk,ifxk+αkdk∈intΩ,θkαkdk,otherwise,whereθk∈θl,1, for some 0<θl<1 and θk−1=Odk.(5)
Ifαk≥η1, then set xk+1=xk+hk and(10)Δk+1=γ−1Δk,ifgkTΓk1/2<η2Δk,minγΔk,Δmax,iffgkTΓk1/2≥η2Δk.(6)
Ifαk<η1, then set xk+1=xk and Δk+1=γ−1Δk. Increase k by one and go to step 1.Remark 1.
Each realization of the algorithm defines a sequence of realizations for the corresponding random variables, in particular,mk=Mkωk,xk=Xkωk, and Δk=δkωk.Remark 2.
Similarly, the definition ofαk in [15], for the search direction dk, we have(11)αk≤min−aiTxk−biaiTdk|−aiTxk−biaiTdk>0,i=1,2,…,l,
withαk=+∞ if −aiTxk−bi/aiTdk≤0 for all i. A key property of the scalar αk is that an arbitrary step αkdk to the point xk+αkdk does not violate any linear inequality constraints. To see this, we first observe that if −aiTxk−bi/aiTdk≤0 for some i=1,2,…,m, and hence aiTxk−bi>0 implies aiTdk>0. Therefore, for all αk∈0,+∞,(12)aiTxk+αkdk−bi=aiTxk−bi+αkaiTdk>0,which means that the i-th linear strict inequality constraint holds.
Ifαk=aiTxk−bi/aiTdk>0 for some i, hence aiTxk−bi>0 implies aiTdk<0. We have that from αk=aiTxk−bi/aiTdk,(13)aiTxk−bi≤−αkaiTdk.
Hence, (12) and (13) mean that no matter what cases the inequality aiTxk+αkdk>bi for all any i=1,2,…,l holds.


## 3.1.1. Initialization Step

Fix the positive parameterη1,η2,γ,Δmax with γ>1>η2>η1,β∈0,1/2, and ω∈0,1. Select initial k=0,Δ0≤Δmax, and x0∈intΩ.


## 3.1.2. Main Step

(1)
Approximate the functionf in Bxk,Δk by the linear function m.(2)
Compute the affine scaling linear programming subproblem (Lk) and obtain the stepddk.(3)
Chooseαk=1,ω,ω2,…, until the following inequality is satisfied:(8)fxk+αkdk≤fxk+αkβgkTdkwithxk+αkdk∈Ω.(4)
Let(9)hk=αkdk,ifxk+αkdk∈intΩ,θkαkdk,otherwise,whereθk∈θl,1, for some 0<θl<1 and θk−1=Odk.(5)
Ifαk≥η1, then set xk+1=xk+hk and(10)Δk+1=γ−1Δk,ifgkTΓk1/2<η2Δk,minγΔk,Δmax,iffgkTΓk1/2≥η2Δk.(6)
Ifαk<η1, then set xk+1=xk and Δk+1=γ−1Δk. Increase k by one and go to step 1.Remark 1.
Each realization of the algorithm defines a sequence of realizations for the corresponding random variables, in particular,mk=Mkωk,xk=Xkωk, and Δk=δkωk.Remark 2.
Similarly, the definition ofαk in [15], for the search direction dk, we have(11)αk≤min−aiTxk−biaiTdk|−aiTxk−biaiTdk>0,i=1,2,…,l,
withαk=+∞ if −aiTxk−bi/aiTdk≤0 for all i. A key property of the scalar αk is that an arbitrary step αkdk to the point xk+αkdk does not violate any linear inequality constraints. To see this, we first observe that if −aiTxk−bi/aiTdk≤0 for some i=1,2,…,m, and hence aiTxk−bi>0 implies aiTdk>0. Therefore, for all αk∈0,+∞,(12)aiTxk+αkdk−bi=aiTxk−bi+αkaiTdk>0,which means that the i-th linear strict inequality constraint holds.
Ifαk=aiTxk−bi/aiTdk>0 for some i, hence aiTxk−bi>0 implies aiTdk<0. We have that from αk=aiTxk−bi/aiTdk,(13)aiTxk−bi≤−αkaiTdk.
Hence, (12) and (13) mean that no matter what cases the inequality aiTxk+αkdk>bi for all any i=1,2,…,l holds.


## 4. Global Convergence

In this section, we will introduce the global convergence of our algorithm, which means that the algorithm proposed in this paper converges to the first-order critical points. For the purpose of proving the global convergence of the algorithm, we assume that the initial iterate pointx0 is given. Then, all the subsequent iterates belong to the level set(14)Lx0=x∈Rn|fx≤fx0.However, the failed iterates may lie outside this set. In the setting considered in this paper, all potential iterates are restricted to the region(15)Lenlx0=Lx0∪x∈Lx0Bx,Δmax,=∪x∈Lx0Bx,Δmax,where Δmax is the upper bound on the size of the interpolation radius, as imposed by the algorithm, and we define conv (Lenlx0) to be the convex hull of Lenlx0.Next, we also give the following assumption for the global convergence of algorithm.Assumption 1.
Supposex0 and Δmax are given. We assume that f is continuously differentiable in an open set containing the set Lenlx0 and that ∇f is Lipschitz continuous on conv (Lenlx0) with constant κLf, i.e., for ∀x,y∈Lnelx0, we have that(16)∇fx−∇fy≤κLfx−y.
We also assume thatf is bounded from below on Lx0.Assumption 2.
Suppose the level setLx0 is bounded.
Based on solving the augmented linear programming subproblemLk, the following lemma shows that gkTdk<−gkTΓk1/2Δk holds.Lemma 1.
For every realization of the algorithm, letdk be the solution of subproblem Lk; then,(17)gkTdk<−gkTΓk1/2Δk.Proof.
Let(18)sk=ΓkgkTdk1/2Δk.
Since(19)Dk−1/2AΓkigkTΓk1/2Δk≤Dk−1/2AΓkgkTΓk1/2Δk.
BecausegkTΓk1/2=Γk2+Dk−1/2AΓk2, we have that(20)Dk−1/2AΓkigkTΓk1/2Δk≤Δk.
Hence, we have that−Δk≤Dk−1/2AΓki≤Δk. Because Γki/gkTΓk1/2Δk=Γki/Γk2+Dk−1/2AΓk2Δk≤Δk, we have that sk is a feasible solution of subproblem Lk.
Becausedk is the solution of subproblem Lk, we have that(21)φ^dk≤φ^sk.
The value ofφ^sk is(22)ϕ^sk=gkTsk=gkTΓkgkTΓk1/2Δk.
According to the definition ofΓk, we have that(23)gkTΓk=−gkTgk−ATμk=−gkTgk+ATAAT+Dk−1−ATμk=−gkTΓk,where μk=−AAT+Dk−1Agk.
By (21)–(23), we have that(24)ϕ^dk=gkTdk≤−gkTΓkΔk.
The following lemma shows that every iterative point of the algorithm is a feasible solution of problem (P), i.e.,xk+αkdk∈Ω.Lemma 2.
Letxk be a sequence generated by the algorithm; then, we have that there exists αk>0 such that xk+αkdk∈Ω.Proof.
According to the definition ofDk and −Δk≤Dk−1/2Adki≤Δk, we have that(25)−b−AxkiΔk≤Adki≤b−AxkiΔk.
BecauseAxk−bi=aiTxk−bi and Adki=aiTdk, we have that(26)−b−aiTxkiΔk≤aiTdk≤b−aiTxkiΔk.
IfaiTdk≤0, we have that aiTxk+dk−bi≤0 holds, i.e., αk=1>0. Otherwise, if aiTdk>0, then from (26), for i∈1,2,…,l, we have that(27)−aiTxk−biaiTdk≥1Δk≥1Δmax>0.
By (11) and (27), we have that there exists αk>0 such that Axk+αkdk≤b holds in the k-th iteration, i.e., there exists αk>0 such that xk+αkdk∈Ω.
The following lemma states that the interpolation radius converges to zero regardless of the realization of the model sequenceMk made by the algorithm.Lemma 3.
Suppose Assumptions1 and 2 hold, for every realization of the algorithm, we have that(28)limk⟶∞Δk=0.Proof.
First, we can see that ifαk<η1 holds as k⟶∞, then (28) holds according to the definition of Δ¯k and Step 6 of the algorithm. If αk≥η1 and gkTΓk1/2<η2Δk hold as k⟶∞, equation (28) holds according to the definition of Δ¯k and Step 5 of the algorithm. Otherwise, according to Step 3 of the algorithm, we have that(29)fxk−fxk+αkdk≥−αkβgkTdk.
By (17), we have that(30)fxk−fxk+αkdk≥αkβgkTΓk1/2Δk.
Assumption2 shows that fx is a bounded function; hence, fxk−fxk+αkdk⟶0 as k⟶∞. By (30), we know that αkβgkTΓk1/2Δk⟶0 as k⟶∞. If αk⟶0, then there exists a positive index k0 such that αk<η1 for all k>k0, we have that (28) holds. Otherwise, we have that αk>0 for all k, if limk⟶∞Δk=0, then the conclusion holds. Otherwise, we assume that Δk>0 for all k⟶∞, then (30) means that limk⟶∞gkTΓk1/2=0 holds; according to the modified method of Δk in Step 5, we have that Δk will decrease and tend to zero, which is a contradiction. Hence, we have that (28) holds.
The following lemma shows that in the presence of sufficient model accuracy, a successful step will be achieved, provided that the interpolation radius is sufficiently small relative to the size of the model gradient.Lemma 4.
Suppose Assumptions1 and 2 hold, if m is κeg,κef fully linear on Bxk,Δk and(31)Δk≤ω1−βgkTΓk1/2η1κLf+κeg,where κΓ is a positive constant such that Γk≤κΓ, then at the k-th iteration αk≥η1.Proof.
SinceA−Dk1/2 has full row rank in the compact set Lx0, μk is bounded, there exists a constant κμ>0 such that μk≤κμ. By the definition of Γk and Lemma 2, we have that(32)Γk=−gk−ATμk≤gk+Aμk≤κg+κμA≤κg+κμκA.
LetκΓ=κg+κμκA, we have that Γk≤κΓ.
By accepting the rule of the stepαkdk, we have that(33)fxk+αkωdk−fxk>βαkω∇mxkTdk.
Using the mean value theorem, we have the following equality:(34)fxk+αkωdk−fxk=αkω∇fxk+ξkαkωdkTdk=αkω∇fxk+ξkαkωdk−∇fxkTdk+αkω∇fxkTdk=αkω∇fxk+ξkαkωdk−∇fxkTdk+αkω∇fxk−∇mxkTdk+αkω∇mxkTdk.
withξk∈0,1. Since ∇fx is Lipschitz continuously differentiable with constant κLf, we have that(35)∇fxk+ξkαkωdk−∇fxkTdk≤κLfξkαkωdk2≤κLfαkωdk2≤κLfαkωΔk2.
Becausemxk is a κeg,κef fully linear model of f on Bxk,Δk and αk/ω<1, mxk is a κeg,κef fully linear model of f on Bxk,αk/ωΔk. According to Definition 1, we have that(36)∇fxk−∇mxkTdk≤αkωκegΔk2.
By (33), (34), (35), and (36), we have that(37)−1−βαkω∇mxkTdk<αkω2κLf+κegΔk2.
Dividingαk/ω in the two sides of (37) and noting 1−β>0, we have that(38)αk>−ω1−βκLf+κegΔk2gkTdk≥ω1−βgkTΓk1/2κLf+κegΔk.
By (33), we have that αk≥η1.
Next, we assume that the model used in the algorithm is probabilistically fully linear, and we will state an auxiliary result from the martingale literature that will be useful in our analysis.Theorem 1 (see [14]).
LetGk be a submartingale, i.e., a sequence of random variables which, for every k, are integrable EGk<∞ and EGk|Fk−1G≥Gk−1, where Fk−1G=σG0,…,Gk−1 is the σ-algebra generated by G0,…,Gk−1, and EGk|Fk−1G denotes the conditional expectation of Gk given the past history of events Fk−1G. We assume further that Gk−Gk−1≤M<∞ for every k. We consider the random events C = {limk⟶∞Gk exists and is finite} and D=limsupk⟶∞Gk=∞. Then, PC∪D=1.Proof.
The theorem is a simple extension of Theorem 5.3.1 in [16].
Now, we will give the global convergence of the algorithm.Theorem 2.
Suppose that the model sequenceMk is probabilistically κeg,κef fully linear for some positive constants κeg and κef. Let xk be a sequence of random iterates generated by the algorithm. Then, almost surely,(39)liminfk⟶∞∇fxkTΓ^k=0.Proof.
According to the definition ofS^k in Definition 2, we use the the following random walk in [14]:(40)Wk=∑i=0k21S^i−1,where 1S^i=1 if S^i occurs, 1S^i=0 otherwise. It is easily seen that Wk is a submartingale. The proof can be found in Theorem 2 [14]. By Theorem 4.6, we can also see that the event D=limsupk⟶∞Wk=∞ holds almost surely. Suppose there exist positive constant ε and index K1 such that(41)P∇fxkTΓ^k≥ε>ε.
For allk≥K1.
Letxk and Δk be any realization of Xk and δk, respectively. By Lemma 3, there exists K2 such that for all k≥K2,(42)Δk<b:=ω1−βεη1κLf+κeg,Δmaxγ.
From the definition of fully linear models, we have that(43)∇fxk−gk≤κegΔk⟶0ask⟶∞.
According to Assumption2, we have that AAT+Dk is a reversible matrix, and by AAT+Dkμk=−Agk and AAT+Dkμ^k=−A∇fxk, we have that(44)μk−μ^k=AAT+Dk−1Agk−∇fxk≤AAT+Dk−1Agk−∇fxk≤κ1κegΔk,where κ1 is an upper bound of AAT+Dk−1A.
From (43) and (44), we know that μk−μ^k⟶0. Hence, we have that(45)gkTΓk−∇fxkTΓ^k=gkTΓk−gkTΓ^k+gkTΓ^k−∇fxkTΓ^k≤gkΓk−Γ^k+Γ^kgk−∇fxk≤gkgk−∇fxk+Aμk−μ^k+Γ^kgk−∇fxk.
SinceADk1/2 has full row rank in the compact set Lx0, μk is bounded, there exists a constant κμ^ > 0 such that μ^k≤κμ^. Similarly, to prove (32), we have that there exists a positive constant κΓ^ such that Γ^≤κΓ^. Combining (28), (43), (44), and (45), we can obtain that(46)gkTΓk−∇fxkTΓ^k⟶0,as k⟶∞. By (41) and (46), we obtain that(47)gkTΓk−∇fxkTΓ^k≤ε2.
Hence, we have that(48)gkTΓk1/2≥ε2.
Using Lemma4 and (28), we obtain that αk≥η1 and(49)gkTΓk1/2≥ε2≥η2Δk.
By the construction of the algorithm, and the fact thatΔk≤Δmax/γ, we have that Δk+1=γΔk. Similar to the proof of Theorem 2 in [14], we have that(50)liminfk⟶∞∇fxkTΓ^k1/2=0,almost surely.
Next, we show that(51)limk⟶∞∇fxkTΓ^k1/2=0,almost surely. Before stating and proving the main theorem, we will give the following lemmas.Lemma 5 (see [14]).
LetZk be a sequence of nonnegative uniformly bounded random variables and Bk be a sequence of Bernoulli random variables such that(52)PBk=1|σB1,B2,…,Bk−1·σZ1,Z2,…,Zk−1≥12.LetP be the set of natural numbers k such that Bk=1 and N = N\P, where P and N are the random sequences. Then,(53)Prob∑i∈PZi<∞∩∑i∈NZi=∞=0.Proof.
The proof is similar to Lemma1 in [14].Lemma 6 (see [14]).
Letxk and Δk be sequence of random iterates and random trust region radii generated by the algorithm. We fix ε>0 and define the sequence Ki consisting of the natural numbers k for which ∇fxkTΓ^k>ε, where Ki is a sequence of random variables. Then,(54)∑k∈KiΔk<∞,almost surely.Proof.
The proof is similar to Lemma2 in [14].
Next, we will prove thatlimk⟶∞∇fxkTΓ^k=0 holds.Theorem 3.
Suppose Assumptions1 and 2 hold, we suppose that the model sequence Mk is probabilistically κeg,κef fully linear for some positive constants κeg and κef. Let Xk be a sequence of random iterates generated by the algorithm; then,(55)limk⟶∞∇fxkTΓ^k=0.Proof.
We suppose thatlimk⟶∞∇fxkTΓ^k=0 does not hold almost surely. Then, there exists ε>0 such that P∇fxkTΓ^k>2ε>ε holds for infinitely many k. Let Ki be a subsequence of the iterations for which ∇fxkTΓ^k>2ε. According to Theorem 3, we have that there exists a pair of integers k′,k″ such that 0<k′<k″ and ∇fXk′TΓ^k′≤ε, ∇fXk+1′TΓ^k+1′>ε and ∇fXk+1″TΓ^k+1″≤2ε, hence, for any k∈k′,k″, we have ε<∇fXkTΓ^k≤2ε. Now, we consider the sequence of these intervals kl′,kl″. For which ∇fxkTΓ^k>ε for k∈ki, where xk and Δk are the realizations of Xk and δk, respectively. Then, for any l, we have that(56)ε<∇fxk′TΓ^k′−∇fxk″TΓ^k″=∑j=k′k″−1∇fxjTΓ^j−∇fxj+1TΓ^j+1≤∑j=k′k″−1∇fxjTΓ^j−∇fxj+1TΓ^j+1≤∑j=k′k″−1∇fxjTΓ^j−∇fxjTΓ^j+1−∇fxjTΓ^j+1−∇fxj+1TΓ^j+1≤∑j=k′k″−1∇fxjΓ^j−Γ^j+1+Γ^j+1∇fxj−∇fxj+1≤∑j=k′k″−1κfΓ^j−Γ^j+1+κΓ^∇fxj−∇fxj+1.
Since∇fx is Lipschitz continuous, we have that(57)∇fxj−∇fxj+1≤κLfxj−xj+1≤αjdj≤αjnΔj.
According to the definition ofΓ^x, we have that(58)Γ^j−Γ^j+1=∇fxj−∇fxj+1−AAT+Dk−1A∇fxj−∇fxj+1≤∇fxj−∇fxj+1+AAT+Dk−1A∇fxj−∇fxj+1≤1+κ1∇fxj−∇fxj+1≤1+κ1αjnΔj.
Combining (56), (57), and (58), we have that(59)ε≤αjnκ2Δkl′+∑j=kl′kl″−1Δj,where κ2=κf1+κ1+κΓ^.
From (28), then, for any l large enough, Δkl′≤ε/4αjκDκLf. Similar to the proof of (38), we have that αj>0 for all index j∈ki; hence, ∑j=kl′+1kl″−1Δj>ε/2>0, which give us ∑j∈kiΔj=∞, which contradicts Lemma 6; hence, the conclusion is true.


## 5. Local Convergence

In this section, we will show the local convergence of the algorithm, i.e., we prove that our algorithm superlinearly converges to the first-order stable point of problem (P). LetXk be a sequence of random iterates generated by the algorithm and let xk be realizations of Xk. For the local convergence of the algorithm, we will give the following assumption.Next, we will discuss the convergence rate for the proposed algorithm. For this purpose, it first shows that for large enoughk, the step size αk≡1, i.e., hk=dk for large enough k.Lemma 7.
Suppose Assumptions1 and 2 hold, if the nondegenerate property of problem (P) holds at every limit point x∗ of xk. Moreover, if m is the fully linear model of f, for sufficiently large k, then the step αk≡1.Proof.
First, we suppose that the conclusion is not true, i.e., for everyk, the following inequality holds(60)0<fxk+dk−fxk−βgkTdk=∇fxkTdk+odk=∇fxk−gkTdk+1−βgkTdk+odk≤1−βgkTdk+odk+∇fxk−gkdk.
Because∇fxk−gk≤κegΔ¯k and limk⟶∞Δ¯k=0, then ∇fxk−gk⟶0 as k⟶∞. Hence, we divide dk on both sides of inequality in (60), and we have that(61)1−βgkTdkdk+odkdk>0.
By the preserving character of the limit and (61), we have that(62)limk⟶∞gkTdkdk>0,which is contrary with gkTdk≤0 in Lemma 1. Hence, the assumption is not true, and for sufficiently large k, then the step αk≡1.
The following theorem shows that the algorithm is local superlinear convergence to optimal solutionx∗.Theorem 4.
Suppose Assumptions1 and 2 hold, and let x∗ be the limit point of the sequence xk which is generated by the algorithm, if xk+1−x∗ATDk−1/2Axk−x∗≤0 and xk−x∗=OΔkτ with 0<τ<1, then(63)limk⟶∞xk+1−x∗xk−x∗=0,i.e., the algorithm is superlinear convergence.Proof.
According to the definition ofDk and the Step 3 of the algorithm, we have that xk is a strict interior point of problem (P). Hence, there exists a positive constant ρ such that(64)ρxk+1−x∗2≤xk+1−x∗TATDk−1/2Axk+1−x∗.
By (64), we have that(65)ρxk+1−x∗2≤xk+1−x∗TATDk−1/2Axk+1−x∗=xk+1−x∗TATDk−1/2Axk−x∗+dk=xk+1−x∗TATDk−1/2Axk−x∗+xk+1−x∗TATDk−1/2Adk≤xk+1−x∗TATDk−1/2Adk≤Dk−1/2AdkAxk+1−x∗≤nκAΔkxk+1−x∗.
By (65), we obtain that(66)xk+1−x∗≤nκAρΔk.
By (66), xk−x∗=OΔkτ and Δk1−τ⟶0 as k⟶∞, we have that(67)limk⟶∞xk+1−x∗xk−x∗=0,
i.e., the algorithm is superlinear convergence.


## 6. Numerical Examples

In this section, we present the numerical experiments of our algorithm. In order to show the effectiveness of our algorithm, we choose the same test problem and compare the calculation results with algorithms COBYLA [17] and SDPEN [18]. First, the parameters used in the algorithm are as follows:(68)γ=2,η1=0.1,η2=0.7,Δ0=5andΔmax=10.We test the algorithm for “HS” problems which come from [19]; the test problems are reported in Table 1. The program code is written in MATLAB and run in MATLAB 7.0 environment. In Table 1, n denotes the dimension of the problem, and Ineq. denotes the number of inequality constraints. In this section, we choose 46 test problems and compare the number of function evaluations used in our method with COBYLA [17] and SDPEN [18]. We report our results using performance profiles. Such profiles compare the number of function evaluations needed by each solver to achieve the desired accuracy in the objective function value. We use three different levels of accuracy: 2, 3, and 5 significant figures in fx∗.Table 1
Description of the test problems.
ProblemnIneq.ProblemnIneq.ProblemnIneq.HS1021HS3332HS9041HS1121HS3432HS9151HS1221HS3531HS9261HS1321HS3631HS9362HS1522HS3732HS9464HS1622HS4343HS9564HS1722HS5721HS9664HS1822HS6431HS9764HS1922HS6531HS9864HS2121HS6632HS10074HS2221HS7041HS10176HS2321HS7242HS10276HS2931HS7643HS10376HS3031HS8821HS10486HS3131HS8931HS10581--------------------------HS117155The interpolation sample set is important for the algorithm. In this paper, we generate the sample set by randomly selecting the sample points, i.e., we suppose thatxk is a current iterate point; we choose n random points in a ball of radius Δk around xk. The derivative-free algorithm uses the trust region method which constructs quadratic interpolation function using (n + 1) (n + 2)/2 sample point. According to Section 2.1, we know that the dimension of the interpolation basis matrix MΦ,Y is (n + 1) (n + 2)/2. For taking the quadratic interpolation function, we need to compute the inverse matrix of MΦ,Y; hence, the computation cost is very expensive, especially when the dimension of the problem is very large. Now, we use linear interpolation function, and the dimension of the interpolation basis matrix MΦ,Y is only (n + 1); the computation cost is also far less than the construction of the quadratic interpolation function. The comparison results of CPU time are recorded, as shown in Figure 1; from Figure 1, we can see that the CPU time of building quadratic interpolation function is far more than linear interpolation function.Figure 1
The CPU time of interpolating 100 times.Using the performance profile of Dolan and Moré [20], we can demonstrate the overall behavior of the present algorithm and get more insights about their performance. Let S be the set of all algorithms and P be the set of test problems, with ns=S and np=P. For each problem p and solver s, let τp,s denote the quantities we want to compare for problem p and solver s. The performance ratio is defined as(69)rp,s=τp,sminτp,u:u∈S.In order to compare the performance of solvers on problem p with the best performance by any solver S on P. For any factor t≥1, the overall performance of algorithm s is given by(70)ρst=1npsizep∈P:rp,s≤t.The values on the left side of the figures represent the efficiency of each solver, and the values on the right side give information about the robustness of the solvers. This means that the best solver is the highest on the figures.Figure2 shows the comparison results of using our algorithm, COBYLA, and SDPEN to solve the test problems in Table 1 under termination accuracy level 2 (ε=10−2), respectively. It is obvious from the figure that the speed of the red curve (representing our algorithm) tending to 1 is faster than the other two curves, and when t=1, the value of the red curve is obviously higher than the other two curves. This shows that when the termination level is 2, the number of iterations spent by our algorithm is less than the other two algorithms when we solve the test problems in Table 1, and our algorithm has high computational efficiency under the condition of the low termination level.Figure 2
Performance proﬁle for the number of iterates in level 2.Figure3 records the results of the number of iterations spent by the three algorithms in solving the test problem in Table 1 when the termination level is 3 (ε=10−3). It can be clearly seen from Figure 2 that when the termination level is increased to 3, the performance rate of the three algorithms is lower than that of accuracy level 2. However, when t=1, the value of the red curve is still significantly higher than the other two curves, and the speed of red curve tending to 1 is significantly faster than the other two curves. This shows that after the accuracy level is improved, the number of iterations spent by our algorithm in solving the test problem is still less than that of the other two algorithms, and its computational efficiency is still high.Figure 3
Performance proﬁle for the number of iterates in level 3.Figure4 records the calculation results of three algorithms when the termination level is 5 (ε=10−5). It can be seen from Figure 3 that although the computational efficiency of the three algorithms has decreased significantly at a high accuracy level, the red curve still shows a high level. Whether from the value when t=1 or the speed tending to 1, the red curve shows good characteristics, that is, under the condition of high accuracy level, the number of iterations spent by our algorithm in solving the test problem in Table 1 is still less than that of the other two algorithms.Figure 4
Performance proﬁle for the number of iterates in level 5.Figure1 reports the results of the CPU time of constructing interpolation functions in 100 time iteration using (n + 1) and (n + 1) (n + 2)/2 interpolation points, respectively. From the result, we can see that the interpolation function can be established faster by using our algorithm. Lemma 1 shows that the linear interpolation model has the same drop as the two interpolation trust region model. Hence, our algorithm is more effective than the algorithm using two interpolation function.In order to test the effectiveness of our algorithm in solving problems with large dimensions, we selected 6 problems in Table2, calculated them with our algorithm, COBYLA, and SDPEN under different dimensions, and recorded the calculation results of level 5.Table 2
Description of the large-dimension test problems.
ProblemnBoundedProblemnBoundedRastrigin’s function50[−5.12, 5.12] 50Egg holder50[−512, 512] 50100[−5.12, 5.12] 100100[−512, 512] 100150[−5.12, 5.12] 150150[−512, 512] 150200[−5.12, 5.12] 200200[−512, 512] 200Schwefel’s function50[−512, 512] 50Rana’s function50[−500, 500] 50100[−512, 512] 100100[−500, 500] 100150[−512, 512] 150150[−500, 500] 150200[−512, 512] 200200[−500, 500] 200Stretch V sine wave function50[−10, 10] 50Keane’s function50[0, 10] 50100[−10, 10] 100100[0, 10] 100150[−10, 10] 150150[0, 10] 150200[−10, 10] 200200[0, 10] 200When the accuracy level is 5 (the termination accuracy isε=10−5), we use our algorithm to calculate the test problems in Table 2. At the same time, we use COBYLA and SDPEN to calculate the test problems in Table 2, compare the number of iterations spent by the three algorithms in calculating the test problems, and record the comparison results in Figure 5.Figure 5
Performance proﬁle for the number of iterates.From Figure5, we can see that the red curve representing our algorithm tending to 1 is the fastest, that is, the red curve tends to 1 faster than the curves representing the other two algorithms. Moreover, we can also see that when t=1, the computational efficiency of our algorithm has reached 0.38, while COBYLA rectification has reached 0.25 and SDPEN has only reached 0.11, which shows that our algorithm not only spends less the number of iterations than COBYLA and SDPEN in solving low-dimensional problems but also shows high computational efficiency. Moreover, the number of iterations is still lower than that of COBYLA and SDPEN. In short, our algorithm has high computational efficiency when calculating the actual test problems.


## 7. Conclusion

This paper has described the derivative-free linear programming methods based on probabilistic models for linear inequality constrained optimization. The algorithm is based on methods built by linear polynomial interpolations, which use onlyn + 1 interpolation point and far less than (n + 1) (n + 2)/2 interpolation point in quadratic polynomial interpolations and reduce the computation cost of building interpolation function. The affine scaling linear methods which use probabilistic or random models are built to obtain the descent direction. In this paper, we use the backtracking line search technique to ensure that the step length along the direction of descent is in the feasible region, and the objective function is reduced. Under some reasonable conditions, the global and local fast convergence of the algorithm is shown and the results of numerical experiments are reported to show the effectiveness of the proposed algorithm [21–26].



---

*Source: 5474937-2022-03-27.xml*