Bundle Enrichment Method for Nonsmooth Difference of Convex Programming Problems

Abstract

The Bundle Enrichment Method (BEM-DC) is introduced for solving nonsmooth difference of convex (DC) programming problems. The novelty of the method consists of the dynamic management of the bundle. More specifically, a DC model, being the difference of two convex piecewise affine functions, is formulated. The (global) minimization of the model is tackled by solving a set of convex problems whose cardinality depends on the number of linearizations adopted to approximate the second DC component function. The new bundle management policy distributes the information coming from previous iterations to separately model the DC components of the objective function. Such a distribution is driven by the sign of linearization errors. If the displacement suggested by the model minimization provides no sufficient decrease of the objective function, then the temporary enrichment of the cutting plane approximation of just the first DC component function takes place until either the termination of the algorithm is certified or a sufficient decrease is achieved. The convergence of the BEM-DC method is studied, and computational results on a set of academic test problems with nonsmooth DC objective functions are provided.

1. Introduction

Optimization approaches are essential for solving a wide range of practical problems. There are various problems based on these approaches including unconstrained and constrained problems, problems with linear and nonlinear as well as smooth and nonsmooth objective and/or constraint functions, and problems with continuous and integer decision variables [1]. The majority of optimization problems from applications have special structures (for example, convexity) which can be exploited to design efficient and accurate methods for their solutions. Difference of convex (DC) optimization problems are among such problems, where the objective and/or constraint functions are represented as a difference of two convex functions.

In this research work, we consider the unconstrained nonsmooth DC programming problem

$$\underset{x\in {\mathbb{R}}^n}{min}f\left(x\right),$$

(1)

where $f:{\mathbb{R}}^n\to \mathbb{R}$ is, in general, nonsmooth and is expressed as a difference of two convex functions $f_1,f_2:{\mathbb{R}}^n\to \mathbb{R}$:

$$f\left(x\right)=f_1\left(x\right)-f_2\left(x\right).$$

Here, $f_1-f_2$ is called a DC representation (decomposition) of f while $f_1$ and $f_2$ are DC components of the function $f.$ The DC components $f_1$ and $f_2$ are, in general, nonsmooth [2,3,4,5].

Nonsmooth DC programming is an important subclass of DC optimization problems [6], and many practical problems are modeled as a DC programming problem. They include the bridge location problem, the design centering problem [7], the packing problem [8], the production–transportation planning problem [9], the location planning problem [10], the edge detection problem [11], the conic programming problem [12], cluster analysis [13], and regression analysis [14]. Recently, DC optimization problems with uncertain data has become an interesting topic, and the results from robust optimization, in particular those obtained in [15,16,17], can be extended to robust DC optimization.

DC optimization problems have been studied in the context of both local and global problems, and various methods have been developed for solving these problems globally [7,18]. To the best of our knowledge, the first local search algorithm for solving DC optimization problems is the difference of convex algorithm (DCA) introduced in [19] and further explored, for example, in [8,20]. Since then, the development of local DC optimization methods for solving Problem (1) has attracted remarkable scholarly attention. Next, we provide a short description of such methods and give references for more details. These methods can be classified into three categories:

• The first category consists of the DCA and its modifications. The basic idea of the DCA is to linearize the concave part $-f_2$ around the current iterate by using its subgradient while keeping the convex part $f_1$ as it is in the minimization process. To improve the convergence of the DCA, various modifications have been developed, for instance, in [21,22,23]. The boosted DC algorithm (BDCA), proposed in [21,22], accelerates the convergence of the DCA by using an additional line search step. The inertial DCA, introduced in [23], defines trial points whose sequence of functional values is not necessarily monotonically decreasing. This controls the algorithm from converging to a critical point that is not d-stationary. In [24], the BDCA is combined with a simple derivative-free optimization method. This allows one to force the d-stationarity (lack of descent direction) at the obtained point. To avoid the difficulty of solving the DCA’s subproblem, in [25], the first DC component is replaced with a convex model, and the second DC component is used without any approximation;
• The methods in the second category, which we refer to as DC-Bundle, are various extensions of the bundle methods for convex problems. The piecewise linear underestimates or subgradients of both DC components are utilized to compute search directions. The methods include the codifferential method [26], the proximal bundle method with the concave affine model [27,28], the proximal bundle method for DC optimization [29,30], the proximal point algorithm [31], the proximal linearized algorithm [32], the double bundle method [33], and the nonlinearly constrained DC bundle method [34];
• The methods in the third category are those that use the convex piecewise linear model of the first DC component and one subgradient of the second DC component to compute search directions at each iteration [35,36]. They differ from those in the first category as at each iteration of these methods, the model of the first DC component is updated and the new subgradient of the second component is calculated. They also differ from those in the second category as they use only one subgradient of the second DC component at each iteration whereas the DC-Bundle methods may use more than one subgradient of this component to build its piecewise linear model. Note that some overlapping is unavoidable when classifying methods into different categories: for instance, the proximal bundle method for DC optimization [30] may be classified into both the second and the third categories. The other methods in the third category include the aggregate subgradient method for DC optimization [35] and the augmented subgradient method [36]. In the former method, the aggregate subgradient of the first DC component and one subgradient of the second component are utilized to compute search directions. In the latter method, augmented subgradients of convex functions are defined and used to model the first DC component. Then, this model and one subgradient of the second DC component are used for computing search directions.

The proposed BEM-DC method belongs to the DC-bundle family whose main features, with respect to the DCA, are

• developing a model function based on cutting-plane approximations of $f_1$ and, possibly, $f_2$;
• adding a regularization term (typically a proximity one) to the model for stabilization purposes.

We assume the familiarity of the reader with the basic notions of bundle methods. We refer to the following books and articles from the vast literature available [37,38,39,40,41,42,43,44,45]. In our iterative scheme, we adopt two cutting-plane approximations (based, as usual in bundle methods, on information coming from previous iterations) for $f_1$ and $f_2$, respectively. Consequently, we have a model which is still DC, being the difference of two convex piecewise affine functions. Similarly to  [29,33], we tackle the (global) minimization of the model by solving a set of convex problems whose cardinality depends on the number of linearizations adopted to approximate the function $f_2$. The novelty of the approach consists of the dynamic management of the bundle. This enables us to achieve a parsimonious use of the information available, in view of reducing the computational burden to minimize the model function at each iteration. The adopted strategies are based on the following:

• The information coming from the previous iterates contributes to the approximation of exactly one of the DC components $f_1$ and $f_2$, depending on the sign of the linearization error relative to the current iterate. During the iterative process, every time a new point is generated, subgradients of both component functions $f_1$ and $f_2$ are calculated, but only one of them will enter the calculation of the search direction: the choice being driven, at each of the successive iterations, by the sign of the linearization error with respect to the current point. Such sign may change at each iteration. Therefore, we define a dynamic bundling strategy, which implies a consistent reduction (approximately one-half) of the size of the auxiliary problems to be solved;
• If the displacement suggested by the model minimization provides no sufficient decrease of the objective function (the null step in bundle parlance), the temporary enrichment exclusively of the cutting-plane approximation of $f_1$ takes place.

Note that a significant difference between the proposed approach with those introduced in [29,33] is in the reduction of the size of the bundle of the function $f_2$. This has a strong impact on the solution of the auxiliary problem to be solved at each iteration.

The structure of the paper is organized as follows. Section 2 provides necessary notations and some preliminaries. Section 3 presents the new model function formulation. Section 4 describes the new method, and its convergence is discussed in Section 5. Section 6 reports the results of numerical experiments. Section 7 provides some concluding remarks.

2. Notations and Background

