RPCGB Method for Large-Scale Global Optimization Problems

Abstract

In this paper, we propose a new approach for optimizing a large-scale non-convex differentiable function subject to linear equality constraints. The proposed method, RPCGB (random perturbation of the conditional gradient method with bisection algorithm), computes a search direction by the conditional gradient, and an optimal line search is found by a bisection algorithm, which results in a decrease of the cost function. The RPCGB method is designed to guarantee global convergence of the algorithm. An implementation and testing of the method are given, with numerical results of large-scale problems that demonstrate its efficiency.

1. Introduction

Non-convex optimization is a type of mathematical optimization problem in which the objective function to be optimized is not convex. Unlike convex optimization problems, non-convex problems can have multiple local optima, which can make it difficult to find the global optimum. Non-convex optimization has many applications in various fields, including finance (portfolio optimization, risk management, and option pricing) [1,2,3], computer vision (image segmentation, object recognition) [4,5], signal processing (compressed sensing, channel estimation, and equalization) [6,7,8], engineering (control systems, optimization of structures) [9,10], machine learning [11,12,13,14], and damage characterization [15] based on deep neural networks and the YUKI algorithm [16].

To solve non-convex optimization problems, two broad classes of techniques have been developed: deterministic and stochastic methods [17]. Deterministic methods include gradient-based methods, which rely on computing gradients of the objective function, and which can be sensitive to the choice of initialization and can converge to local optima. On the other hand, stochastic methods use randomness to explore the search space and can be less sensitive to initialization and more likely to find the global optimum.

Several reasons make stochastic methods more appropriate for non-convex optimization problems. Stochastic methods avoid getting stuck in local optima or saddle points, as they explore the search space more thoroughly. Additionally, complex non-convex optimization problems often have a large number of variables, making gradient-based methods computationally expensive. In contrast, stochastic methods can scale better to high-dimensional problems. Stochastic methods can also be more robust to noise and uncertainty in the problem formulation.

Scientific studies have shown the effectiveness of stochastic methods in solving non-convex optimization problems. One of the most widely studied deterministic methods that has been extended with random perturbations is gradient descent. For example, a study by Pogu and Souza de Cursi (1994) [18] compared the performance of deterministic gradient descent and stochastic gradient descent on a variety of non-convex optimization problems and found that stochastic gradient descent was more robust and could converge to better solutions. Another study by Mandt et al. (2016) [19] investigated the use of random perturbations in the context of Bayesian optimization and found that it could lead to better exploration of the search space and improved optimization performance. In addition to stochastic gradient descent, other deterministic methods have also been extended with random perturbations. For example, a study by Nesterov and Spokoiny (2017) [20] proposed a variant of the conjugate gradient method that adds random noise to the search direction at each iteration and showed that it could improve the convergence rate and solution quality compared to the standard conjugate gradient method. Another study by Songtao Lu et al. (2019) [21] proposed a variant of the projected gradient descent method that added random perturbations to the method, and they demonstrated its effectiveness in solving non-convex optimization problems.

We consider non-convex optimization problems with linear equality or inequality constraints of the form

