The Boosted DC Algorithm for linearly constrained DC programming

TL;DR

This paper extends the Boosted DC Algorithm (BDCA) to linearly constrained DC programs, proving convergence properties and demonstrating improved performance over DCA through numerical experiments on various challenging problems.

Contribution

The paper introduces an extension of BDCA for linearly constrained DC problems, with convergence proofs and empirical performance improvements.

Findings

Cluster points are KKT points under Slater condition.

Sequences are bounded and converge R-linearly for quadratic objectives.

BDCA outperforms DCA in numerical tests on multiple problems.

Abstract

The Boosted Difference of Convex functions Algorithm (BDCA) has been recently introduced to accelerate the performance of the classical Difference of Convex functions Algorithm (DCA). This acceleration is achieved thanks to an extrapolation step from the point computed by DCA via a line search procedure. In this work, we propose an extension of BDCA that can be applied to difference of convex functions programs with linear constraints, and prove that every cluster point of the sequence generated by this algorithm is a Karush--Kuhn-Tucker point of the problem if the feasible set has a Slater point. When the objective function is quadratic, we prove that any sequence generated by the algorithm is bounded and R-linearly (geometrically) convergent. Finally, we present some numerical experiments where we compare the performance of DCA and BDCA on some challenging problems: to test the copositivity of a given matrix, to solve one-norm and infinity-norm trust-region subproblems, and to solve piecewise quadratic problems with box constraints. Our numerical results demonstrate that this new extension of BDCA outperforms DCA.

Equations