First, we provide some notations and definitions that we will use throughout the paper. The inner product in ${\mathbb{R}}^n$ is $\langle u,v\rangle ={\sum }_{i=1}^n{u}_i{v}_i$, and $\parallel ·\parallel$ is the associated norm. For $x\in {\mathbb{R}}^n$ and $\epsilon >0$, $B(x;\epsilon )$ is an open ball of the radius $\epsilon >0$ centered at x. The function $f:{\mathbb{R}}^n\to \mathbb{R}$ is locally Lipschitz continuous on ${\mathbb{R}}^n$ if for every $x\in {\mathbb{R}}^n$, there exists a Lipschitz constant $L>0$ and $\epsilon >0$ such that $\left|f\right(y)-f(z\left)\right|\le L\parallel y-z\parallel$ for all $y,z\in B(x;\epsilon ).$

For a convex function $f:{\mathbb{R}}^n\to \mathbb{R}$, its subdifferential at a point $x\in {\mathbb{R}}^n$ is [38,41]

$$\partial f\left(x\right)=\left\{\xi \in {\mathbb{R}}^n:f\left(y\right)-f\left(x\right)\ge \langle \xi ,y-x\rangle \mathrm{for}\mathrm{all}y\in {\mathbb{R}}^n\right\},$$

and for $\epsilon >0$, its $\epsilon$-subdifferential is

$${\partial }_{\epsilon }f\left(x\right)=\left\{{\xi }_{\epsilon }\in {\mathbb{R}}^n:f\left(y\right)-f\left(x\right)\ge \langle {\xi }_{\epsilon },y-x\rangle -\epsilon \mathrm{for}\mathrm{all}y\in {\mathbb{R}}^n\right\}.$$

Each vector $\xi \in \partial f\left(x\right)$(${\xi }_{\epsilon }\in {\partial }_{\epsilon }f\left(x\right)$) is called a subgradient ($\epsilon$-subgradient) of f at x.

A point $x^{*}\in {\mathbb{R}}^n$ is called a critical point of Problem (1) if

$$\partial f_1\left(x^{*}\right)\cap \partial f_2\left(x^{*}\right)\ne \varnothing ,$$

and a point ${\overline{x}}^{*}\in {\mathbb{R}}^n$ is said to be an $\epsilon$-critical point if [31]

$$\begin{array}{c}{\partial }_{\epsilon }{f}_1\left({\overline{x}}^{*}\right)\cap {\partial }_{\epsilon }{f}_2\left({\overline{x}}^{*}\right)\ne \varnothing .\end{array}$$

Next, we recall the basic idea of the standard cutting-plane model for any convex nonsmooth function $f.$ Let $x_k\in {\mathbb{R}}^n$ be the current iteration point, $x_j\in {\mathbb{R}}^n$ be some auxiliary points (from past iterations), ${\xi }_j\in \partial f\left(x_j\right)$ be the subgradients of the function f computed at the point $x_j\in {\mathbb{R}}^n$ for $j\in J^k$, and $J^k$ is a nonempty subset of $\{1,\dots ,k\}$. The cutting-plane model for the function f can be given by

$${\widehat{f}}^k\left(x\right)=\underset{j\in J^k}{max}{\overline{f}}_j\left(x\right),$$

where ${\overline{f}}_j\left(x\right)=f\left(x_j\right)+\langle {\xi }_j,x-x_j\rangle .$ The linearization error

$${\alpha }_j^k=f\left(x_k\right)-{\overline{f}}_j\left(x_k\right)\mathrm{for}\mathrm{all}j\in J^k$$

defines how well ${\overline{f}}_j$ approximates the function f at the current iteration point $x_k.$

In this paper, we assume that f is the DC function, and due to the lack of its convexity, the straightforward application of the cutting-plane approach is meaningless. Nevertheless, it can be separately applied to model the DC components $f_i,i=1,2.$ Thus, we have

$${\widehat{f}}_i^k\left(x\right)=\underset{j\in J^k}{max}\left\{f_i\left(x_j\right)+\langle {\xi }_{i,j},x-x_j\rangle \right\}\mathrm{for}x\in {\mathbb{R}}^n.$$

Here, ${\xi }_{i,j}\in \partial f_i\left(x_j\right),i=1,2$ are the subgradients of the DC components $f_i$ computed at the auxiliary point $x_j\in {\mathbb{R}}^n$ for $j\in J^k.$ This approximation can be rewritten as

$${\widehat{f}}_i^k\left(x\right)=\underset{j\in J^k}{max}\left\{f_i\left(x_k\right)+\langle {\xi }_{i,j},x-x_k\rangle -{\alpha }_{i,j}^k\right\},$$

(2)

where $x_k\in {\mathbb{R}}^n$ is the current iteration point, and ${\alpha }_{i,j}^k$, $i=1,2$ are the linearization errors associated with the j-th first order expansion of $f_i,i=1,2$ rooted at the point $x_j,$ given by

$${\alpha }_{i,j}^k=f_i\left(x_k\right)-f_i\left(x_j\right)-\langle {\xi }_{i,j},x_k-x_j\rangle \mathrm{for}\mathrm{all}j\in J^k.$$

Suppose that information coming from some auxiliary points $x_j\in {\mathbb{R}}^n$, $j\in J^k$, along with $\left(x_k,{\xi }_{1,k}\in \partial f_1\left(x_k\right),{\xi }_{2,k}\in \partial f_2\left(x_k\right)\right)$, is available. We condense all such information into a bundle set $B^k$ defined as a set of tuples, one for each point $x_j.$ That is,

$$b\left(x_j\right)=\left(x_j,f_1\left(x_j\right),f_2\left(x_j\right),{\xi }_{1,j}\in \partial f_1\left(x_j\right),{\xi }_{2,j}\in \partial f_2\left(x_j\right),{\alpha }_{1,j}^k,{\alpha }_{2,j}^k\right).$$

In general, we assume that for some appropriate index j, the tuple

$$b\left(x_k\right)=\left(x_k,f_1\left(x_k\right),f_2\left(x_k\right),{\xi }_{1,k}\in \partial f_1\left(x_k\right),{\xi }_{2,k}\in \partial f_2\left(x_k\right),0,0\right)$$

(3)

associated with $x_k$ is in the bundle too.

3. The New Model Function

In this section, we formulate our new model function. First, we distribute the bundle index set $J^k$ into the subsets

$$\begin{array}{c}{J}_1^k=\left\{j:j\in J^k,{\alpha }_j^k\le 0\right\}\mathrm{and}{J}_2^k=\left\{j:j\in J^k,{\alpha }_j^k\ge 0\right\},\end{array}$$

(4)

where $J^k=J_1^k\cup J_2^k.$ It is worth noting that this does not properly define a partition of $J^k$ as the indexes corresponding to ${\alpha }_j^k=0$ are in both subsets $J_1^k$ and $J_2^k$. In addition, ${\alpha }_j^k={\alpha }_{1,j}^k-{\alpha }_{2,j}^k$ can take any sign due to the nonconvexity of the function $f.$ Based on the sign of ${\alpha }_j^k$, we extract elements from the set $J^k$ and modify the definition of the cutting-plane models, given in (2), accordingly. It is clear that ${\alpha }_j^k\le 0$ corresponds to ${\alpha }_{1,j}^k\le {\alpha }_{2,j}^k$. That is, the linearization error at $x_k$ associated with the linearization of the function $f_1$ rooted at $x_j$ is not bigger than that associated with the function $f_2$. This means that the information provided by $x_j$ is more suited to approximate $f_1$ than $f_2$ around $x_k$. The reverse holds for the case ${\alpha }_j^k\ge 0.$

Remark 1. 

The structures of the subsets $J_i^k,i=1,2$ depend on the point $x_k$. Thus, they should be updated every time a new iterate $x_{k+1}$ is calculated.

Consider the subsets $J_i^k,i=1,2.$ Aiming at a parsimonious use of the available information, we restrict the definition of the cutting-plane functions, given in (2), as

$${\widehat{f}}_i^k\left(x\right)=\underset{j\in J_i^k}{max}\left\{f_i\left(x_k\right)+\langle {\xi }_{i,j},x-x_k\rangle -{\alpha }_{i,j}^k\right\},i=1,2,$$