$$\left\{\begin{array}{cc}\hfill min\hfill & \hfill f\left(x\right)\hfill \\ \hfill s.t\hfill & \hfill Ax\le b\hfill \\ \hfill \hfill & \hfill \nu \le x\le \mu \hfill \end{array}\right.$$

where $f:{\mathbb{R}}^n\to \mathbb{R}$ is a continuously differentiable function, A is an $m\times n$ matrix with rank m, b is an m-vector, and the lower and upper bound vectors, $\nu$ and $\mu$, may contain some infinite components; and

$$\left\{\begin{array}{cc}\hfill min\hfill & \hfill f\left(x\right)\hfill \\ \hfill s.t\hfill & \hfill Ax=b\hfill \end{array}\right.$$

(1)

where $A\in {\mathbb{R}}^{m\times n}$, $b\in {\mathbb{R}}^m$, and $f:{\mathbb{R}}^n\to$$\mathbb{R}$ is a continuously differentiable, non-convex objective function.

One possible numerical method to solve problem (1) is the conditional gradient with bisection (CGB) method. This method generates a sequence of feasible points ${\left\{x^t\right\}}_{t\ge 0}$, starting with an initial feasible point $x^0$. A new feasible point $x^{t+1}$ is obtained from $x^t$ for each $t>0$, using an operator $Q_t$ (details can be found in Section 3). The iterations can be expressed as follows:

$$x^{t+1}=Q_t\left(x^t\right),\forall t\ge 0.$$

In this paper, we present a new approach for solving large-scale non-convex optimization problems by using a modified version of the conditional gradient algorithm that incorporates stochastic perturbations. The main contribution of this paper is to propose the RPCGB algorithm, which is an extension of a method previously presented in [22] that was designed for small- and medium-scale problems. The RPCGB algorithm was developed to deal with large-scale global optimization problems and aims to determine the global optimum.

This method involves replacing the sequence of vectors ${\left\{x^t\right\}}_{t\ge 0}$ with a sequence of random vectors ${\left\{{\mathbf{X}}^t\right\}}_{t\ge 0}$, and the iterations are modified as follows:

$${\mathbf{X}}^{t+1}=Q_t\left({\mathbf{X}}^t\right)+P_t,\forall t\ge 0,$$

where $P_t$ is a random variable that is chosen appropriately, which is commonly known as the stochastic perturbation. It is important that the sequence ${\left\{P_t\right\}}_{t\ge 0}$ converges to zero at a rate slow enough to avoid the sequence ${\left\{{\mathbf{X}}^t\right\}}_{t\ge 0}$ convergence to local minima. For more details, refer to Section 4.

The paper is structured as follows: Section 3 revisits the principle of the conditional gradient with bisection method, while Section 4 provides details on the random perturbation of the CGB method. Notations are introduced in Section 2, and in Section 5, the results of numerical experiments for non-convex optimization tests with linear constraints are presented for large-scale problems.

2. Notations and Assumptions

We denote by $\mathbb{R}$ the set of the real numbers, and $E={\mathbb{R}}^n$ is the n-dimensional real Euclidean space. $x^T$ denotes the transpose of x. We denote by $∥x∥=\sqrt{x^{T}x}={(x_1^2+\cdots +x_n^2)}^{1/2}$ the euclidean norm of x, and let

$$M=\{x\in E|Ax=b,x\ge 0\}$$

and ${\eta }^{*}={\displaystyle \underset{M}{min}}f$, the lower bound of f on M. Let us introduce

$$M_{\phi }=N_{\phi }\cap M;\mathrm{where}{N}_{\phi }=\{x\in E|f\left(x\right)\le \phi \}.$$

Supposing

$$\forall \phi >{\eta }^{*}:meas\left(M_{\phi }\right)>0,$$

(2)

$$\forall \phi >{\eta }^{*}:M_{\phi }\mathrm{is}\mathrm{not}\mathrm{empty},\mathrm{closed},\mathrm{and}\mathrm{bounded},$$

(3)

$$f\mathrm{is}\mathrm{continuously}\mathrm{differentiable}\mathrm{on}E,$$

(4)

where $meas\left(M_{\phi }\right)$ is the measure of $M_{\phi }$.

As the space E is of finite dimensions, condition (3) holds true if M is bounded or if f is coercive, i.e., ${\displaystyle \underset{∥x∥\to +\infty }{lim}}f\left(x\right)=+\infty$. Assumption (3) is satisfied when M comprises a series of neighborhoods of an optimal point $x^{*}$ that possesses a strictly positive measure, meaning $x^{*}$ can be approximated by a sequence of points from the interior of M.

We see that the results of assumptions (3) and (4) are

$$M=\bigcup _{\phi >{\eta }^{*}}{N}_{\phi },i.e.,\forall x\in M:\exists \phi >{\eta }^{*}\mathrm{such}\mathrm{that}x\in M_{\phi }.$$

From (3) and (4), one has:

$${\delta }_1=sup\left\{∥\nabla f\left(x\right)∥:x\in M_{\phi }\right\}<+\infty .$$

Consequently, one deduces

$${\delta }_2=sup\left\{∥d∥:x\in M_{\phi }\right\}<+\infty ,$$

where d is the direction of conditional gradient method.

Thus,

$$\rho \left(\phi ,\epsilon \right)=sup\left\{∥y-(x+\alpha d)∥:(x,y)\in M_{\phi }\times M_{\phi },0\le \alpha \le \epsilon \right\}<+\infty ,$$

(5)

where $\alpha$ and $\epsilon$ are positive real numbers.

3. Conditional Gradient Method

3.1. Conditional Gradient Algorithm

The conditional gradient method, also known as the Frank–Wolfe algorithm, is an iterative optimization algorithm used to find the minimum of a convex function over a convex set. It was introduced by Philip Wolfe and Marguerite Frank in 1956 [23] and is one of the oldest nonlinear constrained optimization techniques. It has recently gained renewed interest due to its projection-free iterations and low memory requirement. This algorithm enables the approximation of a function during each iteration by utilizing the first-order Taylor series expansion.

The algorithm starts with an initial point in the feasible set and iteratively moves towards a direction that minimizes the gradient of the objective function over the feasible set. At each iteration, the algorithm solves a linear optimization problem over the feasible set to find the direction that minimizes the gradient.

The conditional gradient algorithm has several advantages over other optimization methods, including its ability to handle large-scale problems and its ability to find sparse solutions. However, it may converge slowly and may not always find the global minimum. Moving forward, we focus on a problem related to nonlinear programming that involves constraints in the form of linear equalities or inequalities of the form

$$\left\{\begin{array}{c}\mathrm{minimize}f\left(x\right)\\ \\ \mathrm{subject}\mathrm{to}x\in M\end{array}\right.$$

(6)

The search direction is $d_t:=s_t-x^t$, with $s_t$ being the optimal solution of a linear programming problem and

$$x^{t+1}=Q_t\left(x^k\right)=x^t+{\alpha }_k{d}_t.$$

(7)

We ascertain the optimal step by selecting the value of ${\alpha }_t$ that satisfies

$$f\left(x^t+{\alpha }_t{d}_t\right)=\underset{0\le \alpha \le 1}{min}\left\{f\left(x^t+\alpha d_t\right)\right\}.$$

(8)

The conditional gradient algorithm can be summarized as follows (Algorithm 1):

Algorithm 1 Conditional gradient algorithm.

1 Choose an initial point $x^{\left(0\right)}\in M$ in the feasible set M.

2:$\mathbf{for}t=0,1,2,\dots ,T\mathbf{do}$

3: Compute ${\mathbf{s}}_t:=$LMO${}_M$$(\nabla f\left(x^{\left(t\right)}\right)):=\underset{s\in M}{argmin}$$\nabla f{\left(x^{\left(t\right)}\right)}^{\top }s$

($\mathbf{LMO}$: Linear minimization oracle)

4: Let ${\mathbf{d}}_t:={\mathbf{s}}_t-x^{\left(t\right)}$                (Conditional gradient direction)

5: Compute $g_t:=〈-\nabla f\left(x^{\left(t\right)}\right),{\mathbf{d}}_t〉$          (Conditional gradient gap)

6: $\mathbf{if}$ $g_t<\epsilon \mathbf{then}\mathbf{return}$ $x^{\left(t\right)}$

7: optimal line search step size

${\alpha }_k\in \underset{\alpha \in [0,1]}{argmin}f(x^{\left(k\right)}+\alpha {\mathbf{d}}_k)$

8: Update $x^{(t+1)}:=x^{\left(t\right)}+{\alpha }_k{\mathbf{d}}_t$

9: $\mathbf{end}\mathbf{for}$

10: $\mathbf{return}{x}^{\left(T\right)}$

For non-convex objectives, the conditional gradient algorithm may not converge to a global minimum, but it can still converge to a stationary point under certain conditions. Simon Lacoste-Julien and his colleagues have shown that the Frank–Wolfe algorithm can converge to a stationary point for non-convex objectives, as shown in [24].

3.2. Bisection Algorithm

In this paper, we employ the bisection algorithm to tackle the unconstrained optimization problem with one variable (8). The method is described in [25]. We refer to the recursive bisection procedure as bis $(h,{\theta }_1,{\theta }_2,ϵ)$, which takes as inputs the h calculation procedure, the $[{\theta }_1,{\theta }_2]$ interval, and the precision $ϵ$. The outputs of this procedure are an approximation of $x_m$ for the minimizer $x^{*}$ and $h_m$ for the minimum value of the h function over the $[{\theta }_1,{\theta }_2]$ interval.

The recursive procedure iteration involves the application of the subsequent steps.

Step 0: If ${\theta }_2-{\theta }_1\ge ϵ,$ go to step $1,$ otherwise stop.

Step 1: Compute

$${\theta }_3=\frac{{\theta }_1+{\theta }_2}{2},{\theta }_1^{\prime }=\frac{{\theta }_2+{\theta }_3}{2},{\theta }_2^{\prime }=\frac{{\theta }_3+{\theta }_2}{2},h\left({\theta }_3\right),h\left({\theta }_1^{\prime }\right),h\left({\theta }_2^{\prime }\right).$$

Step 2: If $h\left({\theta }_1^{\prime }\right)\le h\left({\theta }_3\right)\le h\left({\theta }_2^{\prime }\right),$ set ${\theta }_2={\theta }_2^{\prime }$ If $h\left({\theta }_1^{\prime }\right)\ge h\left({\theta }_3\right)\ge h\left({\theta }_2^{\prime }\right),$ set ${\theta }_1={\theta }_1^{\prime }$.

$$\mathrm{If}h\left({\theta }_3\right)\le min\left\{h\left({\theta }_1^{\prime }\right),h\left({\theta }_2^{\prime }\right)\right\},\mathrm{set}{\theta }_1={\theta }_1^{\prime },{\theta }_2={\theta }_2^{\prime }.$$

Step 3: Execute bis $(h,{\theta }_1,{\theta }_2,ϵ)$ with new inputs.

4. RPCGB Method

From [23], when it comes to objective functions that are non-convex, optimization algorithms based on gradients (CGB) cannot guarantee the discovery of the global minimum. Convex functions are the only ones for which CGB methods can find the global minimum. To deal with this issue, we suggest utilizing a suitable random perturbation method. Next, we will demonstrate how RPCGB can converge to a global minimum for non-convex optimization problems.

The sequence of real numbers ${\left\{x^t\right\}}_{t\ge 0}$ is replaced by a sequence of random variables ${\left\{{\mathbf{X}}^t\right\}}_{t\ge 0}$ involving a random perturbation $P_t$ of the deterministic iteration (7). We have ${\mathbf{X}}^0=$$x^0$;

$$\forall t\ge 0{\mathbf{X}}^{t+1}=Q_t\left({\mathbf{X}}^t\right)+P_t={\mathbf{X}}^t+{\alpha }_k{d}^t+P_t={\mathbf{X}}^t+{\alpha }_t(d^t+{\displaystyle \frac{P_t}{{\alpha }_t}}),$$

(9)

$$P_t\mathrm{is}\mathrm{independent}\mathrm{from}({\mathbf{X}}^{t-1},\dots ,{\mathbf{X}}^0),\forall t\ge 1,$$

where ${\alpha }_t\ne 0$ satisfies Step 7 in the conditional gradient algorithm (Algorithm 1), and

$$\mathbf{X}\in M\Rightarrow Q_t\left(\mathbf{X}\right)+P_t\in M.$$

Equation (9) can be considered a perturbation of the upward direction $d^t$, which is substituted with a new direction $D_k=d^t$$+\frac{P_t}{{\alpha }_t}$. As a result, iterations (9) become:

$${\mathbf{X}}^{t+1}={\mathbf{X}}^t+{\alpha }_t{D}_t.$$

In the literature [18,26,27], general properties can be found to select a sequence suitable for perturbation ${\left\{P_t\right\}}_{t\ge 0}$. Typically, perturbations that satisfy these features are produced using sequences of Gaussian laws.

We define a random vector $Z_t$ and use the symbols ${\mathsf{\Phi }}_t$ and ${\varphi }_t$ to represent its cumulative distribution function and probability density function, respectively.

The conditional probability density function of ${\mathbf{X}}^{t+1}$ is represented by $f_{t+1}$, and the conditional cumulative distribution function is designated as $F_{t+1}(y$| ${\mathbf{X}}^t=x)$.

$$F_{t+1}\left(y\right|{\mathbf{X}}^t=x)=P({\mathbf{X}}^{t+1}<y|{\mathbf{X}}^t=x).$$

We define a sequence of n-dimensional random vectors ${\left\{Z_t\right\}}_{t\ge 0}\in M$. Additionally, we also take into account ${\left\{{\xi }_t\right\}}_{t\ge 0}$, a decreasing sequence of positive real numbers that steadily approaches 0, where ${\xi }_0$ is less than or equal to 1.

Let $P_t={\xi }_t{Z}_t$

$$F_{t+1}\left(y\right|{\mathbf{X}}^t=x)=P({\mathbf{X}}^{t+1}<y|{\mathbf{X}}^t=x).$$

It follows that

$$F_{t+1}\left(y\right|{\mathbf{X}}^t=x)=P\left(Z_t<{\displaystyle \frac{y-Q_t\left(x\right)}{{\xi }_t}}\right)={\mathsf{\Phi }}_t\left({\displaystyle \frac{y-Q_t\left(x\right)}{{\xi }_t}}\right).$$

Therefore, we have

$$f_{t+1}\left(y\right|{\mathbf{X}}^t=x)={\displaystyle \frac{1}{{\xi }_t^n}}{\varphi }_t\left({\displaystyle \frac{y-Q_t\left(x\right)}{{\xi }_t}}\right),y\in M.$$

(10)

Relation (5) shows that

$$∥y-Q_t\left(x\right)∥\le \rho (\phi ,\epsilon )\mathrm{for}(x,y)\in M_{\phi }\times M_{\phi }.$$

We suppose that $t\mapsto h_t\left(t\right)>0$ is a decreasing function defined on ${\mathbb{R}}^{+}$ such that

$$y\in M_{\phi }\Rightarrow {\varphi }_t\left({\displaystyle \frac{y-Q_t\left(x\right)}{{\xi }_t}}\right)\ge h_t\left({\displaystyle \frac{\rho (\phi ,\epsilon )}{{\xi }_t}}\right).$$

(11)

For simplicity, let

$$Z_t={\mathbf{1}}_M\left(Z_t\right)Z_t,$$

(12)

and $Z\sim \mathbf{N}(0,1),$ where Z is a random variable.

The procedure generates a sequence $V_t=f\left({\mathbf{X}}^t\right).$ By construction, this sequence is increasing and upper-bounded by ${\eta }^{*}.$

$$\forall t\ge 0:{\eta }^{*}\ge V_{t+1}\ge V_t.$$

(13)

Thus, there exists V ≤ ${\eta }^{*}$ such that

$$V_t\to Vfort\to +\infty .$$

Lemma 1.

Let $P_t={\xi }_t{Z}_t$ and $\gamma =f\left(x^0\right)$ if $Z_t$ is given by (12). Then, there exists $\ell >0$ such that

$$P(V_{t+1}>\omega |V_t\le \omega )\ge {\displaystyle \frac{meas(M_{\gamma }-M_{\omega })}{{\xi }_t^n}}{h}_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right)>0\forall \omega \in ({\eta }^{*},{\eta }^{*}+\ell ],$$

where $n=dim\left(E\right)$.

Proof. 

Let $M_{\omega }=\left\{x\in M\left|f\right(x)<\omega \right\},$ for $\omega \in ({\eta }^{*},{\eta }^{*}+\ell ].$

Since $M_{\phi }\subset {\widehat{M}}_{\omega },$${\eta }^{*}<\phi <\omega ,$ it can be deduced from (2) that ${\widehat{M}}_{\omega }$ is non-empty and has a strictly positive measure.

If $meas(M-{\widehat{M}}_{\omega })=0$ for any $\omega \in ({\eta }^{*},{\eta }^{*}+\ell ],$ the result is immediate, since we have $f\left(x\right)={\eta }^{*}$ on $M.$

Let us assume that there exists $\epsilon >0$ such that $meas(M-{\widehat{M}}_{\omega })>0.$ For $\omega \in$$({\eta }^{*},{\eta }^{*}+\epsilon ],$ we have ${\widehat{M}}_{\omega }\subset {\widehat{M}}_{\epsilon }$ and $meas(M-{\widehat{M}}_{\omega })>0.$

$P({\mathbf{X}}^t\notin {\widehat{M}}_{\omega })=P({\mathbf{X}}^t\in S-{\widehat{M}}_{\omega })={\int }_{M-{\widehat{M}}_{\omega }}P({\mathbf{X}}^t\in dx)>0$ for any $\omega \in ({\eta }^{*},{\eta }^{*}+\epsilon ]$ and, since the sequence ${\left\{{\mathbf{V}}_i\right\}}_{i\ge 0}$ is increasing, we also have

$${\left\{{\mathbf{X}}^i\right\}}_{i\ge 0}\subset M_{\gamma }.$$

(14)

Thus

$$P({\mathbf{X}}^t\notin {\widehat{M}}_{\omega })=P({\mathbf{X}}^t\in N-{\widehat{M}}_{\omega })={\int }_{M_{\gamma }-{\widehat{M}}_{\omega }}P({\mathbf{X}}^t\in dx)>0\mathrm{for}\mathrm{any}\omega \in ({\eta }^{*},{\alpha }^{*}+\epsilon ].$$

Letting $\omega \in$$({\eta }^{*},{\alpha }^{*}+\epsilon ],$ we have from (13)

$$P(V_{t+1}>\omega |V_t\le \omega )=P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^i\notin {\widehat{M}}_{\omega },i=0,\dots ,t).$$

However, the Markov chain produces

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^i\notin {\widehat{M}}_{\omega },i=0,\dots ,t)=P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega }).$$