$$\left\{\begin{array}[]{ll}\displaystyle\min_{x\in\mathbb{R}^{n}}&\phi(x):=g(x)-h(x)\vspace{1ex}\\ \text{s.t.}&\langle a_{i},x\rangle\leq b_{i},~{}i=1,\ldots,p,\end{array}\right.$$

$$(+\infty )-(+\infty )=+\infty ,$$

$$(+\infty )-\lambda =+\infty \text{and}\lambda -(+\infty )=-\infty ,\forall \lambda \in ]-\infty ,+\infty [.$$

$$\mathcal{F}:=\left\{x\in {\mathbb{R}}^n\mid ⟨a_i,x⟩\le b_i,i=1,\dots ,p\right\}.$$

$$\partial f(x):=\left\{w\in {\mathbb{R}}^n\mid f(y)\ge f(x)+⟨w,y-x⟩,\forall y\in {\mathbb{R}}^n\right\}$$

$$f^{\prime }(x;d):=\underset{t\searrow 0}{\lim }\frac{f(x+td)-f(x)}{t}.$$

$$f^{\ast }(u)=\underset{x\in {\mathbb{R}}^n}{\sup }\{⟨x,u⟩-f(x)\}.$$

$$u\in \partial f(x)\iff f(x)+f^{\ast }(u)=⟨x,u⟩\iff x\in \partial f^{\ast }(u).$$

$$\iota_{C}(x):=\left\{\begin{array}[]{ll}0,&\text{if }x\in C,\\ +\infty,&\text{otherwise;}\end{array}\right.$$

$$N_{C}(x):=\left\{\begin{array}[]{ll}\{u\in\mathbb{R}^{n}:\langle u,c-x\rangle\leq 0,\quad\forall c\in C\},&\text{if }x\in C,\\ \emptyset,&\text{otherwise.}\end{array}\right.$$

$$\underset{x\in \mathcal{F}}{\inf }\varphi (x)>-\infty .$$

$$\underset{u\in {\mathbb{R}}^n}{\inf }{h}^{\ast }(u)-(g+{\iota }_{\mathcal{F}}{)}^{\ast }(u),$$

$$\underset{u\in {\mathbb{R}}^n}{\inf }\{h^{\ast }(u)-(g+{\iota }_{\mathcal{F}}{)}^{\ast }(u)\}$$

$$=\underset{x\in {\mathbb{R}}^n}{\inf }\{g(x)+{\iota }_F(x)-\underset{u\in {\mathbb{R}}^n}{\sup }\{⟨u,x⟩-h^{\ast }(u)\}\}$$

$$=\underset{x\in {\mathbb{R}}^n}{\inf }\{g(x)+{\iota }_F(x)-h(x)\}$$

$$=\underset{x\in \mathcal{F}}{\inf }\{g(x)-h(x)\}=\underset{x\in \mathcal{F}}{\inf }\varphi (x).$$

$$\left\{\begin{array}[]{l}0\in\nabla{g}(\overline{x})-\partial{h}(\overline{x})+\sum_{i=1}^{p}\mu_{i}a_{i},\vspace{1ex}\\ 0=\mu_{i}(\langle a_{i},\overline{x}\rangle-b_{i}),~{}~{}i=1,\ldots,p,\vspace{1ex}\\ \mu_{i}\geq 0,~{}~{}\langle a_{i},\overline{x}\rangle\leq b_{i},~{}~{}i=1,\ldots,p.\end{array}\right.$$

$$D(\overline{x}):=\left\{d\in {\mathbb{R}}^n\mid \exists \epsilon >0\text{such that}\overline{x}+td\in \mathcal{F},\forall t\in \left[0,\epsilon \right]\right\}.$$

$$N_{\mathcal{F}}(\overline{x})=A(\overline{x}):=\mathrm{cone}\left\{a_i,i\in I(\overline{x})\right\},\text{for }\overline{x}\in \mathcal{F},$$

$$I(\overline{x})=\{i\in \{1,\dots ,p\}\mid ⟨a_i,\overline{x}⟩=b_i\}.$$

$$D(\overline{x})=\left\{d\in {\mathbb{R}}^n\mid ⟨a_i,d⟩\le 0,i\in I(\overline{x})\right\}.$$

$$~{}\left\{\begin{array}[]{ll}\displaystyle\min_{x\in\mathbb{R}^{n}}&\phi_{k}(x):=g(x)-\langle u_{k},x\rangle\vspace{1ex}\\ \textrm{s.t.}&\langle a_{i},x\rangle\leq b_{i},~{}i=1,\ldots,p.\end{array}\right.$$

$$\varphi \left(x_k+{\beta }^l{d}_k\right)\le \varphi (x_k)-\alpha {\beta }^l\Vert d_k{\Vert }^2,$$

$${\overline{\lambda }}_k\le {\hat{\lambda }}_k:=\min \left\{\frac{b_i-⟨a_i,y_k⟩}{|⟨a_i,d_k⟩|}:i\notin I(x_k)\text{ with }⟨a_i,d_k⟩\ne 0\right\};$$

$$I(y_k)\subseteq I(x_k)\iff d_k:=y_k-x_k\in D(y_k)\iff d_k\perp a_i,\forall i\in I(y_k).$$

$$⟨a_i,d_k⟩$$

$$\varphi \left(y_k+\lambda d_k\right)\le \varphi (y_k)-\alpha {\lambda }^2\Vert d_k{\Vert }^2,\text{for all }\lambda \in [0,{\delta }_k].$$

$${\varphi }^{\prime }(y_k;d_k)$$

$$=\underset{t\downarrow 0}{\lim }\frac{g(y_k+t{d}_k)-g(y_k)}{t}-\underset{t\downarrow 0}{\lim }\frac{h(y_k+t{d}_k)-h(y_k)}{t}$$

$$\le ⟨\nabla g(y_k),d_k⟩-⟨v,d_k⟩,$$

$$\left\{\begin{array}[]{l}\nabla{g}(y_{k})+\sum_{i=1}^{p}\mu_{k,i}a_{i}=u_{k}\in\partial h(x_{k}),\vspace{1ex}\\ \mu_{k,i}(\langle a_{i},y_{k}\rangle-b_{i})=0,~{}~{}\mu_{k,i}\geq 0,~{}~{}\langle a_{i},y_{k}\rangle\leq b_{i},~{}~{}i=1,\ldots,p.\end{array}\right.$$

$$⟨u_k-v,x_k-y_k⟩\ge \rho \Vert x_k-y_k{\Vert }^{\text{2}}.$$

$$⟨\nabla g(y_k)-v,d_k⟩$$

Full text

Dedicated to Professor Miguel A. Goberna on the occasion of his 70th birthday

Keywords: Difference of convex functions; boosted difference of convex functions algorithm; global convergence; constrained DC programs; copositivity problem; trust region subproblem.

1 Introduction

In this paper, we are interested in solving the following DC (difference of convex functions) optimization problem:

[TABLE]

where $g:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$, $h:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ are proper, closed, and convex functions with $g$ being smooth, $a_{i}\in\mathbb{R}^{n}$, $b_{i}\in\mathbb{R}$ for $i=1,\ldots,p$, and $\langle\cdot,\cdot\rangle$ denotes an inner product. We use the conventions:

[TABLE]

Observe that we can rewrite problem ($\mathcal{P}$) as an unconstrained nonsmooth DC optimization problem, whose objective function is $g+\iota_{\mathcal{F}}-h$, where $\iota_{\mathcal{F}}$ denotes the indicator function of the feasible set

[TABLE]

For solving this problem, one can apply the classical DC Algorithm (DCA) in [28, 15]. DC programming and the DCA have been developed and studied for more than 30 years [15]. The DCA has been successfully applied in different fields such as machine learning, financial optimization, supply chain management, and telecommunication, see, e.g. [17, 11, 24]. Nowadays, DCA has become a useful method to solve nonconvex problems.

To accelerate the convergence of DCA, which can be slow for some problems, a new method called Boosted DC Algorithm (BDCA) has been recently proposed in [1, 3]. The key idea of BDCA is to perform an extrapolation step via a line search procedure at the point computed by DCA at each iteration. This step allows the algorithm to take longer steps than the classical DCA, achieving in this way a larger reduction of the objective value per iteration. In addition to accelerating its convergence, BDCA may have better chances to escape from bad local optima thanks to the line search procedure, see [3, Example 3.3]. Therefore, BDCA is not only faster than DCA but also can provide better solutions. Extensive numerical experiments in diverse applications such as biochemistry [1], machine learning [36, 22] and data science [3] have been performed where BDCA clearly outperforms DCA. Note that all these applications are modeled as DC problems with $g$ smooth. However, it is important to emphasize that, for unconstrained DC programs, the BDCA proposed in [1, 3] is not applicable when the function $g$ in ($\mathcal{P}$) is nonsmooth (see [3, Example 3.4]).

The aim of this paper is to show that BDCA can still be applied if the nonsmooth function $g$ is the sum of a smooth convex function and the indicator function of a polyhedral set. More precisely, we will show that it is possible to use BDCA for solving DC programs with linear constraints of the form ($\mathcal{P}$). To compare the performance of DCA and BDCA, we provide numerical experiments on various challenging problems: to test the copositivity of a given matrix, to solve $\ell_{1}$ and $\ell_{\infty}$ trust-region subproblems, and to find the minimum of a piecewise quadratic function with box constraints. The three first problems have a quadratic objective function and are known to be NP-hard [23] (the first was already heuristically investigated in [8] using DCA), while the last problem has a nonsmooth objective function. Our results confirm that BDCA significantly outperforms DCA in these applications.

The rest of this paper is organized as follows. Section 2 recalls some preliminary results. In Section 3, we propose a new variant of BDCA for solving ($\mathcal{P}$) and prove that the objective value of the sequence generated by the algorithm is monotonically decreasing, and that any limit point of the sequence is a KKT point of the problem. The R-linear convergence of any sequence generated by BDCA in the special case of quadratic objective functions is derived in Section 4. In Section 5, we provide some numerical experiments for testing the copositivity of a given matrix, for solving $\ell_{1}$ and $\ell_{\infty}$ trust-region subproblems and for solving piecewise quadratic problems with box constraints, where we compare BDCA and DCA. Finally, some conclusions and future research are briefly discussed in Section 6.

2 Preliminaries

In this section, we state our assumptions imposed on ($\mathcal{P}$). We also recall some preliminary and basic results which will be used in the sequel.

Let $f:\mathbb{R}^{n}\to\mathbb{R}\cup\left\{+\infty\right\}$ be a proper extended real-valued convex function. The set $\mathrm{dom}\,f:=\left\{x\in\mathbb{R}^{n}\mid f(x)<+\infty\right\}$ denotes its (effective) domain, and

[TABLE]

denotes the subdifferential of $f$ at $x$. If $f$ is differentiable at $x$, then $\partial{f}(x)=\{\nabla f(x)\}$, where $\nabla f(x)$ stands for the gradient of $f$ at $x$. The one-side directional derivative of $f$ at $x$ with respect to the direction $d\in\mathbb{R}^{n}$ is

[TABLE]

Recall that $f$ is said to be strongly convex with strong convexity parameter $\rho>0$ if the function $f-\frac{\rho}{2}\|\cdot\|^{2}$ is convex. The conjugate of $f$ is the function $f^{\ast}:\mathbb{R}^{n}\to\mathbb{R}\cup\left\{+\infty\right\}$ defined at $u\in\mathbb{R}^{n}$ as

[TABLE]

It holds that $f^{\ast\ast}:=(f^{\ast})^{\ast}=f$. Moreover, for all $x,u\in\mathbb{R}^{n}$,

[TABLE]

Given a set $C\subseteq\mathbb{R}^{n}$, its interior, its relative interior and the convex cone generated by $C$ are denoted by $\mathrm{int}\,C$, $\mathrm{ri}\,C$ and $\mathop{\rm cone}C$, respectively. The function

[TABLE]

denotes the indicator function of $C$, and $C$ is convex if and only if $\iota_{C}$ is so. Moreover, its subdifferential $\partial\iota_{C}$ is the normal cone to $C$, which is given by

[TABLE]

For proving our convergence results, we will make use of the following assumptions.

Assumption 1*.*

Both $g$ and $h$ are strongly convex on their domain with the same strong convexity parameter $\rho>0$.

Assumption 2*.*

The function $h$ is subdifferentiable at every point in $\mathrm{dom}\,h$; i.e., $\partial h(x)\neq\emptyset$ for all $x\in\mathrm{dom}\,h$. The function $g$ is continuously differentiable on an open set containing $\mathrm{dom}\,h$ and

[TABLE]

Assumption 3*.*

The feasible set $\mathcal{F}$ has a Slater point, that is, there exists $\hat{x}\in\mathbb{R}^{n}$ such that $\langle a_{i},\hat{x}\rangle<b_{i}$, for all $i=1,\ldots,p$.

Remark 2.1*.*

Assumption 1 is not restrictive, in the sense that any DC decomposition of $\phi$ as $\phi=g-h$, can be expressed as $\phi=(g+\frac{\rho}{2}\|\cdot\|^{2})-(h+\frac{\rho}{2}\|\cdot\|^{2})$ for any $\rho>0$. Observe that $\partial h(x)\neq\emptyset$ holds for all $x\in\mathrm{ri}\,\mathrm{dom}\,h$ (by [34, Theorem 23.4]), so the first part of Assumption 2 is clearly satisfied if $\mathrm{dom}\,h=\mathbb{R}^{n}$. A key point of our method is the smoothness of $g$ in Assumption 2, which cannot be in general omitted (see [3, Example 3.4]).

Associated with problem ($\mathcal{P}$), we can construct its dual as

[TABLE]

which is also a DC program with the same optimal value. Indeed,

[TABLE]

We say that $\overline{x}\in\mathbb{R}^{n}$ (resp. $\overline{u}\in\mathbb{R}^{n}$) is a *critical point *of ($\mathcal{P}$) (resp. ($\mathcal{D}$)) if $\nabla g(\overline{x})\in\partial(h+\iota_{\mathcal{F}})(\overline{x})$ (resp. $\partial h^{\ast}(\overline{u})\cap\partial(g+\iota_{\mathcal{F}})^{\ast}(\overline{u})\neq\emptyset$). Recall from [21, Theorem 5.19] that $\overline{x}$ is called a KKT point of ($\mathcal{P}$) if there exist $\mu_{1},\mu_{2},\ldots,\mu_{p}\in\mathbb{R}$ such that

[TABLE]

The DC algorithm is a primal-dual method in the sense that it is aimed to find solutions for both ($\mathcal{P}$) and ($\mathcal{D}$), simultaneously. Our goal then is to design a BDCA variant that allow us to find KKT points of ($\mathcal{P}$) and critical points of ($\mathcal{D}$).

To finish this section, we need to introduce the following notions regarding the geometry of the feasible set $\mathcal{F}$. The cone of feasible directions at $\overline{x}\in\mathcal{F}$ is denoted by

[TABLE]

Due to the polyhedral structure of $\mathcal{F}$, its normal cone coincides with the active cone; i.e.,

[TABLE]

where $I(\overline{x})$ stands for the set of active constraints at $\overline{x}$, i.e.,

[TABLE]

Since we deal with affine constraints, we have (see e.g. [2, Proposition 4.14])

[TABLE]

3 The Boosted DC Algorithm and its convergence

For solving ($\mathcal{P}$), we propose the following method, Algorithm 1, which can tackle more general problems than the Boosted DC Algorithm proposed in [3].

Let us make some comments on Algorithm 1.

• (i)

Lines 1 to 6 of Algorithm 1 correspond to the classical DCA for solving ($\mathcal{P}$). Thus, when $\overline{\lambda}_{k}=0$ for all $k$, Algorithm 1 coincides with DCA.

• (ii)

Lines 7 to 15 present the boosting step. It first checks if $d_{k}$ is a feasible direction at $y_{k}\in\mathcal{F}$. If so, it then performs a line search step along the direction $d_{k}$ which maintains feasibility to improve the objective value $\phi$. Otherwise, the boosting step is skipped and we simply use the DCA point $y_{k}$.

• (iii)

In terms of per-iteration complexity, the boosting step requires to check the feasibility of direction $d_{k}$, which can be done by comparing the sets of active constraints at $x_{k}$ and $y_{k}$ (see Lemma 3.1). It also requires evaluating the objective function and checking the feasibility of the trial step $y_{k}+\lambda_{k}d_{k}$. The computational effort of this task will depend on the particular structure of $\phi$ and $\mathcal{F}$. In the case of box constraints, the largest step-size that makes $y_{k}+\lambda_{k}d_{k}$ feasible can be efficiently computed in Line 8, see Remark 3.1 below. In Section 5 we show some computational results for the trust-region subproblems with $\ell_{1}$ norm constraints (where feasibility needs to be checked) and $\ell_{\infty}$ norm constraints (where the largest step-size can be readily computed).

• (iv)

When $h$ is differentiable, the algorithms introduced by Fukushima and Mine in [10, 20] can be applied in our setting. On the one hand, the algorithm in [10] performs a line search which is similar to the one in Lines 9-11 of Algorithm 1, but with the significant difference that it is performed at the point $x_{k}$, instead of doing it at the DCA point $y_{k}$; that is, it searches for the smallest non-negative integer $l$ such that

[TABLE]

where $0<\beta<1$. Thus, the largest step-size allowed by their algorithm is $1$, which corresponds with the point determined by the DCA, since $x_{k}+d_{k}=y_{k}$. On the other hand, the algorithm defined in [20] performs an exact line search in the direction $d_{k}$, which may be unaffordable in many practical applications.

Remark 3.1*.*

An explicit formula to compute the largest step-size that makes $y_{k}+\lambda_{k}d_{k}$ feasible in Line 8 is provided in [9]. More precisely, feasibility of the trial step size is guaranteed if $\overline{\lambda}_{k}$ is chosen so that

[TABLE]

see [9, Lemma 5.4(ii)]. In principle, this permits to avoid the extra time consumed for checking feasibility in Line 8. However, the computation of $\widehat{\lambda}_{k}$ in (4) may be even more expensive and inefficient than checking feasibility of $\lambda_{k}$ in some applications. This is the case, for instance, of the $\ell_{1}$-norm trust region subproblem, whose reformulation as a linearly constrained DC program requires an exponential number of constraints (see Section 5.2).

The next auxiliary lemma shows the equivalence between Line 7 of Algorithm 1 and checking the feasibility of the direction generated by DCA.

Lemma 3.1**.**

If $x_{k}$ and $y_{k}$ are generated by Algorithm 1, then

[TABLE]

Proof.

Observe that, for any $i\in I(y_{k})$, it holds that

[TABLE]

Hence, the result easily follows by taking into account (3). ∎

In the following proposition, we collect some key inequalities which are useful in the sequel for the convergence analysis of Algorithm 1.

Proposition 3.1**.**

Under Assumptions 1 and 2, for all $k\in\mathbb{N}$, the next statements hold:

1. (i)

$\phi(y_{k})\leq\phi(x_{k})-\rho\|d_{k}\|^{2}$; 2. (ii)

$\phi^{\prime}(y_{k};d_{k})\leq-\rho\|d_{k}\|^{2}$; 3. (iii)

if the condition at Line 7 of Algorithm 1 holds, then there exists some $\delta_{k}>0$ such that $y_{k}+\lambda_{k}d_{k}\in\mathcal{F}$ and

[TABLE]

Consequently, the backtracking step at Lines 9–11 of Algorithm 1 terminates after a finite number of iterations.

Proof.

The proof of (i) is similar to the one of [1, Proposition 3] and is therefore omitted. To prove (ii), pick any $v\in\partial h(y_{k})$. Note that the one-sided directional derivative $\phi^{\prime}(y_{k};d_{k})$ is given by

[TABLE]

by convexity of $h$. Since $y_{k}$ is the unique solution of the strongly convex problem ($\mathcal{P}_{k}$), we can write down the KKT conditions (see, e.g., [2, Theorem 4.20]) of this problem as

[TABLE]

The fact that $h$ is strongly convex with a parameter $\rho$ implies, by [35, Exercise 12.59], that $\partial h$ is strongly monotone with constant $\rho$. Therefore, since $v\in\partial h(y_{k})$ and $u_{k}\in\partial h(x_{k})$, we have

[TABLE]

Hence, combining these expressions, together with the fact that $x_{k}\in\mathcal{F}$, we can derive

[TABLE]

and the result follows by combining the last inequality with (3).

Having in mind the condition at Line 7 of Algorithm 1 and Lemma 3.1, we observe that the proof of (iii) is similar to the one of [3, Proposition 3.1], so we omit it for brevity. ∎

Remark 3.2* (General convex constraints).*

Consider a generalized version of ($\mathcal{P}$) where the feasible set is formed by arbitrary convex constraints, i.e.,

[TABLE]

where $g$ and $h$ satisfy Assumptions 1 and 2 and $c_{i}:\mathbb{R}^{n}\to\mathbb{R}$ are smooth, proper, closed and convex functions, for $i=1,\ldots,p.$ Note that problem ($\mathcal{P}$) is a particular instance of ($\mathcal{P}^{\prime}$) with $c_{i}(x):=\langle a_{i},x\rangle-b_{i}$, for $i=1,\ldots,p.$ The assertion in Proposition 3.1(ii) still holds true for the more general problem ($\mathcal{P}^{\prime}$); that is, the direction generated by DCA remains a descent direction provided that $x_{k}$ is feasible for ($\mathcal{P}^{\prime}$). To confirm this, one can easily check that the proof can be rewritten by replacing the linearity of the gradients by the inequality

[TABLE]

However, Line 7 of Algorithm 1 is no longer useful to verify if $d_{k}$ is a feasible direction, as the equality in (3) only holds for affine constraints. For general convex constraints, we have the inclusion

[TABLE]

Therefore, one possibility would be to run the boosting step whenever $\langle\nabla c_{i}(y_{k}),d_{k}\rangle<0$ for all $i\in I(y_{k})$. Nevertheless, this will never be the case because $x_{k}$ is feasible for ($\mathcal{P}^{\prime}$). Indeed, from (7), we obtain that

[TABLE]

In fact, it can be proved that if $y_{k}+\lambda d_{k}\in\mathcal{F}$ for some particular $\lambda>0$, then the points in the segment $[x_{k},y_{k}+\lambda d_{k}]$ must be active for all $i\in I(y_{k})$.

We are now in the position to establish the main convergence result of Algorithm 1.

Theorem 3.1**.**

Under Assumptions 1, 2 and 3, for any $x_{0}\in\mathcal{F}$, either BDCA returns a KKT point of ($\mathcal{P}$) or it generates an infinite sequence such that the following statements hold.

1. (i)

$\phi(x_{k})$* is monotonically decreasing and hence convergent to some $\bar{\phi}$.* 2. (ii)

Suppose that $\{x_{k}\}$ and $\{u_{k}\}$ are bounded. Then any limit point of $\left\{x_{k}\right\}$ is a KKT point of ($\mathcal{P}$) and any limit point of $\{u_{k}\}$ is a critical point of ($\mathcal{D}$). 3. (iii)

We have $\sum_{k=0}^{+\infty}\|d_{k}\|^{2}<+\infty$. Moreover, if there is some $\overline{\lambda}$ such that $\lambda_{k}\leq\overline{\lambda}$ for all $k$, then $\sum_{k=0}^{+\infty}\|x_{k+1}-x_{k}\|^{2}<+\infty$.

Proof.

If Algorithm 1 is terminated at Line 5 and returns $x_{k}$, then $x_{k}=y_{k}$. From (2) and (6), it is clear that $x_{k}$ is a KKT point of ($\mathcal{P}$). Otherwise, by Proposition 3.1 and Line 15 of Algorithm 1, we have

[TABLE]

where $\lambda_{k}\geq 0$. Therefore, the sequence $\left\{\phi(x_{k})\right\}$ converges to some $\bar{\phi}$, since it is monotonically decreasing and bounded from below, by (1). As a consequence, we obtain

[TABLE]

which implies $\|d_{k}\|^{2}=\|y_{k}-x_{k}\|^{2}\to 0$, by (8).

Now, if $\overline{x}$ is a limit point of $\left\{x_{k}\right\}$, then there exists a subsequence $\left\{x_{k_{j}}\right\}$ converging to $\overline{x}$. Then, as $\|y_{k_{j}}-x_{k_{j}}\|\to 0$, we have $y_{k_{j}}\to\overline{x}$. From (6), we obtain

[TABLE]

Since the sequence $\{u_{k}\}$ is bounded by assumption, without lost of generality we may assume that $u_{k_{j}}\to\overline{u}$. It also holds that $\nabla{g}(y_{k_{j}})\to\nabla{g}(\overline{x})$ by continuity of $\nabla g$. The sequence of Lagrange multipliers $\{\mu_{k_{j}}\}\in\mathbb{R}^{p}$ must also be bounded. Indeed, suppose to the contrary that $\|\mu_{k_{j}}\|\to\infty$ and assume without loss of generality that $\lim_{j\to\infty}\frac{\mu_{k_{j}}}{\|\mu_{k_{j}}\|}=\mu^{*}$, with $\mu^{*}\in\mathbb{R}^{p}_{+}$ (the non-negative orthant) satisfying $\|\mu^{*}\|=1$. Then, dividing the first equality in (9) by $\|\mu_{k_{j}}\|$ and letting $j\to\infty$, we obtain

[TABLE]

Performing the same procedure in the second equality in (9) gives

[TABLE]

so we deduce that $\mu_{i}^{*}=0$ for all $i\not\in I(\overline{x})$. Thus,

[TABLE]

Thanks to the Slater Assumption 3, we know that there exists $\hat{x}\in\mathbb{R}^{n}$ such that $\langle a_{i},\hat{x}\rangle<b_{i}$ for all $i\in I(\overline{x})$. Hence,

[TABLE]

which implies $\mu_{i}^{*}=0$ for all $i\in I(\overline{x})$, since $\mu^{*}\in\mathbb{R}^{p}_{+}$. This implies that $\mu^{*}=0_{p}$, a contradiction with the fact that $\|\mu^{*}\|=1$.

Therefore, by extracting subsequences if necessary, we can assume that

[TABLE]

Taking the limit as $j\to\infty$ in (9), thanks to the closedness of the graph of $\partial{h}$ (see [34, Theorem 24.4]), we obtain

[TABLE]

which means that $\overline{x}$ is a KKT point of ($\mathcal{P}$). From (11) we derive that $\overline{x}\in\partial h^{\ast}(\overline{u})$ and also that

[TABLE]

which is equivalent to $\overline{x}\in\partial(g+\iota_{\mathcal{F}})^{\ast}(\overline{u})$ and, hence, $\overline{u}$ is a critical point of ($\mathcal{D}$). The proof of (iii) is similar to that of [1, Proposition 5(iii)] and is thus omitted. ∎

Remark 3.3*.*

The assertion in Theorem 3.1(ii) remains valid for any limit point of $\{x_{k}\}$ in the interior of $\mathrm{dom}\,h$ without requiring the dual sequence $\{u_{k}\}$ to be bounded. Indeed, let $\overline{x}\in\mathrm{int}\,\mathrm{dom}\,(h)$ be a cluster point of $\{x_{k}\}$ and let $\{x_{k_{j}}\}$ be a subsequence converging to $\overline{x}$. We can apply [34, Corollary 24.5.1] to obtain that for any $\varepsilon>0$, there exists a positive integer $j_{0}\geq 1$ such that

[TABLE]

where $\mathbb{B}$ denotes the Euclidean unit ball of $\mathbb{R}^{n}$. Hence, since $\partial h(\overline{x})$ is a bounded set, the dual subsequence $\{u_{k_{j}}\}$ is also bounded and the proof of Theorem 3.1(ii) stands.

Remark 3.4*.*

Similar to [1, Theorem 1] and [3, Theorem 4.3 and Theorem 4.9], if we further assume that the function $\phi$ satisfies the Kurdyka–Łojasiewicz property, then it can be proved that the sequence $\left\{x_{k}\right\}$ converges to a KKT point of ($\mathcal{P}$). Moreover, convergence rates can also be deduced depending on the Łojasiewicz exponent. Especially, when the objective function $\phi$ is quadratic (e.g., in some of our numerical experiments), it was proved [18, Theorem 4.2] that the function $\phi+\iota_{\mathcal{F}}$ satisfies the Kurdyka–Łojasiewicz property with exponent $\frac{1}{2}$. Combining this with the technique in [1, Theorem 1], it is a routine task to derive the linear convergence of the sequence $\left\{x_{k}\right\}$ when the sequence has a cluster point. The purpose of the next section is to prove, using a similar technique to [32, Theorem 2.1], that the latter condition is not needed: for quadratic functions, BDCA always generates a sequence which is linearly convergent to a KKT point of the problem without requiring the existence of a cluster point.

4 Linear convergence for quadratic objective functions

In this section we prove the R-linear convergence of BDCA when the objective function $\phi$ of ($\mathcal{P}$) is quadratic, that is, for problems of the form

[TABLE]

where $Q\in\mathbb{R}^{n\times n}$ is a symmetric matrix, $q\in\mathbb{R}^{n}$, and $a_{i}\in\mathbb{R}^{n}$ and $b_{i}\in\mathbb{R}$, for $i=1,\ldots,p$. In this setting, $\overline{x}$ is a KKT point of ($\mathcal{P}_{Q}$) if there exist multipliers $\mu_{1},\mu_{2},\ldots,\mu_{p}\in\mathbb{R}$ such that

[TABLE]

We denote by $\overline{\mathcal{F}}$ the set of all KKT points of ($\mathcal{P}$).

As $Q$ is not required to be positive semidefinite, ($\mathcal{P}_{Q}$) is a nonconvex problem. However, the matrix $Q$ can be easily decomposed as $Q=Q_{1}-Q_{2}$, with $Q_{1}$ and $Q_{2}$ positive definite. Indeed, one can take

[TABLE]

for $\sigma>\max\{0,\lambda_{\max}(Q)\}$, where $\lambda_{\max}(Q)$ is the largest eigenvalue of $Q$ and $I$ denotes the identity matrix. Thus, ($\mathcal{P}_{Q}$) can be equivalently written in the form of ($\mathcal{P}$) as

[TABLE]

with

[TABLE]

Observe that both functions $g$ and $h$ are strongly convex with parameters $\sigma$ and $\sigma-\lambda_{\max}(Q)$, respectively. Thus, Assumption 1 holds for

[TABLE]

while Assumption 2 trivially holds.

By using the indicator function $\iota_{\mathcal{F}}$ of the feasible set $\mathcal{F}$, we can rewrite problem ($\mathcal{P}_{Q}$) as an unconstrained nonsmooth DC optimization problem, which can be tackled by the DCA. In this case, the DCA becomes the projected gradient method from convex programming [30]

[TABLE]

where $P_{\mathcal{F}}$ denotes the projection mapping onto the feasible set $\mathcal{F}$.

To derive the R-linear convergence of the iterative sequence generated by the BDCA, we will use the following lemmas. The first one is a classical result regarding the connected components of the KKT set $\overline{\mathcal{F}}$ and can be found in [19, Lemma 3.1].

Lemma 4.1**.**

Let $\mathcal{F}_{1},\mathcal{F}_{2},\ldots,\mathcal{F}_{r}$ be the connected components of the KKT set $\overline{\mathcal{F}}$. Then we have

[TABLE]

and the following properties hold:

1. (i)

each $\mathcal{F}_{i}$ is the union of finitely many polyhedral convex sets; 2. (ii)

the sets $\mathcal{F}_{i}$, $i=1,2,\ldots,r$, are properly separated from each others, that is

[TABLE] 3. (iii)

$\phi$ is constant on each $\mathcal{F}_{i}$.

The second one is a local error bound result originally stated in [19, Theorem 2.3] and later extended in [32, Lemma 2.1].

Lemma 4.2**.**

There exist scalars $\epsilon>0$ and $\tau>0$ such that

[TABLE]

for all $x\in\mathcal{F}$ with

[TABLE]

We are now in a position to establish the R-linear (aka geometric) convergence of the iterative sequence $\{x_{k}\}$ generated by BDCA. Geometric convergence is a special type of R-linear convergence, which is implied by Q-linear convergence (see, e.g. [2, Section 5.2.1]). The next result extends [32, Theorem 2.1] to the BDCA. Its proof employs similar techniques, which are originally based on [19].

Theorem 4.1**.**

If ($\mathcal{P}_{Q}$) has a solution and Assumption 3 holds, then the sequence $\{x_{k}\}$ generated by BDCA converges geometrically to a KKT point $\overline{x}$ of ($\mathcal{P}_{Q}$), that is, there exist some constants $C>0$ and $\eta\in{]0,1[}$ such that

[TABLE]

Proof.

First, observe that by (8) we have

[TABLE]

By Theorem 3.1(i), we know that the right-hand side of this inequality converges to zero as $k\to\infty$. Thus,

[TABLE]

which implies

[TABLE]

Let $\epsilon>0$ and $\tau>0$ be such that (16) and (17) hold. Since $y_{k}=P_{\mathcal{F}}\left(x_{k}-\frac{1}{\sigma}(Qx_{k}+q)\right)$, there exists some $k_{0}\in\mathbb{N}$ such that

[TABLE]

Hence, we obtain

[TABLE]

Due to the fact that $\overline{\mathcal{F}}$ is nonempty and closed, there exists $z_{k}\in\overline{\mathcal{F}}$, for each $k\in\mathbb{N}$, such that ${d}(x_{k},\overline{\mathcal{F}})=\|x_{k}-z_{k}\|$. Thus

[TABLE]

which implies

[TABLE]

Since

[TABLE]

it follows from (19) and (21) that

[TABLE]

Now, let $\mathcal{F}_{1},\mathcal{F}_{2},\ldots,\mathcal{F}_{r}$ be the connected components of $\overline{\mathcal{F}}$. By Lemma 4.1(ii) and (22) there exists $\mathcal{F}_{0}\in\{\mathcal{F}_{1},\mathcal{F}_{2},\ldots,\mathcal{F}_{r}\}$ and $k_{1}\geq k_{0}$ such that $z_{k}\in\mathcal{F}_{0}$ for all $k\geq k_{1}$. The last assertion of Lemma 4.1 implies that

[TABLE]

Note that, since $z_{k}$ is a KKT point of ($\mathcal{P}_{Q}$), we have $\left\langle Qz_{k}+q,x_{k}-z_{k}\right\rangle\geq 0$. Then,

[TABLE]

From Theorem 3.1(i) we know that $\lim_{k\to\infty}\phi(x_{k})=\bar{\phi}$. Hence, for all $k\geq k_{1}$,

[TABLE]

We prove now that $\bar{\phi}\leq c$. Indeed, on the one hand, from (8) and (23), we have for all $k\geq k_{1}$ that

[TABLE]

On the other hand, since $y_{k}=P_{\mathcal{F}}\left(x_{k}-\frac{1}{\sigma}(Qx_{k}+q)\right)$ and $z_{k}\in\mathcal{F}$, we deduce (see, e.g., [4, Theorem 3.16]) that

[TABLE]

Therefore,

[TABLE]

Combining (25) and (26), we obtain for all $k\geq k_{1}$ that

[TABLE]

From (20), we have

[TABLE]

Hence, we deduce from (27) that

[TABLE]

where $\beta:=(1+\tau)\left(\sigma+\|Q\|\tau+\frac{1+\tau}{2}\|Q\|\right)$. Passing (28) to the limit as $k\to\infty$, we get

[TABLE]

The latter together with (24) imply that $\bar{\phi}=c$. Therefore, it follows from (28) and (8) that

[TABLE]

or, equivalently,

[TABLE]

Since $\left\{\phi(x_{k})\right\}$ is monotonically decreasing to $\bar{\phi}$, we deduce from the last inequality that

[TABLE]

where $M_{0}:=\left(\phi(x_{k_{1}})-\bar{\phi}\right)(\rho+\beta)^{k_{1}}\beta^{-k_{1}}$ and $\eta:=\sqrt{\frac{\beta}{\rho+\beta}}\in{]0,1[}$. Consider now

[TABLE]

Hence, from (8) and (29), we obtain

[TABLE]

which implies

[TABLE]

with $M:=\sqrt{M_{1}M_{0}}$. Then, for all $k\geq k_{1}$ and $m\geq 1$, we have

[TABLE]

This implies that $\left\{x_{k}\right\}$ is a Cauchy sequence and hence converges to a point $\overline{x}\in\mathcal{F}$. According to Theorem 3.1(ii) and Remark 3.3, combined with the fact that $\mathrm{dom}\,h=\mathbb{R}^{n}$ is an open set, we must have that $\overline{x}$ is a KKT point of ($\mathcal{P}_{Q}$). Moreover, passing the inequality (30) to the limit as $m\to\infty$, we obtain

[TABLE]

which concludes the proof. ∎

5 Numerical experiments

The purpose of this section is to compare the performance of BDCA (Algorithm 1) against the classical DCA. We refrain from comparing the algorithms against other methods, as our only aim here is to show that it is more advantageous to apply BDCA than DCA whenever this is possible, not to show that these methods are the most efficient for the problems that we consider. Consequently, we tested both algorithms in the three different settings that can occur, depending on whether the feasible set is unbounded, bounded with general constraints, and bounded with box constraints (so the feasibility in the line search step can be ensured). We applied both algorithms for testing copositivity [8] (unbounded) and for solving the $\ell_{1}$ and $\ell_{\infty}$ trust-region subproblem [7, 12] (bounded, with box constraints in the case of $\ell_{\infty}$). In our last experiment we test the performance on piecewise quadratic problems with box constraints, whose objective function is thus nonsmooth. All these problems clearly satisfy Assumptions 1-3 and admit a DC decomposition with $g$ being a quadratic function. Thus, subproblem $(\mathcal{P}_{k})$ can be solved by computing a projection onto the feasible set $\mathcal{F}$, analogous to that stated in equation (15). The projector onto the feasible set in our applications can be directly computed.

All the codes were written in Python 3.7 and the tests were run on an Intel Core i7-4770 CPU 3.40GHz with 32GB RAM, under Windows 10 (64-bit).

5.1 Testing copositivity

Recall that a given $n\times n$ matrix $A$ is said to be copositive if

[TABLE]

where $\mathbb{R}^{n}_{+}$ stands for the non-negative orthant. Copositivity has recently attracted considerable attention in mathematical optimization, see e.g. [5, 6, 8, 25]. The problem of determining whether a given matrix is not copositive is known to be NP-complete [23] and it can be recast as the following non-convex optimization problem

[TABLE]

The copositivity of $A$ is now equivalent to $\min_{x\in\mathbb{R}^{n}_{+}}\,\phi(x)=0$. In [8], the authors reformulated (31) as a DC problem according to the decomposition in (12)-(13), and applied DCA as a heuristic for testing whether a matrix is not copositive. To be more specific, given $\sigma>\max\left\{\lambda_{\max}(A),0\right\}$, the reformulation of problem (31) as a DC problem becomes

[TABLE]

with

[TABLE]

Under this decomposition, DCA is applied as a heuristic to determine the copositivity of a given matrix as follows: if at some iterate $\phi(x_{k})<0$, then the matrix is non-copositive; otherwise, if a critical point is reached, the instance is undecidable.

Copositive matrices play an important role in graph theory. The size of the largest complete subgraph contained in a given graph $G$, denoted by $\gamma(G)$, is known as the clique number of $G$. If $A$ and $E$ are the adjacency matrix of $G$ and the matrix of all ones, respectively, it can be shown (see [14, Corollary 2.4]) that

[TABLE]

Therefore, the matrix $\mu(E-A)-E$ will be copositive if $\mu\geq\gamma(G)$ and non-copositive otherwise. Furthermore, in the latter case, the matrix will be closer to the copositive cone as $\mu$ approaches $\gamma(G)$ from the left.

In our tests we considered matrices constructed as follows. Let $G$ be the cycle graph of $n$ nodes whose adjacency matrix, $A_{\text{cycle}}=(a_{ij})\in\mathbb{R}^{n\times n}$, is given component-wise by

[TABLE]

Its clique number is clearly $\gamma(G)=2$. Hence, the matrix

[TABLE]

is copositive for all $\mu\geq 2$ and non-copositive for $\mu<2$. In fact, when $\mu=2$ it coincides with the so-called Horn matrix $H_{n}$ (see, e.g., [13, Section 4]). For instance, the Horn matrix $H_{5}$ takes the form

[TABLE]

Experiments

In our numerical tests we used the parameter setting as

[TABLE]

The trial step-size $\overline{\lambda}_{k}$ in the boosting step of BDCA (Line 8 of Algorithm 1) was chosen to be self-adaptive as in [3]. This technique proceeds as follows. At the first iteration, choose any $\overline{\lambda}_{0}>0$. Then, for $k\geq 1$, if the line search has never been used, we take $\overline{\lambda}_{k}=\overline{\lambda}_{0}$. Otherwise, if the two previous trial step-sizes have been directly accepted (without being reduced by the backtracking step), then the last accepted positive $\lambda$ is scaled by a factor of $\gamma>1$ and used as the current trial step-size. If that is not the case, the trial step-size is set as the last positive value of $\lambda$ accepted in previous iterations. In our tests we used

[TABLE]

In our first numerical experiment, we considered Horn matrices of different sizes, $H_{n}$, for $n\in\{1000,1250,\ldots,5000\}$. For each size, DCA and BDCA were run from the same $100$ starting points randomly generated in the intersection of the non-negative orthant with the unit ball. We stopped the algorithms when $\left\|d_{k}\right\|\leq 10^{-9}$ for the first time. The results are shown in Figure 1, where we can observe that, on average, BDCA was more than 15 times faster than DCA for all sizes. As expected, since Horn matrices are copositive, both algorithms converged to critical points with a positive objective value very close to 0. It is worth to mention that the objective function at the points found by BDCA was usually smaller than at the ones found by DCA. We also show in Figure 2 the percentage of iterations at which the boosting step was activated, that is, when condition in line 7 of Algorithm 1 was fulfilled. We observed that, for all sizes, the linesearch was performed in around the 45% of the iterations. In Figure 3 we show the behavior of both algorithms in a particular instance for testing the copositivity of $H_{1000}$.

In our second experiment we considered matrices of the form $Q_{n}^{\mu}$ as defined in (32). In order to generate hard instances (those which are close to be copositive) we took $\mu:=1.9$. For each size $n\in\{1000,1250,\ldots,5000\}$, DCA and BDCA were run from the same $100$ random starting points generated as in our previous experiment. In this case, we let the algorithms run until they find a negative objective value (which exists because of the non-copositivity of the matrices). We used two stopping criteria, whose results are depicted in Figure 4: on the left, the algorithms were stopped when any negative objective value was found; on the right, the objective value was required to be smaller than $-10^{-4}$. We do not show any results on the second criterion for $n$ greater than $2000$ because DCA becomes extremely slow (for $n=2000$, the instances solved by DCA required more than 5 minutes on average). This time the advantage of BDCA with respect to DCA increased with the size $n$, and was significantly greater than in the previous experiment when the second criterion was used.

Remark 5.1*.*

Clearly, BDCA iterations will be more time consuming than those of DCA due to the need of checking condition in Line 7 and the linesearch procedure. Note that the linesearch requires of function evaluations of $f$, which may be expensive at some applications. For this reason, in our experiments we compare the total CPU running time of the algorithm, as it includes the possibly wasted extra time of the boosting step.

5.2 Solving the $\ell_{1}$ and $\ell_{\infty}$ trust-region subproblem

The trust-region subproblem (TRSP) arises in trust-region methods, which are optimization algorithms that consist in replacing the original objective function by a model which is a good approximation of the original function at the current iterate. The models are usually defined to be a quadratic function, in which case the task consists in solving nonconvex quadratic optimization problems of the type

[TABLE]

where $A$ is an $n\times n$ real symmetric matrix, $b\in\mathbb{R}^{n}$, $r$ is a positive number and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ is any given norm. Of course, the choice of the norm has a serious impact on the difficulty for solving the subproblems.

The application of the DCA for solving the Euclidean trust-region subproblem (the most common choice for the norm) was first proposed in [29] and its convergence was further studied in [16, 31, 33]. Thanks to the structure of the trust-region subproblem, the implementation of the DCA is very simple and efficient, as the iterations are given by (15) and only require matrix-vector products. On the other hand, as observed in [7, Section 7.8] and [12], the choice of the Euclidean norm for the definition of the trust-region has several drawbacks, especially when the problem at hand is box-constrained, in which case the intersection of the Euclidean ball of the trust region with the feasible set has a complicated structure. The choice of the $\ell_{\infty}$ norm $\|\cdot\|_{\infty}$ for the trust-region when the problem involves bounds of the type $l\leq x\leq u$ permits to ensure feasibility by simply requiring similar bounds on the trust-region subproblem. Unfortunately, the $\ell_{1}$ and $\ell_{\infty}$ trust-region subproblem are NP-hard [23], a class of problems for which it is commonly believed that there are no polynomial time algorithms.

In this subsection we compare the performance of DCA and BDCA on the challenging $\ell_{1}$ and $\ell_{\infty}$ trust-region subproblems. To this aim, we consider the DC decomposition (12)-(13) of the (possibly) nonconvex part of the objective function. Thus, given $\sigma>\max\left\{\lambda_{\max}(A),0\right\}$, the subproblems for the $\ell_{1}$ norm take the form

[TABLE]

with

[TABLE]

while for the $\ell_{\infty}$ norm these problems are

[TABLE]

Observe that the feasible set of ($\mathcal{P}_{1}$) can be equivalently written in terms of $2^{n}$ linear constraints of the type $\pm x_{1}\pm x_{2}\pm\cdots\pm x_{n}\leq r$, so the problem is a particular instance of ($\mathcal{P}$). In principle, as mentioned in Remark 3.1, the largest step-size derived in [9] ensuring feasibility could be computed. However, as the feasible set has $2^{n}$ linear constraints, computing (4) would be very time-consuming. Nonetheless, given a point $y_{k}$ and a feasible direction $d_{k}$, an upper bound of the maximum value of the step-size $\overline{\lambda}_{k}$ for ($\mathcal{P}_{1}$) is given by the Euclidean norm, since

[TABLE]

where

[TABLE]

so one must always choose $\overline{\lambda}_{k}\leq\widehat{\lambda}_{k}$ to avoid extra time in checking feasibility. On the other hand, for the $\ell_{\infty}$ norm, the situation is more favorable, as feasibility can be guaranteed whenever

[TABLE]

which coincides with (4), so no time will be wasted in Step 8 of Algorithm 1 if one takes $\overline{\lambda}_{k}$ satisfying the latter inequalities.

Experiments

We used the same parameter setting as in the previous section for Algorithm 1, except for the value of $\gamma$, which was increased to $20$. The matrix $A$ and the vector $b$ were generated with coordinates randomly and uniformly chosen in ${]-1,1[}$, while the value of $r$ was randomly and uniformly chosen in the range ${]0,\sqrt{n}/4[}$ for ($\mathcal{P}_{1}$) and ${]0,1/4[}$ for ($\mathcal{P}_{\infty}$). For each size $n\in\{1000,1500,\ldots,5000\}$, DCA and BDCA were run from the same $100$ starting points randomly picked inside the trust-region. We stopped the algorithms when $\left\|d_{k}\right\|/\|x_{k}\|\leq 10^{-8}$ for the first time. The time comparison for both norms are summarized in Figure 5. BDCA was consistently faster than DCA. On average, it was 3.8 times faster for solving ($\mathcal{P}_{1}$) and 3.65 times faster for ($\mathcal{P}_{\infty}$).

We also show in Figure 6 the percentage of iterations at which the boosting step was activated. We observe that, for all sizes, the linesearch was performed on average in around 80% of the iterations for the case of the $\ell_{1}$ norm, and 40% for the $\ell_{\infty}$ norm.

When applied to ($\mathcal{P}_{1}$), both algorithms obtained the same value in 875 instances, BDCA reached a smaller objective value in 18 instances, while DCA did so in 7. When applied to ($\mathcal{P}_{\infty}$), BDCA obtained a smaller objective value in 461 instances, while DCA did so in 417 (the value was the same in 22 instances). In Figure 7 we show the behavior of both algorithms for a particular instance with $n=1000$, which was chosen so that both algorithms attained the same objective value.

We conclude this section with some comments about the rate of convergence of DCA and BDCA. According to Theorem 4.1, both methods are R-linearly convergent, so by (18),

[TABLE]

for some $\eta\in{]0,1[}$. In Figure 8, for a particular random instance, we have represented the sequence $\{\|x_{k}-\overline{x}\|^{1/k}\}$, which was computed for both algorithms after obtaining the limit point $\overline{x}$ with a higher precision. We observe that the value of $\eta$ is very close to $1$ for DCA, while line-searches clearly helped improving this slow rate. A possible explanation about the bad R-linear convergence rate of DCA can be deduced from the proof of Theorem 4.1. The upper bound obtained there is $\frac{\beta}{\rho+\beta}$, with $\beta>\sigma$. When $\lambda_{\max}(Q)>0$, Assumption 1 holds for $\rho$ as in (14), which would then be equal to $\sigma-\lambda_{\max}(Q)$. Therefore,

[TABLE]

where the last inequality holds since the left term is an increasing function with respect to $\beta$. In our numerical tests we took $\sigma-\lambda_{\max}(Q)=0.01$ (we numerically observed how larger values slowed down both DCA and BDCA). In the random instances shown in Figure 8 $\lambda_{\max}(Q)$ was equal to $25.77$ for the $\ell_{1}$ norm and $25.89$ for the $\ell_{\infty}$ norm, which would give upper bounds on the R-linear convergence rate of $0.9996$ for both problems. According to the figures, these bounds seem to be tight, particularly for the $\ell_{\infty}$ norm.

5.3 Piecewise quadratic problems with box constraints

Finally, we test BDCA on optimization problems with nonsmooth objective function. To this aim, we consider piecewise quadratic problems with linear constraints of the form

[TABLE]

for given $l=(l_{1},\ldots,l_{n}),u=(u_{1},\ldots,u_{n})\in\mathbb{R}^{n}$ and $c^{1},\ldots,c^{m}\in\mathbb{R}^{n}$.

As discussed in [26, § 6.3], the objective function of this problem admits the DC decomposition $\phi(x)=g(x)-h(x)$ with

[TABLE]

Experiments

For our numerical tests, we generated problems of the form (33) as follows. First, the vector $l\in\mathbb{R}^{n}$ was generated with coordinates randomly chosen in ${]-5,5[}$. Then, we randomly generated a vector $e\in\mathbb{R}^{n}$ with coordinates in ${]0,5[}$ and we set $u:=l+e$. Finally, the points $c_{1},\ldots,c_{m}\in\mathbb{R}^{n}$ were generated so that they became unfeasible for all constraints as follows: each component $i\in\{1,\ldots,n\}$ of each of these vectors was randomly picked in one of the intervals ${]l_{i}-10,l_{i}[}$ or ${]u_{i},u_{i}+10[}$. In our experiment, for each $n\in\{100,150,\ldots,1000\}$ and each $m\in\{100,150,\ldots,1000\}$, DCA and BDCA were run from the same $100$ starting points randomly picked inside the feasible set. We used the same parameter setting and stopping criterion for the algorithms as in the previous TRSP experiments. The time comparison is shown in Figure 9, where we represent the median ration between DCA and BDCA among the 100 instances. Detailed results for all 100 runs are shown in Figure 10(a) for fixed $n=500$ and Figure 10(b) for fixed $m=500$. From the results of this comparison we infer that the superiority of BDCA increases as $m$ does, while it stabilizes as $n$ increases. As in previous experiments, we show in Figure 11 the percentage of iterations at which the boosting step was activated. The objective value attained by DCA and BDCA were the same in all instances.

6 Concluding remarks

We have extended the Boosted DC Algorithm for solving linearly constrained DC programming. The algorithm is proved to provide KKT points of the constrained problem. In addition, we have shown why this approach cannot be extended to more general convex constraints. The theoretical results were confirmed by some numerical experiments for testing the copositivity of matrices and for solving $\ell_{1}$ and $\ell_{\infty}$ trust-region subproblems. For copositive matrices, BDCA was on average fifteen times faster than DCA; for non-copositive ones, this advantage was much more superior; and for trust-region subproblems, BDCA was more than three times faster than DCA. We also considered piecewise quadratic problems with box constraints, which thus have nonsmooth objective functions. We observed again that BDCA was faster than DCA and, further, that the advantage was more noticeable as the number of pieces of the objective function increased. Future research includes investigation of alternative approaches to derive a Boosted DCA that permits to address any type of constrained DC programs. It would also be interesting to combine BDCA with the inertial technique employed in [27].

Acknowledgments

FJAA and RC were partially supported by the Ministry of Science, Innovation and Universities of Spain and the European Regional Development Fund (ERDF) of the European Commission (PGC2018-097960-B-C22), and by the Generalitat Valenciana (AICO/2021/165). PTV was supported by Vietnam Ministry of Education and Training Project hosting by the University of Technology and Education, Ho Chi Minh City Vietnam (2023-2024).