which reduces the number of affine pieces defining the convex approximations of the functions $f_i,i=1,2.$

Next, we introduce the variable $d=x-x_k$ and define the function $h^k\left(d\right)$ (the model of the difference function $f(x_k+d)-f\left(x_k\right)$) as

$$h^k\left(d\right)=h_1^k\left(d\right)-h_2^k\left(d\right)=\underset{j\in J_1^k}{max}\left(\langle {\xi }_{1,j},d\rangle -{\alpha }_{1,j}^k\right)-\underset{j\in J_2^k}{max}\left(\langle {\xi }_{2,j},d\rangle -{\alpha }_{2,j}^k\right),$$

and calculate

$$w^k=\underset{d}{min}{w}^k\left(d\right)=\underset{d}{min}{h}^k\left(d\right)+\frac{1}{2}{\delta }_k{\parallel d\parallel }^2.$$

(5)

Here, ${\delta }_k$ is the proximity parameter used in most bundle methods, and $\frac{1}{2}{\delta }_k{\parallel d\parallel }^2$ is a stabilizing term used to guarantee the existence of the solution. Letting $d_k$ be the solution to (5), we take $x_{k+1}=x_k+d_k$ as the (tentative) new iterate point. Then, we check this point for a possible sufficient decrease in the objective function f as follows:

$$f(x_k+d_k)\le f\left(x_k\right)+\mu v^k,$$

(6)

where $\mu \in (0,1)$ is a given parameter, and $v^k$ is the predicted reduction of the function f provided by the model function $h^k\left(d\right).$ That is, $v^k=v_1^k-v_2^k,$ where

$$\begin{array}{cc}\hfill v_1^k=h_1^k\left(d_k\right)=& \underset{j\in J_1^k}{max}\langle {\xi }_{1,j},d_k\rangle -{\alpha }_{1,j}^k,\mathrm{and}\hfill \\ \hfill v_2^k=h_2^k\left(d_k\right)=& \underset{j\in J_2^k}{max}\langle {\xi }_{2,j},d_k\rangle -{\alpha }_{2,j}^k.\hfill \end{array}$$

Note that we have $v^k\le 0.$ It is sufficient to observe that both functions $h_1^k\left(d\right)$ and $h_2^k\left(d\right)$ are non-positive at $d=0$. If (6) is satisfied (serious step in bundle method parlance), then the point $x_{k+1}=x_k+d_k$ becomes the new estimate of a minimum. Thus, the related information to this point is appended to the bundle, and the bundle index set is updated. Next, the elements of the new index set $J^{k+1}$ are distributed between $J_1^{k+1}$ and $J_2^{k+1}$ according to (4) (see Remark 1), and the procedure is iterated.

If the condition (6) is not satisfied, that is $f(x_k+d_k)>f\left(x_k\right)+\mu v^k,$ then there is no sufficient decrease (null step). Taking into account $\mu \in (0,1)$ and $v^k\le 0$, it follows that

$$\begin{array}{cc}\hfill f_1(x_k+d_k)-f_2(x_k+d_k)& >f_1\left(x_k\right)-f_2\left(x_k\right)+\mu (v_1^k-v_2^k)\hfill \\ & >f_1\left(x_k\right)-f_2\left(x_k\right)+v_1^k-v_2^k,\hfill \end{array}$$

which in turn implies that

$$\left(f_1(x_k+d_k)-f_1\left(x_k\right)\right)-v_1^k>\left(f_2(x_k+d_k)-f_2\left(x_k\right)\right)-v_2^k.$$

(7)

This means that whenever a sufficient decrease does not occur, then the gap between the actual and the predicted reduction for the function $f_1$ is bigger than that of $f_2$. Such an observation is at the basis of our bundle enlargement strategy in the case of the null step. More precisely, whenever a null step occurs, we invest in improving the approximation of $f_1$ more than in that of $f_2$.

The following proposition ensures that in the case of the null step, by inserting the point $x_{k+1}=x_k+d_k$ into the bundle, any couple $\left({\xi }_1^{+},{\alpha }^{+}\right)$ generates a substantial cut of the epigraph of $h_1^k$ and, thus, an improved model of the function $f_1$. Here, ${\xi }_1^{+}\in \partial f_1(x_k+d_k)$ and

$$\begin{array}{c}{\alpha }^{+}=f_1\left(x_k\right)-f_1(x_k+d_k)+\langle {\xi }_1^{+},d_k\rangle .\end{array}$$

Proposition 1. 

Assume that $v^k<-\eta$ for some $\eta >0$, and the sufficient decrease condition (6) is not fulfilled at the point $x_k+d_k$. Then, for any ${\xi }_1^{+}\in \partial f_1(x_k+d_k)$, we have

$$\langle {\xi }_1^{+},d_k\rangle -{\alpha }^{+}>v_1^k+\rho \eta ,$$

where $\rho =1-\mu$.

Proof. 

From the definition of ${\alpha }^{+}$, it follows

$$\langle {\xi }_1^{+},d_k\rangle -{\alpha }^{+}=f_1(x_k+d_k)-f_1\left(x_k\right).$$

Further, from (7), considering the definition of $v_2^k$ and $v^k=v_1^k-v_2^k\le 0$, we have

$$\begin{array}{cc}\hfill \langle {\xi }_1^{+},d_k\rangle -{\alpha }^{+}& =f_1(x_k+d_k)-f_1\left(x_k\right)\hfill \\ & >f_2(x_k+d_k)-f_2\left(x_k\right)+\mu (v_1^k-v_2^k)\hfill \\ & \ge v_2^k+\mu (v_1^k-v_2^k)=v_1^k+\rho (v_2^k-v_1^k)\hfill \\ & =v_1^k-\rho v^k\ge v_1^k+\rho \eta .\hfill \end{array}$$

This completes the proof.    □

4. The Proposed BEM-DC Algorithm

In this section, we describe the BEM-DC algorithm for solving Problem (1) and give its step-by-step format. The algorithm has both inner and outer iterations. The inner iteration contains the evaluation, in terms of the decrease in the objective function, of tentative displacements from the current point. Null steps might occur within the inner iteration whenever a sufficient decrease is not achieved. Once, instead, such reduction is obtained (serious step), the current estimate of the minimum is updated and a new outer iteration takes place. Based on Proposition 1, a subgradient accumulation process takes place any time a null step occurs within the inner iteration. The specific feature of the BEM-DC algorithm is such that a process involves only information about the function $f_1$ which is stored in a temporary bundle of tuples, thus exclusively enriching the bundle $J_1^k$.

Let us denote the outer iteration counter by k. We use l to count for the l-th inner iteration within the k-th outer iteration. Denote the current temporary bundle by $T{B}_l$ and its corresponding index set by $T{J}_l.$ Then, the displacement finding subproblem (5) at the inner iteration l within the k-th iteration takes the following form:

$$w_l^k=\underset{d}{min}{w}_l^k\left(d\right).$$

Define

$$\begin{array}{cc}\hfill d_l^k& =arg\underset{d}{min}{w}_l^k\left(d\right)=arg\underset{d}{min}{h}_l^k\left(d\right)+\frac{1}{2}{\delta }_k{\parallel d\parallel }^2,\mathrm{and}\hfill \\ \hfill v_l^k& =h_l^k\left(d_l^k\right).\hfill \end{array}$$

(8)

Then, considering the presence of two distinct bundles for the function f, we have

$$h_l^k\left(d\right)=h_{1,l}^k\left(d\right)-h_2^k\left(d\right),$$

where