By the conditional probability rule,

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })=\frac{P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega },{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })}{P({\mathbf{X}}^t\notin {\widehat{M}}_{\omega })}.$$

Moreover,

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })={\int }_{M-{\widehat{M}}_{\omega }}P({\mathbf{X}}^t\in dx){\int }_{{\widehat{M}}_{\omega }}{f}_{t+1}\left(y\right|{\mathbf{X}}^t=x)dy.$$

From (14), we have

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })={\int }_{M_{\gamma }-{\widehat{M}}_{\omega }}P({\mathbf{X}}^t\in dx){\int }_{{\widehat{M}}_{\omega }}{f}_{t+1}\left(y\right|{\mathbf{X}}^t=x)dy,$$

and

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })\ge \underset{x\in M_{\gamma }-{\widehat{M}}_{\omega }}{inf}\left\{{\int }_{{\widehat{M}}_{\omega }}{f}_{t+1}\left(y\right|{\mathbf{X}}^t=x)dy\right\}{\int }_{M_{\gamma }-{\widehat{M}}_{\omega }}P({\mathbf{X}}^t\in dx).$$

Thus

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })\ge \underset{x\in M_{\gamma }-{\widehat{M}}_{\omega }}{inf}\left\{{\int }_{{\widehat{M}}_{\omega }}{f}_{t+1}\left(y\right|{\mathbf{X}}^t=x)dy\right\}.$$