$$\begin{array}{cc}\hfill h_{1,l}^k\left(d\right)& =max\{\underset{j\in J_1^k}{max}\left(\langle {\xi }_{1,j}^k,d\rangle -{\alpha }_{1,j}^k\right),\underset{j\in T{J}_l}{max}\left(\langle {\xi }_{1,j}^k,d\rangle -{\alpha }_{1,j}^k\right)\}\hfill \\ & =\underset{j\in J_1^k\cup T{J}_l}{max}\left(\langle {\xi }_{1,j}^k,d\rangle -{\alpha }_{1,j}^k\right),\mathrm{and}\hfill \\ \hfill h_2^k\left(d\right)& =\underset{j\in J_2^k}{max}\left(\langle {\xi }_{2,j}^k,d\rangle -{\alpha }_{2,j}^k\right).\hfill \end{array}$$

The sets $T{B}_l$ and $T{J}_l$ are updated at each inner iteration, while they are reset to the empty sets at the beginning of each outer iteration. Note that the set $J_2^k$ remains unchanged during the inner iteration process.

Let us focus on Problem (8). Note that it is a DC programming problem whose global optimal solution can be found by solving $|J_2^k|$ convex problems (see [29,33]). Consider the following (convex) problem $P_{l,j}^k$, $j\in J_2^k$:

$$\underset{d}{min}{w}_{l,j}^k\left(d\right)$$

with

$$w_{l,j}^k\left(d\right)=h_{1,l}^k\left(d\right)-\left(\langle {\xi }_{2,j}^k,d\rangle -{\alpha }_{2,j}^k\right)+\frac{1}{2}{\delta }_k{\parallel d\parallel }^2.$$

Let $d_{l,j}^k=arg{min}_d{w}_{l,j}^k\left(d\right)$. It is clear that, defining

$$j^{*}=arg\underset{j\in J_2^k}{min}{w}_{l,j}^k\left(d_{l,j}^k\right),$$

the global optimal solution of Problem (8) is

$$\begin{array}{cc}\hfill d_l^k& =d_{l,j^{*}}^k=arg\underset{d}{min}{w}_{l,j^{*}}^k\left(d\right),\mathrm{and}\hfill \\ \hfill v_l^k& =h_l^k\left(d_l^k\right)=h_{1,l}^k\left(d_l^k\right)-\left(\langle {\xi }_{2,j^{*}}^k,d_l^k\rangle -{\alpha }_{2,j^{*}}^k\right).\hfill \end{array}$$

Summing up, $d_l^k$ is a global optimal solution of Problem (8), and the couple $(d_l^k,v_l^k)$ is the unique optimal solution of the convex problem $P_{l,j^{*}}^k$

$$\begin{array}{ccc}\hfill \underset{d\in {\mathbb{R}}^n,v\in \mathbb{R}}{min}& v+\frac{1}{2}{\delta }_k{\parallel d\parallel }^2\hfill & \\ & v\ge \langle {\xi }_{1,j}^k-{\xi }_{2,j^{*}}^k,d\rangle -\left({\alpha }_{1,j}^k-{\alpha }_{2,j^{*}}^k\right),\hfill & j\in J_1^k\cup T{J}_l.\end{array}$$

(9)

Then, applying the standard duality arguments to Problem (9), we have the following primal-dual relations:

$$\begin{array}{cc}{d}_l^k=& -\frac{1}{{\delta }_k}\left({\displaystyle \sum _{j\in J_1^k\cup T{J}_l}}{\lambda }_j{\xi }_{1,j}^k-{\xi }_{2,j^{*}}^k\right),\end{array}$$

(10)

$$\begin{array}{cc}{v}_l^k=& -\frac{1}{{\delta }_k}\parallel {\displaystyle \sum _{j\in J_1^k\cup T{J}_l}}{\lambda }_j{\xi }_{1,j}^k-{\xi }_{2,j^{*}}^k{\parallel }^2-{\displaystyle \sum _{j\in J_1^k\cup T{J}_l}}{\lambda }_j{\alpha }_{1,j}^k+{\alpha }_{2,j^{*}}^k,\end{array}$$

(11)

$$\begin{array}{cc}{w}_l^k=& -\frac{1}{2{\delta }_k}\parallel {\displaystyle \sum _{j\in J_1^k\cup T{J}_l}}{\lambda }_j{\xi }_{1,j}^k-{\xi }_{2,j^{*}}^k{\parallel }^2-{\displaystyle \sum _{j\in J_1^k\cup T{J}_l}}{\lambda }_j{\alpha }_{1,j}^k+{\alpha }_{2,j^{*}}^k.\end{array}$$

(12)

Here, ${\lambda }_j\ge 0$, $j\in J_1^k\cup T{J}_l$, with ${\sum }_{j\in J_1^k\cup T{J}_l}{\lambda }_j=1$, are the optimal variables of the dual of Problem (9). In addition, the definition of the set of problems $P_{l,j}^k$, $j\in J_2^k$ implies that

$$w_l^k=\underset{d}{min}{w}_l^k\left(d\right)=w_{l,j^{*}}^k\left(d_{l,j^{*}}^k\right)\le w_{l,j}^k\left(d_{l,j}^k\right),j\in J_2^k.$$

Since in the bundle index set $J_2^k$, there exists an index, say $j_0$, associated with the tuple $b\left(x_k\right)$ (see (3)) for which ${\xi }_{2,j_0}^k\in \partial f_2\left(x_k\right)$ and ${\alpha }_{2,j_0}^k=0$, we have

$$w_l^k\le w_{l,j_0}^k=w_{l,j_0}^k\left(d_{l,j_0}^k\right)=-\frac{1}{2{\delta }_k}\parallel \sum _{j\in J_1^k\cup T{J}_l}{\lambda }_j^0{\xi }_{1,j}^k-{\xi }_{2,j_0}^k{\parallel }^2-\sum _{j\in J_1^k\cup T{J}_l}{\lambda }_j^0{\alpha }_{1,j}^k,$$

(13)

for ${\lambda }_j^0\ge 0$, $j\in J_1^k\cup T{J}_l$, with ${\sum }_{j\in J_1^k\cup T{J}_l}{\lambda }_j^0=1$. This together with (12) suggests a possible termination criterion for the proposed algorithm. In fact, whenever $w_l^k\ge -\eta$, we have

$$\frac{1}{2{\delta }_k}\left(\parallel \sum _{j\in J_1^k\cup T{J}_l}{\lambda }_j^0{\xi }_{1,j}^k-{\xi }_{2,j_0}^k{\parallel }^2\right)+\sum _{j\in J_1^k\cup T{J}_l}{\lambda }_j^0{\alpha }_{1,j}^k\le \eta ,$$

which indicates that the subgradient ${\xi }_{2,j_0}^k\in \partial f_2\left(x_k\right)$ is at a distance not bigger than $\sqrt{2{\delta }_k\eta }$ from the $\epsilon$-subdifferential of $f_1$ at $x_k$, where

$$\epsilon =\sum _{j\in J_1^k\cup T{J}_l}{\lambda }_j^0{\alpha }_{1,j}^k\le \eta .$$

This property can be interpreted as the approximate satisfaction at the point $x_k$ of criticality. Nevertheless, we implement the termination test of the proposed algorithm in the more common form of $v_l^k\ge -\eta$. In fact, from (10) and (11), $v_l^k\ge -\eta$ implies $w_l^k\ge -\eta$. Further, the fulfillment of the condition $w_{l,j_0}^k\ge -\eta$ provides an alternative termination criterion as

$$v_{l,j_0}^k=-\frac{1}{{\delta }_k}\left(\parallel \sum _{j\in J_1^k\cup T{J}_l}{\lambda }_j^0{\xi }_{1,j}^k-{\xi }_{2,j_0}^k{\parallel }^2\right)-\sum _{j\in J_1^k\cup T{J}_l}{\lambda }_j^0{\alpha }_{1,j}^k\ge -\eta .$$

(14)

Remark 2. 

In the above case, we embed the switching direction technique. That is, we use $d_l^k$ as a tentative displacement, and if the descent failure occurs, then we implement an Armijo-type line search along the direction $d_{l,j_0}^k$. Only in the case of failure of the latter is a null step declared.