Taking (10) into account, we have

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })\ge {\displaystyle \frac{1}{{\xi }_t^n}}\underset{x\in M_{\gamma }-{\widehat{M}}_{\omega }}{inf}\left\{{\int }_{{\widehat{M}}_{\omega }}{\varphi }_t\left({\displaystyle \frac{y-Q_t\left(x\right)}{{\xi }_t}}\right)dy\right\}.$$

Relation (5) shows that

$$∥y-Q_t\left(x\right)∥\le \rho (\gamma ,\epsilon ),$$

and (11) yields

$${\varphi }_t\left({\displaystyle \frac{y-Q_t\left(x\right)}{{\xi }_t}}\right)\ge h_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right).$$

Hence,

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })\ge {\displaystyle \frac{1}{{\xi }_t^n}}\underset{x\in M_{\gamma }-{\widehat{M}}_{\omega }}{inf}{\int }_{{\widehat{M}}_{\omega }}{h}_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right)dy.$$

$$P({\mathbf{X}}^{t+1}\in {\widehat{M}}_{\omega }|{\mathbf{X}}^t\notin {\widehat{M}}_{\omega })\ge {\displaystyle \frac{meas(M_{\gamma }-M_{\omega })}{{\xi }_t^n}}{h}_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right).$$

□

The following result, which follows from Borel–Catelli’s lemma (as described in [18], for example), is a consequence of the global convergence:

Lemma 2.

Let ${\left\{V_t\right\}}_{t\ge 0}$ be a increasing sequence, upper-bounded by ${\eta }^{*}$. Then, there exists V such that $V_t\to V$ for $t\to +\infty$. Assume that there exists $\ell >0$ such that, for any $\omega \in ({\alpha }^{*},{\alpha }^{*}+\ell ]$, there is a sequence of strictly positive real numbers ${\left\{{\mathbf{c}}_t\left(\omega \right)\right\}}_{t\ge 0}$, such that

$$\forall t\ge 0:P(V_{t+1}>\omega |V_t\le \omega )\ge {\mathbf{c}}_t\left(\omega \right)>0and\sum _{t=0}^{+\infty }{\mathbf{c}}_t\left(\omega \right)=+\infty .$$

Then $V={\eta }^{*}$ almost surely.

Proof. 

For instance, see [18,28]. □

Theorem 1.

Assuming $x^0$ belongs to M, and letting $\gamma =f\left(x^0\right)$, let the sequence ${\xi }_t$ be non-increasing, and

$$\sum _{t=0}^{+\infty }{h}_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right)=+\infty .$$

(15)

Then $V={\eta }^{*}$ almost surely.

Proof. 

Let

$${\mathbf{c}}_t\left(\omega \right)={\displaystyle \frac{meas(M_{\gamma }-M_{\omega })}{{\xi }_t^n}}{h}_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right)>0.$$

Since the sequence ${\left\{{\xi }_t\right\}}_{t\ge 0}$ is non increasing,

$${\mathbf{c}}_t\left(\omega \right)\ge {\displaystyle \frac{meas(M_{\gamma }-M_{\omega })}{{\xi }_t^n}}{h}_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right)>0.$$

Thus, Equation (15) shows that

$$\sum _{t=0}^{+\infty }{\mathbf{c}}_t\left(\omega \right)\ge {\displaystyle \frac{meas(M_{\gamma }-M_{\omega })}{{\xi }_t^n}}\sum _{t=0}^{+\infty }{h}_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right)=+\infty .$$

We can conclude that $V={\eta }^{*}$ almost surely by applying Lemmas 1 and 2. □

Theorem 2.

Let $Z_t$ be defined by (12) and ${\xi }_t$ by

$${\xi }_t=\sqrt{\frac{\widehat{b}}{log(t+\widehat{a})}},$$

(16)

where $\widehat{b}>0$, $\widehat{a}>0$, and t is the iteration number. If $x^0\in M$, then for $\widehat{b}$ large enough, $V={\eta }^{*}$ almost surely.

Proof. 

We have