Next, we present the proposed BEM-DC algorithm 1 in the step-by-step format. We denote by $i_{\mathrm{d}s}=1$ when the switching direction technique is embedded.

Algorithm 1: BEM-DC

Require: The stopping tolerance parameter $\eta >0$, the null step parameter $\theta >0$, the proximity threshold ${\delta }_{min}>0$, the sufficient descent parameter $\mu \in (0,1)$, and the step size reduction parameters ${\sigma }_1,{\sigma }_2\in (0,1)$. Ensure: An approximate critical points of Problem (1). Select a starting point $x_1\in R^n.$ Compute ${\xi }_{1,1}\in \partial f_1\left(x_1\right)$ and ${\xi }_{2,1}\in \partial f_2\left(x_1\right)$. Set $B^1=\left\{\left(x_1,f_1\left(x_1\right),f_2\left(x_1\right),{\xi }_{1,1},{\xi }_{2,1},0,0\right)\right\}$, $J_1^1=\left\{1\right\}$, $J_2^1=\left\{1\right\}$, and $k=1$.

Remark 3. 

We update the proximity parameter ${\delta }_k$ using the following formula given in [46]:

$$\begin{array}{cc}\hfill {\overline{\delta }}_{k+1}& ={\delta }_k\left(1-\frac{f(x_k+d_k)-f\left(x_k\right)}{v_k}\right),\hfill \\ \hfill {\delta }_{k+1}& =max\left\{{\overline{\delta }}_{k+1},{\delta }_k/10,{\delta }_{min}\right\}.\hfill \end{array}$$

(16)

Remark 4. 

An outline of the main differences between BEM-DC and PBDC described in [29] is in order.

• In BEM-DC, only one of the two subgradients ${\xi }_{1,j}$ and ${\xi }_{2,j}$, gathered at any iteration j, enters into the calculation of the tentative displacement at each of the successive iterations; the choice is driven by the sign of the linearization error at the current iterate. This is not the case in PBDC;
• In BEM-DC, there exists a temporary bundle whose “birth and death” takes place within the procedure for escaping from the null step and does not increase the bundle size once a descent is achieved, while PBDC adds more information to the bundles at each step and, thus, increases the bundle sizes at every iteration;
• In BEM-DC, a second direction is checked for the descent before the null step is declared (see Remark 2), whereas there is only one possible direction used in each step of PBDC.

5. Termination Property of Algorithm BEM-DC

In this section, we provide the proof of the Algorithm BEM-DC, taking any starting point $x_1\in {\mathbb{R}}^n$ as an input, that returns an approximate critical point $x^{*}$. Assume that the set

$${\mathcal{F}}_1\triangleq \left\{x\in {\mathbb{R}}^n:f\left(x\right)\le f\left(x_1\right)\right\}$$

is bounded and the numbers $L_1$ and $L_2$ are the Lipschitz constants of $f_1$ and $f_2$, respectively, on the set ${\mathcal{F}}_1.$ Note that, whenever the tuple $b\left(y\right)$ is inserted into the temporary bundle $T{B}_l$ associated with the function $f_1$ at the trial point y (see Step 6), then $t\parallel d_l^k\parallel \le \theta$ implies ${\xi }_1^{+}\left(y\right)\in {\partial }_{\epsilon }{f}_1\left(x_k\right)$ for $\epsilon \le 2\theta L_{1\theta }$, where $L_{1\theta }$ is the Lipschitz constant of $f_1$ on the set

$${\mathcal{F}}_{\theta }\triangleq \left\{x\in {\mathbb{R}}^n:dist(x,{\mathcal{F}}_1)\le \theta \right\}.$$

Lemma 1. 

Let $L_{1\theta }$ and $L_2$ be the Lipschitz constants of $f_1$ and $f_2$, respectively, on the set ${\mathcal{F}}_{\theta }.$ Then the following bound holds:

$$\parallel d_l^k\parallel \le \frac{L_{1\theta }+L_2}{{\delta }_{min}}=D.$$

(17)

Proof. 

Throughout the algorithm, we have ${\delta }_k\ge {\delta }_{min}$, and thus, the inequality follows from (10).    □

Remark 5. 

The bound given in (17) is also valid for $\parallel d_{l,j_0}^{*}\parallel$.

Next, we prove that the number of inner iterations of algorithm BEM-DC is finite.

Lemma 2. 

At any given k-th outer iteration, the inner iteration terminates, either fulfilling the sufficient decrease condition (Step 4) or satisfying the stopping condition (Step 2).

Proof. 

Suppose by contradiction that the inner iteration does not terminate. Since (see Remark 5) $\parallel d_{l,j_0}^{*}\parallel$ is bounded and t becomes arbitrarily small, the algorithm cannot loop infinitely many times between Steps 4 and 5. Thus, it loops infinitely many times between Steps 1 and 7, giving rise to an infinite number of null steps.

Observe now that every time a null step occurs and the tuple $b\left(y\right)$ is generated at Step 6, we have

$$\begin{array}{cc}\hfill f(x_k+t_l{d}_l^k)-f\left(x_k\right)& =f_1(x_k+t_l{d}_l^k)-f_2(x_k+t_l{d}_l^k)-\left(f_1\left(x_k\right)-f_2\left(x_k\right)\right)\hfill \\ & \ge t_l\mu v_l^k,\mathrm{for}\mathrm{some}0<t_l<1.\hfill \end{array}$$

Considering

$${\alpha }_{1,l+1}^k=f_1\left(x_k\right)-f_1(x_k+t_l{d}_l^k)+t_l\langle {\xi }_{1,l+1}^k,d_l^k\rangle ,$$

the convexity of $f_2$, and ${\xi }_{2,j_0}^k\in \partial f_2\left(x_k\right)$, it follows

$$\langle {\xi }_{1,l+1}^k-{\xi }_{2,j_0}^k,d_l^k\rangle -{\alpha }_{1,l+1}^k\ge \left(\frac{1}{t_l}-1\right){\alpha }_{1,l+1}^k+\mu v_l^k>\mu v_l^k.$$

(18)

Furthermore, since the sequence ${\left\{w_l^k\right\}}_{l\to \infty }$, is monotonically non-decreasing and bounded from above by zero, it is convergent. Note that the sequence ${\left\{w_{l,j_0}^k\right\}}_{l\to \infty }$ (see (13)) is convergent as well. Consequently, from boundedness of $\parallel d_{l,j_0}^k\parallel$, there exists a convergent subsequence ${\left\{d_{l,j_0}^k\right\}}_{l\in \mathcal{L}}\to {\overline{d}}_0^k$, and thus, the subsequence ${\left\{v_{l,j_0}^k\right\}}_{l\in \mathcal{L}}$ also converges to a limit, say ${\overline{v}}_0^k\le 0$, for some subset of indices $\mathcal{L}\subset$$\{1,2,\dots \}.$

Next, consider two successive indices $p,q\in \mathcal{L}$. From (18) and the definition of Problem $P_{l,j_0}^k$, we obtain

$$\begin{array}{cc}\hfill \langle {\xi }_{1,p+1}^k-{\xi }_{2,j_0}^k,d_p^k\rangle -{\alpha }_{1,p+1}^k& >\mu v_{p,j_0}^k,\mathrm{and}\hfill \\ \hfill \langle {\xi }_{1,p+1}^k-{\xi }_{2,j_0}^k,d_q^k\rangle -{\alpha }_{1,p+1}^k& \le v_{q,j_0}^k.\hfill \end{array}$$

Therefore, we obtain

$$v_{q,j_0}^k-\mu v_{q,j_0}^k>\langle {\xi }_{1,p+1}^k-{\xi }_{2,j_0}^k,d_q^k-d_p^k\rangle .$$