$${\varphi }_t\left(Z\right)={\displaystyle \frac{1}{{\left(\sqrt{2\pi }\right)}^n}}exp(-{\displaystyle \frac{1}{2}}{∥Z∥}^2)=h_t\left(∥Z∥\right)>0,$$

so

$$h_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right)={\displaystyle \frac{1}{{\left(\sqrt{2\pi }\right)}^n{(t+\widehat{a})}^{\rho {(\gamma ,\epsilon )}^2/2\widehat{b}}}}.$$

For $\widehat{b}$, such that

$$0<{\displaystyle \frac{\rho {(\gamma ,\epsilon )}^2}{2\widehat{b}}}<1,$$

we have

$$\sum _{t=0}^{+\infty }{h}_t\left({\displaystyle \frac{\rho (\gamma ,\epsilon )}{{\xi }_t}}\right)=+\infty ;$$

furthermore, as per the previous Theorem 2, it can be deduced that V is almost surely equal to ${\eta }^{*}$. □

5. Numerical Results

In this section, we present numerical results of six examples implemented using the CGB method and the perturbed RPCGB method. Our aim is to compare the performance of these two algorithms.

We begin by applying the algorithm to the initial value, which is ${\mathbf{X}}^0=x^0\in M$. At each step $t\ge 0$, ${\mathbf{X}}^t$ is known, and we calculate ${\mathbf{X}}^{t+1}$.

$k_{\mathrm{sto}}$ denotes the number of perturbations. When $k_{\mathrm{sto}}=0$, the method used is the conditional gradient with bisection, without any perturbations (unperturbed conditional gradient with bisection method).

The Gaussian variates used in our experiments are generated using regular generator calls. Specifically, we use

$${\xi }_t=\sqrt{\frac{\widehat{b}}{log(t+2)}},\mathrm{where}\widehat{b}>0.$$

The definitions of the methods listed in the tables are as follows:

(i)

“CGB”, the method of conditional gradient and bisection;

(ii)

“RPCGB”, the method of random perturbation of conditional gradient and bisection.

The proposed RPCGB algorithm is implemented using the MATLAB programming language. We evaluate the performance of the RPCGB method and compare it with the CGB method for high-dimensional problems. We test the efficacy of these algorithms on several problems [29,30,31,32] with linear constraints, using predetermined feasible starting points $x^0$. The results are presented in

Table 1. The results obtained from the CGB and RPCGB algorithms.

Problem 1		Algorithm
		CGB			RPCGB
$\mathit{n}$	${\mathit{n}}_{\mathit{c}}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{f}}_{\mathit{CGB}}^{*}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{k}}_{\mathit{sto}}$	${\mathit{f}}_{\mathit{RPCGB}}^{*}$
500	1000	9	0.11	−1.06 × ${10}^6$	9	0.68	2	−1.61 × ${10}^7$
900	1800	9	0.13	−3.41 × ${10}^6$	4	0.25	2	−6.32 × ${10}^7$
2000	4000	10	3.27	−1.08 × ${10}^7$	6	5.31	5	−3.53 × ${10}^8$
4000	8000	12	19.33	−5.53 × ${10}^7$	7	27.27	10	−1.48 × ${10}^9$
6000	12,000	19	35.70	−1.34 × ${10}^7$	9	46.54	10	−3.35 × ${10}^9$
9000	18,000	27	87.92	−1.61 × ${10}^7$	13	91.16	10	−7.62 × ${10}^9$

,

Table 2. The results obtained from the CGB and RPCGB algorithms.

Problem 2		Algorithm
		CGB			RPCGB
$\mathit{n}$	${\mathit{n}}_{\mathit{c}}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{f}}_{\mathit{CGB}}^{*}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{k}}_{\mathit{sto}}$	${\mathit{f}}_{\mathit{RPCGB}}^{*}$
500	1000	4	0.07	−50	2	0.02	1	−50
900	1800	5	0.09	−90	2	0.05	1	−90
2000	4000	7	0.11	−199.99	3	0.09	1	−200
4000	8000	9	0.19	−400	4	0.12	1	−400
6000	12,000	10	0.37	−599.99	7	0.15	1	−600
9000	18,000	13	0.42	−900	9	0.21	1	−900

,

Table 3. The results obtained from the CGB and RPCGB algorithms.

Problem 3		Algorithm
		CGB			RPCGB
$\mathit{n}$	${\mathit{n}}_{\mathit{c}}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{f}}_{\mathit{CGB}}^{*}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{k}}_{\mathit{sto}}$	${\mathit{f}}_{\mathit{RPCGB}}^{*}$
500	1000	5	0.05	−498.99	3	0.04	1	−499
900	1800	7	0.07	−898.99	6	0.07	1	−898.76
2000	4000	11	0.12	−1475.44	11	0.19	1	−1998.99
4000	8000	18	0.31	−2951.62	12	0.47	1	−3999
6000	12,000	24	0.79	−4427.81	19	0.74	1	−5998.87
9000	18,000	35	1.03	−6642.08	27	0.96	1	−8998.25

,

Table 4. The results obtained from the CGB and RPCGB algorithms.

Problem 4		Algorithm
		CGB			RPCGB