Passing this to the limit, we have $(1-\mu ){\overline{v}}_0^k\ge 0,$ which, taking into account ${\overline{v}}_0^k\le 0$, implies ${\overline{v}}_0^k=0$. This contradicts that the stopping condition at Step 2 is not satisfied infinitely many times.    □

The next theorem proves the termination of Algorithm BEM-DC.

Theorem 1. 

For any starting point $x_1\in {\mathbb{R}}^n$, the algorithm terminates after finitely many iterations at a point satisfying the stopping criterion at Step 2.

Proof. 

Assume the contrary. That is, the stopping criterion at Step 2 is never fulfilled. This implies, taking into account Lemma 2, that an infinite sequence of serious steps takes place, giving rise to infinitely many outer iterations.

Note that every time a serious step is achieved, considering the failed stopping test at Step 2, we have $v_l^k\le -\eta <0.$ Furthermore, it holds that either $t=1$ or $t\parallel d^k\parallel >\theta$. Consequently, it follows from (15) that in the case of $t=1$, we obtain

$$f(x_k+t{d}^k)-f\left(x_k\right)<-\mu \eta ,$$

(19)

or in the second case, considering (17), we have

$$f(x_k+t{d}^k)-f\left(x_k\right)<-\frac{\mu \eta \theta }{D}.$$

This together with (19) implies that the decrease of the objective function value in both cases is bounded away from zero every time a serious step occurs. This is a contradiction under the assumption that ${\mathcal{F}}_1$ is bounded.    □

Remark 6. 

The “escape procedure” (Algorithm 1, introduced in [33]) can escape from critical points which are not Clarke stationary. Algorithm BEM-DC can be combined with this procedure to design an algorithm for finding Clarke stationary points of Problem (1). More precisely, once an approximate critical point is obtained byBEM-DC, the escape procedure could be applied to approximate the Clarke subdifferential at this point [42,43]. Then, it verifies whether this approximation contains the origin with respect to some tolerance. During this approximation, the procedure generates directions that are used to find elements of the Clarke subdifferential. It is proved that this procedure is finitely convergent. That is, after a finite number of steps, it either confirms that the current critical point is also Clarke stationary or the descent direction is found to escape this point (see [33] for more details about different stationary points and their relationships).

6. Numerical Experiments

To evaluate the performance of the BEM-DC algorithm and to compare it with some existing nonsmooth DC programming algorithms, we carry out numerical experiments using different known academic test problems designed for DC programming. Problems 1–9 are from [29], Problems 11–14 are described in [33], and Problems 15–17 are designed in [36]. We exclude Problem 10 from [29] in our experiments as it has many local minimizers and its usage does not provide an unbiased picture of the performance of local search solvers. For the formulations, optimal values, and initial points of these test problems, we refer to references [29,33,36].

6.1. Solvers and Parameters

We consider the following nonsmooth DC optimization solvers in our comparison:

• Augmented subgradient method for nonsmooth DC optimization (ASM-DC) [36];
• Aggregate subgradient method for nonsmooth DC optimization (AggSub) [35];
• Proximal bundle method for nonsmooth DC optimization (PBDC) [29];
• DC Algorithm (DCA) [8,20];
• Proximal bundle method for nonsmooth DC programming (PBMDC) [30];
• Sequential DC programming with cutting-plane models (SDCA-LP) [25].

The parameters in algorithm BEM-DC are chosen as follows: $\eta ={10}^{-7}$, ${\delta }_{min}={10}^{-5}$, $\theta =0.5$, $\mu =0.2$, ${\sigma }_1=0.2,$ and ${\sigma }_2=0.4.$ The same set of parameters is used with all the test problems. The initial value of the proximity parameter is computed according to the recommendation given in [46]; that is, ${\delta }_1=\parallel {\xi }_{1,1}-{\xi }_{2,1}\parallel .$ We select the parameters of other algorithms considering the recommendations given in their references.

The algorithms BEM-DC, ASM-DC, AggSub, PBDC, and DCA are implemented in Fortran 95 and compiled using the gfortran compiler. For PBMDC and SDCA-LP, we use their MATLAB implementations, available at http://www.oliveira.mat.br/solvers (accessed on 1 August 2023). Since the SDCA-LP algorithm requires the feasible set to be compact, to apply this method to unconstrained problems, we define a large n-dimensional box around the starting point and consider only those points generated by the SDCA-LP that belong to this box.

We set a time limit for all algorithms. For solvers in MATLAB, we consider the limit to be three hours, and for those in Fortran, the limit is half an hour. The source code of algorithm BEM-DC is freely available on GitHub: https://github.com/SnTa2019/Nonsmooth-Optimization (accessed on 1 August 2023) and at http://napsu.karmitsa.fi/bemdc (accessed on 1 August 2023). We carry out our numerical experiments on a computer with Intel(R) Core (TM), i7-9750H, CPU @ 2.60 GHz, 32GB (RAM), under Windows 10, 64Bits. Since the solvers PBMDC and SDCA-LP are implemented in MATLAB, we do not report their CPU time in our comparison.

6.2. Evaluation Measures and Notations

We report the numbers of function and subgradient evaluations and the final value of the objective function found by algorithms. We also provide the computational time (in seconds) required by the algorithms implemented in the same platform (Fortran). The following notations are used:

• Prob. is the label of the problem;
• n is the number of variables;
• $N_f$ is the number of function evaluations for the objective function f;
• $N_{\xi }$ is the number of subgradient evaluations defined as the average of subgradient evaluations for the DC components $f_1$ and $f_2$;
• CPU is the computational time in seconds;
• $f_A^{*}$ is the best value of the objective function obtained by algorithm A.

All methods used in this paper are local search methods. In this case, a method is successful if it can find a stationary point of a problem with the given accuracy. For each test problem, a subset of stationary points is known. This subset contains, in particular, stationary points found in methods used in this paper. We say that an algorithm “A” finds a solution with respect to a tolerance ${\beta }_1$ if

$$0\le \frac{f_A^{*}-f_{opt}}{|f_{opt}|+1}\le {\beta }_1,$$

(20)

where $f_{opt}$ is the value of the objective function at one of the known stationary points. Otherwise, we say that the solver fails. We set ${\beta }_1={10}^{-4}.$ In addition, we apply performance profiles, introduced in [47], to analyze the results.

6.3. Results

In this subsection, we report and discuss the results of the numerical experiments. We present the results obtained by algorithms for Problems 1–9 from [29] and Problems 11–14 from [33]. We consider two cases by using two types of starting points. In the first case, for each problem, we use one starting point, given in [29,33]. In this case, using tables, we report the values of the objective functions and the numbers of function and subgradient evaluations. In the second case, we use 20 randomly generated starting points for each problem and present the results by applying performance profiles. We also report the average CPU time required by the BEM-DC algorithm on these problems for a different number of variables using both 1 and 20 randomly generated starting points. Further, we show the results for three large-scale test problems from [36] by reporting the values of objective functions and the numbers of function and subgradient evaluations.

Table 1. Best values of the objective functions obtained by solvers (one starting point).