$\mathit{n}$	${\mathit{n}}_{\mathit{c}}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{f}}_{\mathit{CGB}}^{*}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{k}}_{\mathit{sto}}$	${\mathit{f}}_{\mathit{RPCGB}}^{*}$
500	1000	23	11.41	−131.81	45	19.53	25	−176.72
900	1800	29	17.06	−214.79	57	21.34	30	−293.51
2000	4000	34	42.66	−417.95	69	71.26	30	−536.38
4000	8000	56	67.18	−768.22	75	96.02	50	−1.06 × ${10}^3$
6000	12,000	73	79.63	−846.01	94	110.63	70	−1.11 × ${10}^3$
9000	18,000	89	99.25	−919.85	124	136.71	90	−1.35 × ${10}^3$

,

Table 5. The results obtained from the CGB and RPCGB algorithms.

Problem 5		Algorithm
		CGB			RPCGB
n	${\mathit{n}}_{\mathit{c}}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{f}}_{\mathit{CGB}}^{*}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{k}}_{\mathit{sto}}$	${\mathit{f}}_{\mathit{RPCGB}}^{*}$
500	500	19	8.04	4.71 × ${10}^{-4}$	14	9.23	70	1.38 × ${10}^{-11}$
900	900	29	14.13	8.51 × ${10}^{-4}$	14	14.25	70	3.05 × ${10}^{-10}$
2000	2000	42	33.09	0.0097	26	45.31	150	4.96 × ${10}^{-10}$
4000	4000	30	53.63	0.0194	17	97.27	200	7.93 × ${10}^{-7}$
6000	6000	59	71.47	0.0291	19	122.54	300	7.93 × ${10}^{-9}$
9000	9000	77	92.55	0.0436	13	153.16	800	7.93 × ${10}^{-5}$

and

Table 6. The results obtained from the CGB and RPCGB algorithms.

Problem 6		Algorithm
		CGB			RPCGB
n	${\mathit{n}}_{\mathit{c}}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{f}}_{\mathit{CGB}}^{*}$	${\mathit{k}}_{\mathit{iter}}$	CPU	${\mathit{k}}_{\mathit{sto}}$	${\mathit{f}}_{\mathit{RPCGB}}^{*}$
500	499	2	0.09	−1.8856	5	0.28	10	−4.0147
900	899	19	1.65	0.7033	7	0.38	10	−4.4432
1000	999	13	0.92	0.1626	12	1.08	10	−4.6572
2000	1999	21	1.70	2.0535	8	17.24	40	−3.1399
3000	2999	23	2.96	0.9158	9	27.54	60	−3.2145
4000	3999	12	1.47	2.0097	11	48.56	100	−4.6168

and Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6, where n denotes the dimension of the problem under consideration and $n_c$ represents the number of constraints. The reported test results include the optimal value $f_{RPCGB}^{*}$ and the number of iterations $k_{iter}$.

The optimal line search process of CGB and RPCGB is found using the bisection method with $ϵ={10}^{-4}$. We terminate the iterative process when either the best solution (global solution) is found or the maximum number of iterations has been reached.

All algorithms were run on a TOSHIBA Intel(R) Core(TM) CPU running at 2.40 GHz with 6 GB of RAM, a Core i7 processor, and the 64-bit Windows 7 Professional operating system. The “CPU” column in the table displays the mean CPU time for one run in seconds.

Problem 1.

The Neumaier 3 Problem (NF3) is a mathematical optimization problem introduced by Arnold Neumaier in 2003 (see [29]). The problem is defined as follows:

$$\left\{\begin{array}{cc}minimize:& \hfill {\displaystyle \sum _{j=1}^n}{(x_j-1)}^2-{\displaystyle \sum _{j=2}^n}{x}_j{x}_{j-1}\hfill \\ subjectto:& \hfill -n^2\le x_j\le n^2,j=1,\dots ,n\hfill \end{array}\right.$$

Problem 2.

The Cosine Mixture Problem (CM) is an optimization problem introduced by Breiman and Cutler in 1993 (see [29]). The problem is defined as follows:

$$\left\{\begin{array}{cc}minimize:& \hfill {\displaystyle \sum _{j=1}^n}{x}_j^2-0.1{\displaystyle \sum _{j=1}^n}cos\left(5\pi x_j\right)\hfill \\ subjectto:& \hfill -1\le x_j\le 1,j=1,\dots ,n\hfill \end{array}\right.$$

Problem 3.

The Inverted Cosine Wave Function or the Cosine Mixture with Exponential Decay Problem. This is a commonly used benchmark problem in global optimization and was introduced by Price et al. in 2006 (see [30]). The problem is defined as follows:

$$\left\{\begin{array}{cc}minimize:& -{\displaystyle \sum _{j=1}^n}exp\left({\displaystyle \frac{-x_j^2-x_{j+1}^2-0.5x_j{x}_{j+1}}{8}}\right)cos\left(4\sqrt{x_j^2+x_{j+1}^2+0.5x_j{x}_{j+1}}\right)\\ subjectto:& \hfill -5\le x_j\le 5,j=1,\dots ,n\hfill \end{array}\right.$$

Problem 4.

The Epistatic Michalewicz Problem (EM) is a type of optimization problem commonly used as a benchmark in evolutionary computation and optimization. It was introduced by Michalewicz in 1996 (see [29]). The problem is defined as follows:

$$\left\{\begin{array}{cc}minimize:& -{\displaystyle \sum _{j=1}^n}sin\left(y_j\right){(sin({\displaystyle \frac{j{y}_j^2}{\pi }}))}^{20}\\ subjectto:& 0\le x_j\le \pi ,j=1,\dots ,n\\ y_j=& \left\{\begin{array}{c}{x}_{j}cos\left({\displaystyle \frac{\pi }{6}}\right)-x_{j+1}sin\left({\displaystyle \frac{\pi }{6}}\right),j=1,3,5,\dots ,\le n\\ \\ x_{j}sin\left({\displaystyle \frac{\pi }{6}}\right)+x_{j+1}cos\left({\displaystyle \frac{\pi }{6}}\right),j=2,4,6,\dots ,\le n\\ \\ x_j,j=n\end{array}\right.\end{array}\right.$$

Problem 5.

The problem is a mathematical optimization problem used in global optimization which comes from [32] and is defined as follows:

$$\left\{\begin{array}{cc}minimize:& \hfill {\displaystyle \sum _{j=1}^n}(x_j^2-10cos\left(2\pi x_j\right)+10)\hfill \\ subjectto:& \hfill {\displaystyle \sum _{j=1}^n}{x}_j=0,\hfill \\ & -5.12\le x_j\le 5.12,j=1,\dots ,n\end{array}\right.$$

Problem 6.

Rastrigin’s function is a non-convex, multi-modal function commonly used as a benchmark problem in optimization. It was introduced by Rastrigin in 1974 (see [31]) and is defined as:

$$\left\{\begin{array}{cc}minimize:& \hfill {\displaystyle \sum _{j=1}^n}cos(2\pi x_{j}sin\left({\displaystyle \frac{\pi }{20}}\right))\hfill \\ subjectto:& \hfill x_j-x_{j+1}=0.4j=1,\dots ,n-1\hfill \end{array}\right.$$

To gain a deeper understanding of the effect of the modifications on the proposed algorithm, we utilized a scatter plot that illustrates the distribution of the algorithm’s solutions in two dimensions for both the CGB and RPCGB algorithms. The goal was to generate scatter plots to depict the distribution of solutions in a 2D space for all problems when using two variables (n = 2). However, we found that for Problems 1 to 3, we were able to obtain the optimal solution value using only one iteration, which made the creation of a scatter plot unnecessary in this case. Thus, we generated a scatter plot of the solution distribution for the case of n = 900. This allowed us to effectively illustrate the distribution of solutions using scatter plots. In each algorithm, the scatter plot was generated at the first iteration and continued up to the limit of the required number of iterations to reach the solution. After analyzing Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6, we concluded that the modified algorithm exhibited a more tightly clustered solution distribution in the scatter plot compared to the original algorithm.

We also present in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6 the results of plotting the objective function values in each iteration and the convergence performance for the CGB and RPCGB methods with 9000 variables. The plots (d) show that the proposed algorithm performs better than the CGB algorithm, with the majority of cases showing that the suggested algorithm achieves convergence in fewer iterations than the CGB algorithm. However, there is an exception observed in Problem 4, as presented in Figure 4, where the CGB algorithm stops early. This demonstrates that the convergence behavior of optimization algorithms can vary based on the problem being solved. It is worth noting that both algorithms terminated their execution before the 30th iteration, which is because a stopping criterion of approximately $ϵ={10}^{-4}$ was met. The algorithms cease their iterations upon reaching the optimal solution (the local or global solution) or upon reaching the maximum number of iterations. We observe that the random perturbation has a significant effect on the convergence. This suggests that the changes made to the algorithm led to an improvement in its performance.

The results presented in Table 1, Table 2, Table 3, Table 4, Table 5 and Table 6 demonstrate that the CGB algorithm is capable of obtaining global solutions in certain instances regardless of the number of dimensions, such as Problem 2 (see Table 2). However, in some cases, the CGB algorithm fails to obtain global solutions as the number of dimensions increases, as seen in Problem 3 (see Table 3). In contrast, our RPCGB algorithm can obtain a global solution for all cases, and the computational results indicate that it performs effectively for these high-dimensional problems. These results also indicate that, for larger-dimensions problems, the CGB method necessitates a greater number of iterations to finalize the optimization process, whereas the RPCGB method does not, as evidenced by Problem 5. This difference can be explained by the $k_{\mathrm{sto}}$ parameter, which denotes the number of perturbations. If the number of perturbations is raised, then the number of iterations needed to achieve the optimal solution is reduced.

When analyzing the obtained results, it is evident that the perturbed conditional gradient method with bisection algorithm (RPCGB) performs well compared to the conditional gradient algorithm (CGB).

6. Conclusions

In this work, we generalized the RPCGB method to solve large-scale non-convex optimization problems. The algorithms mentioned in this paper, specifically the conditional gradient algorithm, are commonly used optimization techniques for solving convex optimization problems. However, in the case of non-convex optimization problems, these algorithms may converge to a local minimum instead of the global minimum. To overcome this problem, the proposed approach introduces a random perturbation to the optimization problem. Specifically, at each iteration of the algorithm, a random perturbation is added to the $Q_t$ operator, which allows the algorithm to escape from local minima and explore the search space more effectively. The bisection algorithm is used to find the optimal step size along the search direction. It involves solving a one-dimensional optimization problem to find the step size that minimizes the objective function. By combining these two algorithms with the random perturbation approach, the proposed method is able to efficiently explore the search space for large-scale non-convex optimization problems under linear constraints and reach a global minimum. The tuning of the parameters $k_{sto}$ and $\widehat{b}$ is related to the main difficulty in applying random perturbation in practice.

The RPCGB algorithm has the ability to solve various problems, such as control systems, as well as optimization problems in machine learning, robotics, and image reconstruction. There are problems that contain a part that is not smooth. Therefore, in the future, we plan to use the random perturbation of the conditional subgradient method with bisection algorithm to solve non-convex, non-smooth (non-differentiable) programming under linear constraints. Additionally, we intend to use the perturbed conditional gradient method to address non-convex optimization problems in support vector machines (SVM).