Prob.	n	BEM-DC	ASM-DC	AggSub	PBDC	DCA	PBMDC	SDCA-LP
1	2	2.00000	2.00000	2.00000	2.00000	2.00000	2.00000	–
2	2	0.00000	0.00000	0.00000	0.00000	1.00000	0.00000	0.00000
3	4	0.00000	0.00001	0.00001	0.00000	0.00000	0.00000	2.00000
4	2	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
4	5	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
4	10	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
4	50	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
4	100	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
4	200	0.00000	–	–	0.00000	–	0.00000	0.00000
5	2	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
5	5	0.00000	0.00000	–	0.00000	0.00000	0.00000	0.00000
5	10	0.00001	0.00000	–	0.00000	0.00000	0.00000	0.00001
5	50	0.00004	0.00001	–	0.00000	0.00011	0.00000	0.00002
5	100	0.00000	0.00001	–	0.00000	0.00031	0.00000	0.00007
5	200	0.00000	0.00001	–	0.00000	0.00003	0.00000	0.00008
6	2	−2.50000	−2.50000	−2.50000	−2.50000	−2.50000	−2.50000	−2.50000
7	2	0.50000	0.50000	0.50000	0.50000	1.00000	–	0.50000
8	3	3.50000	3.50000	3.50000	3.50000	–	3.50000	–
9	4	9.20000	9.20000	9.20000	1.83333	9.20000	1.83333	–
11	3	116.33333	116.33333	116.33377	116.33333	116.33333	116.33330	116.33330
12	2	0.61804	0.61804	0.61804	1.61803	0.61803	–	1.61809
12	5	0.61803	0.61804	0.61804	1.61803	0.61803	–	0.61803
12	10	0.61803	0.61804	0.61804	0.61803	0.61803	–	–
12	50	0.61803	0.61804	0.61804	–	0.61804	–	–
12	100	0.61803	0.61803	0.61804	–	0.61804	–	–
12	200	0.61803	0.61803	0.61804	–	0.61804	–	–
13	10	0.00000	0.00000	0.00000	–	0.00000	–	–
14	2	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
14	5	0.00001	0.00000	–	0.00000	0.00000	0.00001	0.00001
14	10	0.00000	0.00000	–	0.00000	0.00000	0.00001	0.00003
14	50	0.00002	0.00000	–	0.00000	0.00000	0.00002	0.00005
14	100	0.00007	0.00001	–	0.00000	0.00000	–	0.00006
14	200	0.00003	0.00003	–	0.00000	0.00000	–	0.00003

,

Table 2. Number of function and subgradient evaluations required by solvers (one starting point).

Prob.	n	BEM-DC		ASM-DC		AggSub		PBDC		DCA		PBMDC		SDCA-LP
		${\mathit{N}}_{\mathit{f}}$	${\mathit{N}}_{\mathit{\xi }}$	${\mathit{N}}_{\mathit{f}}$	${\mathit{N}}_{\mathit{\xi }}$	${\mathit{N}}_{\mathit{f}}$	${\mathit{N}}_{\mathit{\xi }}$	${\mathit{N}}_{\mathit{f}}$	${\mathit{N}}_{\mathit{\xi }}$	${\mathit{N}}_{\mathit{f}}$	${\mathit{N}}_{\mathit{\xi }}$	${\mathit{N}}_{\mathit{f}}$	${\mathit{N}}_{\mathit{\xi }}$	${\mathit{N}}_{\mathit{f}},{\mathit{N}}_{\mathit{\xi }}$
1	2	105	25	189	52	162	75	22	17	60	41	10	16	–
2	2	175	25	109	41	255	92	18	13	72	54	3	5	7
3	4	87	17	149	53	391	174	23	11	222	181	5	8	16
4	2	75	17	43	19	64	31	6	3	52	28	2	2	6
4	5	40	9	126	44	235	120	13	6	165	124	3	4	13
4	10	70	18	256	88	545	273	16	10	361	309	5	8	26
4	50	1865	874	3430	1124	3206	1597	52	32	2915	2807	25	38	103
4	100	9049	4379	10,878	3539	6824	3405	102	66	6312	6232	50	75	204
4	200	40,713	19,626	–	–	–	–	475	546	–	–	101	151	403
5	2	8	3	231	50	75	31	18	4	37	24	2	3	8
5	5	64	26	118	44	–	–	15	7	474	297	6	8	9
5	10	87	42	335	111	–	–	23	15	54,116	53,287	12	24	34
5	50	526	235	1413	457	–	–	135	124	247,813	246,256	17	24	55
5	100	141	71	1337	438	–	–	47	25	471,198	469,845	20	29	63
5	200	129	65	1885	594	–	–	108	52	46,2821	461,267	18	28	66
6	2	39	10	135	43	105	55	22	14	41	33	17	24	29
7	2	303	48	388	111	285	107	72	47	3368	2471	–	–	49
8	3	125	49	271	75	176	88	60	32	–	–	–	–	–
9	4	9	3	199	63	169	80	86	72	63	56	17	25	–
11	3	146	28	174	61	258	126	10	7	18	17	5	8	14
12	2	130	26	290	79	127	64	20	15	61	60	–	–	37
12	5	202	52	457	122	213	103	58	39	67	65	–	–	696
12	10	393	105	831	257	359	170	1415	1257	76	74	–	–	–
12	50	739	156	3015	771	1411	694	–	–	99	97	–	–	–
12	100	1032	173	6998	2519	1239	907	–	–	117	115	–	–	–
12	200	2757	174	15,302	5090	2850	1397	–	–	338	337	–	–	–
13	10	145	21	26	13	70	38	–	–	20	19	–	–	–
14	2	19	4	50	23	71	34	9	5	56	54	11	16	29
14	5	145	26	206	73	226	100	109	105	63	61	22	31	58
14	10	249	59	357	107	322	155	6169	6151	61	60	34	50	107
14	50	1168	439	672	211	1653	790	991	966	16	14	59	89	280
14	100	2273	708	1779	549	–	–	780	740	18	16	–	–	552
14	200	2591	724	1746	541	–	–	983	961	28	26	–	–	619

,

Table 3. CPU time (in seconds) required by BEM-DC.

n	2	3	4	5	10	50	100
1 starting point	0.000	0.000	0.001	0.012	0.009	0.672	16.328
20 starting points	0.000	0.000	0.002	0.018	0.616	14.613	42.461

,

Table 4. Best values of objective functions obtained by solvers for large-scale problems.

P	n	BEM-DC	ASM-DC	AggSub	PBDC	DCA	PBMDC	SDCA-LP
$P_{L1}$	200	0.00000	0.00006	0.00007	98.56721	153.29050	−	−
$P_{L1}$	500	0.00000	0.00009	0.00007	249.00000	539.99295	−	−
$P_{L1}$	1000	0.00000	0.00032	0.00032	499.00000	510.84581	−	−
$P_{L1}$	1500	0.00000	0.00072	0.00051	749.00000	$1.90755\times {10}^3$	−	−
$P_{L1}$	2000	0.00000	0.00247	0.00083	999.00000	$1.00229\times {10}^3$	−	−
$P_{L2}$	200	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	$3.59539\times {10}^3$
$P_{L2}$	500	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000	$1.08151\times {10}^4$
$P_{L2}$	1000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00001	$2.43987\times {10}^4$
$P_{L2}$	1500	0.00000	0.00000	0.00000	0.00000	0.00000	0.00002	$3.90289\times {10}^4$
$P_{L2}$	2000	0.00000	0.00000	0.00000	0.00000	0.00000	0.00003	$5.43386\times {10}^4$
$P_{L3}$	200	0.00004	0.00000	4.66526	0.00000	0.16124	0.00001	0.00005
$P_{L3}$	500	0.00001	0.00001	7.91275	0.00000	0.08352	0.00002	0.00015
$P_{L3}$	1000	0.00006	0.00037	17.32160	0.00000	0.05060	0.00001	0.00004
$P_{L3}$	1500	0.00002	0.00157	10.50161	0.00000	2.78735	0.00001	0.00011
$P_{L3}$	2000	0.00008	0.00056	12.63033	0.00000	26.17493	0.00002	0.00003

and

Table 5. Number of function and subgradient evaluations for large-scale problems.

Prob.	n	BEM-DC		ASM-DC		AggSub		PBDC		DCA		PBMDC		SDCA-LP
		$N_f$	$N_{\xi }$	$N_f$	$N_{\xi }$	$N_f$	$N_{\xi }$	$N_f$	$N_{\xi }$	$N_f$	$N_{\xi }$	$N_f$	$N_{\xi }$	$N_f,N_{\xi }$
$P_{L1}$	200	49	3	2328	628	1075	277	60	48	2508	2502	−	−	−
$P_{L1}$	500	94	3	2083	456	1400	364	146	115	2006	2002	−	−	−
$P_{L1}$	1000	158	3	2165	497	2127	607	50	37	3511	3503	−	−	−
$P_{L1}$	1500	218	4	2691	664	1734	479	102	80	1003	1001	−	−	−
$P_{L1}$	2000	244	4	3123	793	1299	631	49	34	4526	4515	−	−	−
$P_{L2}$	200	863	248	1835	593	1604	787	974	780	5010	5005	108	177	4509
$P_{L2}$	500	1430	424	2161	705	2342	1157	1355	1057	5010	5005	144	23	2264
$P_{L2}$	1000	3589	1057	3153	1041	2277	1125	1664	1416	5904	5896	219	360	1655
$P_{L2}$	1500	4293	1329	3966	1310	2722	1344	2044	1740	4158	4153	295	482	1509
$P_{L2}$	2000	6433	1817	4433	1466	2936	1459	2244	1914	5511	5505	299	497	955
$P_{L3}$	200	1989	536	1044	333	14639	7151	1928	1845	3006	3003	55	97	97
$P_{L3}$	500	4598	729	1444	469	12869	6252	2539	2513	3507	3503	58	99	3555
$P_{L3}$	1000	13,844	3006	3445	1107	21,456	10,555	2980	2958	3006	3003	73	126	370
$P_{L3}$	1500	12,384	2291	2408	730	41,096	20,239	2763	2747	1503	1501	69	121	1130
$P_{L3}$	2000	17,520	3156	4017	1283	34,758	17,149	3766	3741	1503	1501	73	127	1090

report the results, and the sign “–” is used to show the failure of a method in finding the correct solution.

6.3.1. Results with One Starting Point

Table 1 presents the best value of the objective function obtained by solvers using one starting point. We can see that the BEM-DC algorithm is sufficiently accurate in finding local solutions. Except for Problem 9, it finds global minimizers of DC optimization problems in all other cases. Other solvers fail in one or several problems. More specifically, ASM-DC, AggSub, and DCA fail to find solutions in Problem 4 with $n=200$. Furthermore, AggSub fails in Problems 5 and 14 with $n=5,10,50,100,200$; PBDC fails in Problem 12 with $n=50,100,200$ and in Problem 13; DCA fails in Problem 8; PBMDC fails in Problems 7, 12, 13 and also fails in Problem 14 with $n=100,200$; and finally, SDCA-LP fails in Problems 1, 8, 9, 13 and also fails in Problem 12 with $n=10,50,100,200$.

In Table 2, we report the number of function and subgradient evaluations required by the algorithms. These numbers are computed as an average of the number of function and subgradient evaluations of DC components. Note that in SDCA-LP, the number of function and subgradient evaluations are the same. The results show that, in general, PBMDC requires the least number of function and subgradient evaluations among all algorithms. These numbers are similar for PBDC and SDCA-LP. We can see that the computational effort required by the BEM-DC is reasonable and in some instances it is similar to that of by PBMDC, PBDC, and SDCA-LP. Here, Problem 4 with $n=50,100,200$ is an exception. Three other solvers, ASM-DC, AggSub, and DCA also require a reasonable computational effort in most cases; however, in general, they use (in some cases significantly) more function and subgradient evaluations than BEM-DC PBMDC, PBDC, and SDCA-LP.

6.3.2. Results with 20 Starting Points

Performance profiles using the number of function evaluations, the number of subgradient evaluations, and the computational time (CPU time) required by algorithms are depicted in Figure 1, Figure 2 and Figure 3. Here, we use results obtained by solving Problems 1–9 and 11–14 with 20 random starting points. We report the pairwise comparison of performance profiles as per the note by [48]. The comparison with other methods shows that the BEM-DC is more efficient and robust than the AggSub, DCA, PBMDC, and SDCA-LP methods. Compared with the ASM-DC method, we can see that the BEM-DC uses significantly fewer subgradient evaluations, and the performance of these two methods are similar concerning two other measures. The performance of the BEM-DC and PBDC methods are similar for the number of function and subgradient evaluations; however, the BEM-DC is more efficient than PBDC if one uses computational time. One reason for this can be the fact that the BEM-DC better manages the bundle of the second DC component than PBDC and decreases the number of the solving of the quadratic programming subproblems required to find search directions.

For a given number n of variables, we also report the average CPU time (in seconds) required by the BEM-DC for Problems 1–9 and 11–14 (29 problems considering different variables) with 1 starting point and 20 random starting points. The results given in Table 3 indicate that the CPU time required by this algorithm is very small for problems with $n\le 10$ number of variables.

6.3.3. Results for Large-Scale Problems

Table 4 and Table 5 provide the results obtained by different algorithms for three large-scale problems from [36]: $P_{L1},P_{L2}$, and $P_{L3}$. More specifically, Table 4 contains the best objective function values found by each solver, and Table 5 displays the number of function and subgradient evaluations. As we mentioned above, the number of function and subgradient evaluations are the same for SDCA-LP. Only two algorithms, BEM-DC and ASM-DC, can find approximate solutions to all three problems. AggSub fails in Problem $P_{L3}$, PBDC in Problem $P_{L1}$, DCA in Problems $P_{L1}$ and $P_{L3}$, PBMDC in Problem $P_{L1}$, and SDCA-LP in Problems $P_{L1}$ and $P_{L2}$. Summarizing results from Table 4, we can conclude that the BEM-DC is the most accurate among all algorithms for solving large-scale problems used in numerical experiments.

Results presented in Table 5 show that among all algorithms, BEM-DC uses the least number of function and subgradient evaluations for solving Problem $P_{L1}$, and PBMDC requires the least computational effort for solving Problems $P_{L2}$ and $P_{L3}$. Although the BEM-DC can find accurate solutions for Problems $P_{L2}$ and $P_{L3}$, it nevertheless requires more—and in the case of Problem $P_{L3}$, significantly more—function and subgradient evaluations than other solvers. All other algorithms also require reasonable computational effort in problems where they succeeded to find solutions.

7. Conclusions

In this paper, a new method, named the bundle enrichment method (BEM-DC), is introduced for solving nonsmooth unconstrained difference of convex (DC) optimization problems. This method belongs to the family of bundle-type methods. It exploits cutting plane models of DC components to build the model of the DC objective function. The main difference between the proposed method and other bundle methods for DC optimization is in the dynamic management of the bundle. In the implementation of most bundle methods, the size of the bundle for the cutting plane models of the second DC component is given by the user, and it is restricted to reduce the number of subproblems for finding search directions. However, in the BEM-DC, this size is self-determined by the method. This allows for avoiding solving subproblems that do not provide descent directions.

We prove that the BEM-DC computes approximate critical points of the unconstrained DC optimization problems in a finite number of iterations. The performance of this method is evaluated using two groups of nonsmooth DC optimization test problems: small- and medium-sized sized test problems and test problems with a large number of variables. We consider two types of starting points for problems from the first group: the single starting point available in the literature and 20 randomly generated starting points. The use of randomly generated starting points allows us to investigate the robustness of the proposed method. In addition, we provide a comparison of the BEM-DC with six other DC optimization methods.

Results of numerical experiments show that the BEM-DC is able to find accurate solutions using reasonable computational effort in most test problems. Nevertheless, in some large-scale problems, the number of function and subgradient evaluations required by the BEM-DC algorithm may increase significantly as the number of variables increases. This means that the BEM-DC algorithm may not be applicable to some nonsmooth DC optimization problems with a very large number of variables ($n\ge 5000$). The extension of this method for solving such problems will be the subject of future research.

The remarkable feature of the BEM-DC is that it is able to dynamically manage the number of subgradients of the second DC components and significantly decrease the number of quadratic programming subproblems for finding search directions. Results obtained using many randomly generated starting points allow us to conclude that the BEM-DC is efficient and is among the most robust methods used in numerical experiments. Further, the BEM-DC algorithm is able to efficiently solve relatively large nonsmooth DC optimization problems.