Solving a Continuous Multifacility Location Problem by DC Algorithms

TL;DR

This paper introduces a novel DC algorithm with Nesterov's smoothing to efficiently solve complex multifacility location problems reformulated as continuous optimization tasks.

Contribution

It develops a new DC-based algorithm for multifacility location problems, addressing their nonconvex and nondifferentiable challenges with a reformulation and smoothing techniques.

Findings

Algorithm successfully solves artificial and real data instances.

Reformulation improves computational efficiency.

Demonstrates effectiveness of DC algorithms in location problems.

Abstract

The paper presents a new approach to solve multifacility location problems, which is based on mixed integer programming and algorithms for minimizing differences of convex (DC) functions. The main challenges for solving the multifacility location problems under consideration come from their intrinsic discrete, nonconvex, and nondifferentiable nature. We provide a reformulation of these problems as those of continuous optimization and then develop a new DC type algorithm for their solutions involving Nesterov's smoothing. The proposed algorithm is computationally implemented via MATLAB numerical tests on both artificial and real data sets.

Tables

Table 1. Table 1: Comparison between Algorithm 3 and k 𝑘 k -means

Method	Optimal Center ($V^N$)	Cost Function
$k$-means		22.1637
Algorithm 3		22.1352

Table 2. Table 2: Comparison between Algorithm 3 and k 𝑘 k -means

Method	Optimal Center ($V^N$)	Cost Function
$k$-means		403.3966
Algorithm 3		401.7506

Table 3. Table 3: Comparison between Algorithm 3 (combined with k 𝑘 k -means) and standard k 𝑘 k -means

Method	Optimal Center ($V^N$)	Cost Function
$k$-means		288.8348
Algorithm 3 (combined with $k$-means)		286.6523

Table 4. Table 4: Cost comparison between Algorithm 3 (combined with k 𝑘 k -means) and k 𝑘 k -means

Method	Cost Function
$k$-means	16556
Algorithm 3 (combined with $k$-means)	16460

Equations

$$\displaystyle\begin{array}[]{ll}&\mbox{\rm minimize }\mathcal{F}(\mathbf{U},\mathbf{V}):=\sum_{i=1}^{k}\sum_{j=1}^{n}u_{i,j}^{2}\|a_{j}-v_{i}\|\\ &\mbox{\rm subject to }\mathbf{U}\in\mathcal{U}\;\;\mbox{\rm and }\;\mathbf{V}\in\mathbb{R}^{k\times d}.\end{array}$$

$$\text{mbox}\mathrm{minimize}f(x):=g(x)-h(x),x\in {\mathbb{R}}^d.$$

$$P(x;\Omega):=\big{\{}w\in\Omega\;\big{|}\;\|x-w\|=d(x;\Omega)\big{\}},$$

$$d(x;\Omega):=\inf\big{\{}\|x-w\|\;\big{|}\;w\in\Omega\big{\}}.$$

$$\varphi^{*}(v):=\sup\big{\{}\langle v,x\rangle-\varphi(x)\;\big{|}\;x\in\mathbb{R}^{d}\big{\}},\quad v\in\mathbb{R}^{d}.$$

$$\partial\varphi(\bar{x}):=\big{\{}v\in\mathbb{R}^{d}\;\big{|}\;\langle v,x-\bar{x}\rangle\leq\varphi(x)-\varphi(\bar{x})\;\;\mbox{\rm for all }\;x\in\mathbb{R}^{d}\big{\}}.$$

$$v\in\mbox{\rm argmin}\,\big{\{}\varphi(x)-\langle y,x\rangle\;\big{|}\;x\in\mathbb{R}^{d}\big{\}}.$$

$$w\in\mbox{\rm argmin}\,\big{\{}\varphi^{*}(y)-\langle x,y\rangle\;\big{|}\;y\in\mathbb{R}^{d}\big{\}}.$$

$$0\in \partial \phi (v)-y,$$

$$0\in \partial {\phi }^{\ast }(w)-x,$$

$$\Vert \mathbf{U}{\Vert }_F:=\sqrt{\sum _{i=1}^k\sum _{j=1}^n|u_{i,j}{|}^2}.$$

$$\underset{\Vert x\Vert \to \infty }{\lim }\frac{g(x)}{\Vert x\Vert }=\infty .$$

$$c+⟨w,x⟩\le g(x)\text{mbox}\mathrm{forall}x\in {\mathbb{R}}^d.$$

$$\|x\|\big{(}\|w\|+1\big{)}\leq g(x)\;\;\mbox{\rm whenever }\;\|x\|\geq\eta.$$

$$\sup\big{\{}\langle v,x\rangle-g(x)\;\big{|}\;\|x\|\geq\eta\big{\}}\leq-\|x\|\;\;\mbox{\rm for any }\;v\in\mathbb{R}^{d}.$$

$$\sup\big{\{}\langle v,x\rangle-g(x)\;\big{|}\;\|x\|\leq\eta\big{\}}\leq\sup\big{\{}\langle v,x\rangle-\langle w,x\rangle-c\;\big{|}\;\|x\|\leq\eta\big{\}}<\infty,\quad v\in\mathbb{R}^{d}.$$

$$\partial g(\bar{x})\cap \partial h(\bar{x})\ne \varnothing .$$

$$f(x):=x^4-2x^2+2x-3\text{mbox}\mathrm{for}x\in \mathbb{R}.$$

$$x_{l+1}=x_l-t(4x_l^3-4x_l+2)\text{mbox}\mathrm{for}l=0,1,\dots .$$

$$x_{l+1}=\nabla g^{*}(y_{l})=\Big{(}\frac{y_{l}}{4}\Big{)}^{1/3}=\Big{(}\frac{4x_{l}-2}{4}\Big{)}^{1/3}=\Big{(}\frac{2x_{l}-1}{2}\Big{)}^{1/3},\quad l=0,1,\ldots.$$

$$\text{mbox}\mathrm{minimize}f(x_1,x_2):=x_1^4+x_2^2-2x_1^2-|x_2|\text{mbox}\mathrm{over}x=(x_1,x_2)\in {\mathbb{R}}^2.$$

$$\partial h(x)=\big{[}4x_{1},\mbox{\rm sign}(x_{2})\big{]}^{T}\;\;\mbox{\rm for any }\;x=(x_{1},x_{2})\in\mathbb{R}^{2}.$$

$$\text{mbox}\mathrm{minimize}{\phi }_l(x):=g(x)-⟨y_l,x⟩\text{mbox}\mathrm{over}x\in {\mathbb{R}}^2$$

$$f(x):=\Vert x-a\Vert ,x\in {\mathbb{R}}^d,$$

$$f_{\mu}(x):=\dfrac{1}{2\mu}\|x-a\|^{2}-\dfrac{\mu}{2}\bigg{[}d\bigg{(}\dfrac{x-a}{\mu};\mathbb{B}\bigg{)}\bigg{]}^{2}.$$

$$\nabla f_{\mu}(x)=P\bigg{(}\dfrac{x-a}{\mu};\mathbb{B}\bigg{)}\;\;\mbox{\rm and }\;f_{\mu}(x)\leq f(x)\leq f_{\mu}(x)+\dfrac{\mu}{2},$$

$$\displaystyle\begin{array}[]{ll}\mathcal{F}_{\mu}(\mathbf{U},\mathbf{V})&:=\displaystyle\dfrac{1}{2\mu}\sum_{i=1}^{k}\sum_{j=1}^{n}u_{i,j}^{2}\|a_{j}-v_{i}\|^{2}-\displaystyle\dfrac{\mu}{2}\sum_{i=1}^{k}\sum_{j=1}^{n}u_{i,j}^{2}\bigg{[}d\bigg{(}\dfrac{a_{j}-v_{i}}{\mu};\mathbb{B}\bigg{)}\bigg{]}^{2}\\ &=\mathcal{G}_{\mu}(\mathbf{U},\mathbf{V})-\mathcal{H}_{\mu}(\mathbf{U},\mathbf{V}),\end{array}$$

$$\displaystyle\begin{array}[]{ll}&\mathcal{G}_{\mu}(\mathbf{U},\mathbf{V}):=\displaystyle\dfrac{1}{2\mu}\sum_{i=1}^{k}\sum_{j=1}^{n}u_{i,j}^{2}\|a_{j}-v_{i}\|^{2},\\ &\mathcal{H}_{\mu}(\mathbf{U},\mathbf{V}):=\displaystyle\dfrac{\mu}{2}\sum_{i=1}^{k}\sum_{j=1}^{n}u_{i,j}^{2}\bigg{[}d\bigg{(}\dfrac{a_{j}-v_{i}}{\mu};\mathbb{B}\bigg{)}\bigg{]}^{2}.\end{array}$$

$$\displaystyle\begin{array}[]{ll}&\mbox{\rm minimize }\;\mathcal{F}_{\mu}(\mathbf{U},\mathbf{V}):=\mathcal{G}_{\mu}(\mathbf{U},\mathbf{V})-\mathcal{H}_{\mu}(\mathbf{U},\mathbf{V})\;\;\mbox{\rm as }\;\mu>0\\ &\mbox{\rm subject to }\;\mathbf{U}\in\mathcal{U}=\Delta^{n}\cap\{0,1\}^{k\times n}\;\;\mbox{\rm and }\;\mathbf{V}\in\mathbb{R}^{k\times d},\end{array}$$

$$\mathcal{F}_{\mu}(\mathbf{U},\mathbf{V})=\dfrac{1}{2\mu}\sum_{i=1}^{k}\sum_{j=1}^{n}u_{i,j}\|a_{j}-v_{i}\|^{2}-\dfrac{\mu}{2}\sum_{i=1}^{k}\sum_{j=1}^{n}u_{i,j}\bigg{[}d\bigg{(}\dfrac{a_{j}-v_{i}}{\mu};\mathbb{B}\bigg{)}\bigg{]}^{2}$$

Full text

**Solving a Continuous Multifacility Location Problem by DC Algorithms

** Anuj Bajaj111 Department of Mathematics, Wayne State University, Detroit, Michigan 48202, USA (anuj.bajaj@wayne.edu ). Research of this author was partly supported by the USA National Science Foundation under grant DMS-1808978 and by the USA Air Force Office of Scientific Research grant #15RT04., Boris S. Mordukhovich222Department of Mathematics, Wayne State University, Detroit, Michigan 48202, USA (boris@math.wayne.edu). Research of this author was partly supported by the USA National Science Foundation under grants DMS-1512846 and DMS-1808978, by the USA Air Force Office of Scientific Research grant #15RT04, and by Australian Research Council under grant DP-190100555., Nguyen Mau Nam333Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR 97207, USA (mnn3@pdx.edu). Research of this author was partly supported by the USA National Science Foundation under grant DMS-1716057., and Tuyen Tran444Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR 97207, USA (tuyen2@pdx.edu)..

**Abstract. The paper presents a new approach to solve multifacility location problems, which is based on mixed integer programming and algorithms for minimizing differences of convex (DC) functions. The main challenges for solving the multifacility location problems under consideration come from their intrinsic discrete, nonconvex, and nondifferentiable nature. We provide a reformulation of these problems as those of continuous optimization and then develop a new DC type algorithm for their solutions involving Nesterov’s smoothing. The proposed algorithm is computationally implemented via MATLAB numerical tests on both artificial and real data sets.

Key words. Mixed integer programming, multifacility location, difference of convex functions, Nesterov’s smoothing, the DCA

AMS subject classifications. 49J52, 49J53, 90C31**

1 Introduction and Problem Formulation

In the 17th century, Pierre de Fermat posed a problem of finding for a point that minimizes the sum of its Euclidean distances to three given points in the plane. The problem was soon solved by Evangelista Torricelli, and it is now known as the Fermat-Torricelli problem. This problem and its extended version that involves a finite number of points in higher dimensions are examples of continuous single facility location problems. Over the years several generalized models of the Fermat-Torricelli type have been introduced and studied in the literature with practical applications to facility location decisions; see [8, 12, 14, 16, 18, 19] and the references therein. An important feature of single facility location problems and the problems studied in the aforementioned references is that only one center/server has to be found to serve a finitely many demand points/customers.

However, numerous practical applications lead to formulations of facility location problems in which more than one center must be found to serve a finite number of demand points. Such problems are referred to as multifacility location problems (MFLPs). Given a finite number of demand points $a_{1},\ldots,a_{n}$ in $\mathbb{R}^{d}$, we consider here the facility location in which $k$ centers $v_{1},\ldots,v_{k}$ ($1\leq k\leq n$) in $\mathbb{R}^{d}$ need to be found to serve these demand points by assigning each of them to its nearest center and minimizing the total distances from the centers to the assigned demand points. In the case where $k=1$, this problem reduces to the generalized Fermat-Torricelli problem of finding a point that minimizes the sum of the distances to a finite number of given points in $\mathbb{R}^{d}$.

Let us formulate the problem under consideration in this paper as the following problem of mixed integer programming with nonsmooth objective functions. It is convenient to use a variable $k\times d$-matrix $\mathbf{V}$ with $v_{i}$ as its $i$th row to store the centers to be found. We also use another variable $k\times n$-matrix $\mathbf{U}=[u_{i,j}]$ with $u_{i,j}\in\{0,1\}$ and $\sum_{i=1}^{k}u_{i,j}=1$ for $j=1,\ldots,n$ to assign demand points to the centers. The set of all such matrices is denoted by $\mathcal{U}$. Note that $u_{i,j}=1$ if the center $v_{i}$ is assigned to the demand point $a_{j}$ while $\sum_{i=1}^{k}u_{i,j}=1$ means that the demand point $a_{j}$ is assigned to only one center. Our goal is to solve the constrained optimization problem formulated as follows:

[TABLE]

Taking into account that $u_{i,j}\in\{0,1\}$, it is convenient to use $u^{2}_{i,j}$ instead of $u_{i,j}$ in the definition of the objective function $\mathcal{F}$; See Section 4.

Note that a similarly looking problem was considered by An, Minh and Tao [4] for different purposes. The main difference between our problem (1.3) and the one from [4] is that in [4] the squared Euclidean norm is used instead of the Euclidean norm in our formulation. From the point of applications this difference is significant; namely, using the Euclidean norm allows us to model the total distance in supply delivery, while using the squared Euclidean norm is meaningful in clustering. Mathematically these two problems are essentially different as well. In addition to the challenging discrete nature and nonconvexity that both problems share, the objective function of our multifacility location problem (1.3) is nondifferentiable in contrast to [4]. This is yet another serious challenge from both theoretical and algorithmic viewpoints. Observe also that for $k=1$ our problem becomes the aforementioned generalized Fermat-Torricelli problem that does not have a closed-form solution, while the problem considered in [4] reduces to the standard problem of minimizing the sum of squares of the Euclidean distances to the demand points. The latter has a simple closed form solution given by the mean of the data points.

In this paper we develop the following algorithmic procedure to solve the formulated nonsmooth problem (1.3) of mixed integer programming:

(i) Employ Nesterov’s smoothing to approximate the nonsmooth objective function in (1.3) by a family of smooth functions, which are represented as differences of convex ones.

(ii) Enclose the obtained smooth discrete problems into constrained problems of continuous DC optimization and then approximate them by unconstrained ones while using penalties.

(iii) Solve the latter class of problems by developing an appropriate modification of the algorithm for minimizing differences of convex functions known as the DCA.

As a result of all the three steps above, we propose a new algorithm for solving the class of multifacility local problems of type (1.3), verify its efficiency and implementation with MATLAB numerical tests on both artificial and real data sets.

Recall that the early developments on the DCA trace back to the work by Tao in 1986 with more recent results presented in [2]–[5], [24], [25], and the bibliographies therein. Nesterov’s smoothing technique was introduced in his seminal paper [20] and was further developed and applied in many great publications; see, e.g., [21, 22] for more details and references. The combination of these two important tools provides an effective way to deal with nonconvexity and nondifferentiability in many optimization problems encountered in facility location, machine learning, compressed sensing, and imaging. It is demonstrated in this paper in solving multifacility location problems of type (1.3).

The rest of this paper is organized as follows. Section 2 contains the basic definitions and some preliminaries, which are systematically employed in the text. In Section 3 we briefly overview two versions of the DCA, discuss their convergence, and present two examples that illustrate their performances.

Section 4 is devoted to applying Nesterov’s smoothing technique to the objective function of the multifacility location problem (1.3) and constructing in this way a smooth approximation of the original problem by a family of DC ones. Further, we reduce the latter smooth DC problems of discrete constrained optimization to unconstrained problems by using an appropriate penalty function method. Finally, the obtained discrete optimization problems are enclosed here into the DC framework of unconstrained continuous optimization.

In Section 5 we proposed, based on the above developments, a new algorithm to solve the multifacility location problem (1.3) by applying the updated version of the DCA taken from Section 3 to the smooth DC problems of continuous optimization constructed in Section 4. The proposed algorithm is implemented in this section to solving several multifacility problems arising in practical modeling. Section 6 summarizes the obtained results and discusses some directions of future research.

2 Basic Definitions and Preliminaries

For the reader’s convenience, in this section we collect those basic definitions and preliminaries, which are largely used throughout the paper; see the books [9, 15, 23] for more details and proofs of the presented results.

Consider the difference of two convex functions $g-h$ on a finite-dimensional space and assume that $g\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}:=(-\infty,\infty]$ is extended-real-valued while $h\colon\mathbb{R}^{d}\to\mathbb{R}$ is real-valued on $\mathbb{R}^{d}$. Then a general problem of DC optimization is defined by:

[TABLE]

Note that problem (2.1) is written in the unconstrained format, but—due to the allowed infinite value for $g$—it actually contains the domain constraint $x\in\mbox{\rm dom}\,(g):=\{u\in\mathbb{R}^{d}\;|\;g(u)<\infty\}$. Furthermore, the explicit constraints of the type $x\in\Omega$ given by a nonempty convex set $\Omega\subset\mathbb{R}^{d}$ can be incorporated into the format of (2.1) via the indicator function $\delta_{\Omega}(x)$ of $\Omega$, which equals [math] for $x\in\Omega$ and $\infty$ otherwise. The representation $f=g-h$ is called a DC decomposition of $f$. Note that the class of DC functions is fairly large and include many nonconvex functions important in optimization. We refer the reader to the recent book [13] with the commentaries and bibliographies therein for various classes of nonconvex optimization problems that can be represented in the DC framework (2.1).

Considering a nonempty (may not be convex) set $\Omega\subset\mathbb{R}^{d}$ and a point $x\in\mathbb{R}^{d}$, define the Euclidean projection of $x$ to $\Omega$ by

[TABLE]

where $d(x;\Omega)$ stands for the Euclidean distance from $x$ to $\Omega$, i.e.,

[TABLE]

Observe that $P(x;\Omega)\neq\emptyset$ for closed sets $\Omega$ while being always a singleton if the set $\Omega$ is convex.

Given further an extended-real-valued and generally nonconvex function $\varphi\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}$, the Fenchel conjugate of $\varphi$ is defined by

[TABLE]

If $\varphi$ is proper, i.e., $\mbox{\rm dom}\,(\varphi)\neq\emptyset$, its Fenchel conjugate $\varphi^{*}\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}$ is automatically convex.

The subdifferential of $\varphi\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}$ at $\bar{x}\in\mbox{\rm dom}\,(\varphi)$ is the set of subgradients given by

[TABLE]

If $\bar{x}\notin\mbox{\rm dom}\,(\varphi)$, we let $\partial\varphi(\bar{x}):=\emptyset$. Recall that for functions $\varphi$ differentiable at $\bar{x}$ with the gradient $\nabla\varphi(\bar{x})$ we have $\partial\varphi(\bar{x})=\{\nabla\varphi(\bar{x})\}$.

The following proposition gives us a two-sided relationship between the Fenchel conjugates and subgradients of convex functions.

Proposition 2.1

Let $\varphi\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}$ be a proper, lower semicontinuous, and convex function. Then $v\in\partial\varphi^{*}(y)$ if and only if

[TABLE]

We have furthermore that $w\in\partial\varphi(x)$ if and only if

[TABLE]

Proof. To verify the first assertion, suppose that (2.5) is satisfied and then get $0\in\partial\psi(v)$, where $\psi(x):=\varphi(x)-\langle y,x\rangle$ as $x\in\mathbb{R}^{d}$. It tells us that

[TABLE]

and hence $y\in\partial\varphi(v)$, which is equivalent to $v\in\partial\varphi^{*}(y)$ due to the biconjugate relationship $\varphi^{**}=\varphi$ valued under the assumptions made.

In the opposite way, assuming $v\in\partial\varphi^{*}(y)$ gives us by the proof above that $0\in\partial\psi(v)$, which clearly yields (2.5) and thus justifies the first assertion.

To verify the second assertion, suppose that (2.6) holds and then get $0\in\partial\psi(w)$, where $\psi(y):=\varphi^{*}(y)-\langle x,y\rangle$ as $y\in\mathbb{R}^{d}$. This clearly implies that

[TABLE]

and hence $x\in\partial\varphi^{*}(w)$, which is equivalent to $w\in\partial\varphi(x)$ due to the biconjugate relationship. The proof of the opposite implication in (2.6) is similar to the one given above. $\hfill\square$

Finally in this section, recall that for a given variable matrix $\mathbf{U}\in\mathbb{R}^{k\times n}$ as in the optimization problem (1.3), the Frobenius norm on $\mathbb{R}^{k\times n}$ is defined by

[TABLE]

3 Overview of the DCA and Some Examples

In this section we first briefly overview two algorithms of the DCA type to solve DC problems (2.1) while referring the reader to [24, 25] for more details and further developments. Then we present numerical examples illustrating both algorithms.

Algorithm 1: DCA-1.

[TABLE]

Since the convex function $h\colon\mathbb{R}^{d}\to\mathbb{R}$ in (2.1) is real-valued on the whole space $\mathbb{R}^{d}$, we always have $\partial h(x)\neq\emptyset$ for all $x\in\mathbb{R}^{d}$. At the same time, the other convex function $g\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}$ in (2.1) is generally extended-real-valued, and so the subdifferential of its conjugate $g^{*}$ may be empty. Let us present an efficient condition that excludes this possibility. Recall that a function $g\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}$ is coercive if

[TABLE]

Proposition 3.1

Let $g\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}$ be a proper, lower semicontinuous, and convex function. If in addition $g$ is coercive, then $\partial g^{*}(v)\neq\emptyset$ for all $v\in\mathbb{R}^{d}$.

Proof. Since $g$ is proper, the conjugate function $g^{*}$ takes values in $(-\infty,\infty]$ being convex on $\mathbb{R}^{d}$. Taking into account that $g$ is also and lower semicontinuous and invoking the aforementioned biconjugate relationship, we find $w\in\mathbb{R}^{d}$ and $c\in\mathbb{R}$ such that

[TABLE]

The coercivity property of $g$ ensures the existence of $\eta>0$ for which

[TABLE]

It follows furthermore that

[TABLE]

By using (3.1), we arrive at the estimates

[TABLE]

This tells us that $g^{*}(v)<\infty$, and therefore $\mbox{\rm dom}\,(g^{*})=\mathbb{R}^{d}$. Since $g^{*}$ is a convex function with finite values, it is continuous on $\mathbb{R}^{d}$ and hence $\partial g^{*}(v)\neq\emptyset$ for all $v\in\mathbb{R}^{d}$. $\hfill\square$

To proceed further, recall that a function $\varphi\colon\mathbb{R}^{d}\to\overline{\mathbb{R}}$ is $\gamma$-convex with a given modulus $\gamma\geq 0$ if the function $\psi(x):=\varphi(x)-\frac{\gamma}{2}\|x\|^{2}$ as $x\in\mathbb{R}^{d}$ is convex on $\mathbb{R}^{d}$. If there exists $\gamma>0$ such that $\varphi$ is $\gamma-$convex, then $\varphi$ is called strongly convex on $\mathbb{R}^{d}$.

We also recall that a vector $\bar{x}\in\mathbb{R}^{d}$ is a stationary point of the DC function $f$ from (2.1) if

[TABLE]

The next result, which can be derived from [24, 25], summarizes some convergence results of the DCA. Deeper studies of the convergence of this algorithm and its generalizations involving the Kurdyka-Lojasiewicz (KL) inequality are given in [2, 3].

Theorem 3.2

Let $f$ be a DC function taken from (2.1), and let $\{x_{l}\}$ be an iterative sequence generated by Algorithm 1. The following assertions hold:

1. (a)

The sequence $\{f(x_{l})\}$ is always monotone decreasing. 2. (b)

Suppose that $f$ is bounded from below, that $g$ is lower semicontinuous and $\gamma_{1}$-convex, and that $h$ is $\gamma_{2}$-convex with $\gamma_{1}+\gamma_{2}>0$. If $\{x_{l}\}$ is bounded, then the limit of any convergent subsequence of $\{x_{l}\}$ is a stationary point of $f$.

In many practical applications of Algorithm 1, for a given DC decomposition of $f$ it is possible to find subgradient vectors from $\partial h(x_{l})$ based on available formulas and calculus rules of convex analysis. However, it may not be possible to explicitly calculate an element of $\partial g^{*}(y_{l})$. Such a situation requires either constructing a more suitable DC decomposition of $f$, or finding $x_{l+1}\in\partial g^{*}(y_{l})$ by using the description of Proposition 2.1. This leads us to the following modified version of the DCA.

Algorithm 2: DCA-2.

[TABLE]

Let us now present two examples illustrating the performances of Algorithms 1 and 2. The first example concerns a polynomial function of one variable.

Example 3.3

Consider the function $f\colon\mathbb{R}\to\mathbb{R}$ given

[TABLE]

This function admits the DC representation $f=g-h$ with $g(x):=x^{4}$ and $h(x):=2x^{2}-2x+3$. To minimize $f$, apply first the gradient method with constant stepsize. Calculating the derivative of $f$ is $f^{\prime}(x)=4x^{3}-4x+2$ and picking any starting point $x_{0}\in\mathbb{R}$, we get the sequence of iterates

[TABLE]

constructed by the gradient method with stepsize $t>0$. The usage of the DC Algorithm 1 (DCA-1) gives us $y_{l}=\nabla h(x_{l})=4x_{l}-2$ and then $g^{*}(x)=3(x/4)^{4/3}$ with $\nabla g^{*}(x)=(x/4)^{1/3}$. Thus the iterates of DCA-1 are as follows:

[TABLE]

Figure 3.1 provides the visualization and comparison between the DCA-1 and the gradient method. It shows that for $x_{0}=0$ and $t=0.01$ the DCA-1 exhibits much faster convergence.**

The next two-dimensional example illustrates the performance of the DCA-2.

Example 3.4

Consider the nonsmooth optimization problem defined by

[TABLE]

The graph of the function $f$ is depicted in Figure 3.3. Observe that this function has four global minimizers, which are $(1,0.5),(1,-0.5),(-1,0.5)$, and $(-1,-0.5)$. It is easy to see that $f$ admits a DC representation $f=g-h$ with $g(x_{1},x_{2}):=x_{1}^{4}+x_{2}^{2}$ and $h(x_{1},x_{2}):=2x_{1}^{2}+|x_{2}|$. We get the gradient $\nabla g(x)=[4x_{1}^{3},2x_{2}]^{T}$ and the Hessian

$\nabla^{2}g(x)=\begin{bmatrix}12{x_{1}}^{2}&0\\ 0&2\end{bmatrix}$,

while an explicit formula to calculate $\partial g^{*}(y_{l})$ is not available. Let us apply the DCA-2 to solve this problem. The subdifferential of $h$ is calculated by

[TABLE]

Having $y_{l}$, we proceed with solving the subproblem

[TABLE]

by the classical Newton method with $\nabla^{2}\varphi_{l}(x)=\nabla^{2}g(x)$ and observe that the DCA-2 shows its superiority in convergence with different choices of initial points. Figure 3.3 presents the results of computation by using the DCA-2 with the starting point $x_{0}=(-2,2)$ and employing the Newton method with $\varepsilon=10^{-8}$ to solve subproblem (3.2).**

4 Smooth Approximation by Continuous DC Problems

In this section we first employ and further develop Nesterov’s smoothing technique for the case of multifacility location problem (1.3). Then we enclose the family of DC mixed integer programs obtained in this way into a class of smooth DC problems of continuous optimization. The suggested procedures are efficiently justified by deriving numerical estimates expressed entirely via the given data of the original problem (1.3).

We begin with recalling the following useful result on Nesterov’s smoothing related to the problem under consideration, which is taken from [17, Proposition 2.2].

Proposition 4.1

Given any $a\in\mathbb{R}^{d}$ and $\mu>0$, a Nesterov smoothing approximation of the function $f\colon\mathbb{R}^{d}\rightarrow\mathbb{R}$ defined by

[TABLE]

admits the smooth DC representation

[TABLE]

Furthermore, we have the relationships

[TABLE]

where $\mathbb{B}\subset\mathbb{R}^{d}$ is the closed unit ball, and where $P$ stands for the Euclidean projection (2.2).

Using Proposition 4.1 allows us to approximate the objective function $\mathcal{F}$ in (1.3) by a smooth DC function $\mathcal{F}_{\mu}$ as $\mu>0$ defined as follows:

[TABLE]

where $\mathcal{G}_{\mu},\mathcal{H}_{\mu}\colon\mathbb{R}^{k\times n}\times\mathbb{R}^{k\times d}\to\mathbb{R}$ are given by

[TABLE]

This leads us to the construction of the following family of smooth approximations of the main problem (1.3) defined by

[TABLE]

where $\Delta^{n}$ is the the $n$th Cartesian degree of the $(k-1)$-simplex $\Delta:=\{y\in[0,1]^{k}\;|\;\sum_{i=1}^{k}y_{i}=1\}$, which is a subset of $\mathbb{R}^{k}$.

Observe that for each $\mu>0$ problem (4.5) is of discrete optimization, while our intention is to convert it to a family of problems of continuous optimization for which we are going to develop and implement a DCA-based algorithm in Section 5.

The rest of this section is devoted to deriving two results, which justify such a reduction. The first theorem allows us to verify the existence of optimal solutions to the constrained optimization problems that appear in this procedure. It is required for having well-posedness of the algorithm construction.

Theorem 4.2

Let $(\overline{\mathbf{U}},\overline{\mathbf{V}})$ be an optimal solution to problem (4.5). Then for any $\mu>0$ we have $\overline{\mathbf{V}}\in\mathcal{B}$, where $\mathcal{B}:=\prod_{i=1}^{k}B_{i}$ is the Cartesian product of the $k$ Euclidean balls $B_{i}$ centered at $0\in\mathbb{R}^{d}$ with radius $r:=\sqrt{\sum_{j=1}^{n}\|a_{j}\|^{2}}$ that contain the optimal centers $\bar{v}_{i}$ for each index $i=1,\ldots,k$.

Proof. We can clearly rewrite the objective function in (4.5) in the form

[TABLE]

due to interchangeability between $u_{i,j}^{2}$ and $u_{i,j}$. Observe that $\mathcal{F}_{\mu}(\mathbf{U},\mathbf{V})$ is differentiable on $R^{k\times n}\times\mathbb{R}^{k\times d}$. Employing the classical Fermat rule in (4.5) with respect to $V$ gives us $\nabla_{\mathbf{V}}\mathcal{F}_{\mu}(\overline{\mathbf{U}},\overline{\mathbf{V}})=0$. To calculate this partial gradient, we need some clarification for the second term in (4.6), which is differentiable as a whole while containing the nonsmooth distance function (2.3). The convexity of the distance function in the setting of (4.6) allows us to apply the subdifferential calculation of convex analysis (see, e.g., [15, Theorem 2.39]) and to combine it with an appropriate chain rule to handle the composition in (4.6). Observe that the distance function square in (4.6) is the composition of the nondecreasing convex function $\varphi(t):=t^{2}$ on $[0,\infty)$ and the distance function to the ball $\mathbb{B}$. Thus the chain rule from [15, Corollary 2.62] is applicable. Thus, we can show that $d^{2}(\cdot;\mathbb{B})$ is differentiable with

[TABLE]

Using (4.7), we consider the following two cases:

Case 1: $(a_{j}-\bar{v}_{i})/\mu\in\mathbb{B}$ for the fixed indices $i\in\{1,\ldots,k\}$ and $j\in\{1,\ldots,n\}$. Then

[TABLE]

which gives us

[TABLE]

for the corresponding partial derivatives of $\mathcal{F}_{\mu}$.

Case 2: $(a_{j}-\bar{v}_{i})/\mu\notin\mathbb{B}$ for the fixed indices $i\in\{1,\ldots,k\}$ and $j\in\{1,\ldots,n\}$. In this case we have

[TABLE]

Thus in both cases above it follows from the stationary condition $\nabla_{\mathbf{V}}\mathcal{F}_{\mu}(\overline{\mathbf{U}},\overline{\mathbf{V}})=0$ that

[TABLE]

since we have $\sum\limits_{j=1}^{n}\bar{u}_{i,j}>0$ due to the nonemptiness of the clusters. Then the classical Cauchy-Schwarz inequality leads us to the estimates

[TABLE]

which therefore verify all the conclusions of this theorem. $\hfill\square$

Our next step is to enclose each discrete optimization problem (4.5) into the corresponding one of continuous optimization. For the reader’s convenience if no confusion arises, we keep the same notation ${\mathbf{U}}$ for all the $k\times n$-matrices without the discrete restrictions on their entries. Define now the function $\mathcal{P}\colon\mathbb{R}^{k\times n}\to\mathbb{R}$ by

[TABLE]

and observe that this function is concave on $\mathbb{R}^{k\times n}$ with $\mathcal{P}(\mathbf{U})\geq 0$ whenever $\mathbf{U}\in\Delta^{n}$. Furthermore, we have the representations

[TABLE]

for the set of feasible $k\times n$-matrices $\mathcal{U}$ in the original problem (1.3). Employing further the standard penalty function method allows us to eliminate the most involved constraint on ${\mathbf{U}}$ in (4.8) given by the function ${\mathcal{P}}$. Taking the penalty parameter $\alpha>0$ sufficiently large and using the smoothing parameter $\mu>0$ sufficiently small, consider the following family of continuous optimization problems:

[TABLE]

Observe that Theorem 4.2 ensures the existence of feasible solutions to problem (4.11) and hence optimal solutions to this problem by the Weierstrass theorem due to the continuity of the objective functions therein and the compactness of the constraints sets $\Delta^{n}$ and ${\mathcal{B}}$.

Let us introduce yet another parameter $\rho>0$ ensuring a DC representation of the objective function in (4.11) as follows:

[TABLE]

where the function $\mathcal{G}(\mathbf{U},\mathbf{V}):=\dfrac{\rho}{2}\|(\mathbf{U},\mathbf{V})\|^{2}$ is obviously convex, and

[TABLE]

Since $\mathcal{H}_{\mu}(\mathbf{U},\mathbf{V})-\alpha\mathcal{P}(\mathbf{U})$ is also convex as $\alpha>0$, we are going to show that for any given number $\mu>0$ it is possible to determine the values of the parameter $\rho>0$ such that the function $\dfrac{\rho}{2}\|(\mathbf{U},\mathbf{V})\|^{2}-\mathcal{G}_{\mu}(\mathbf{U},\mathbf{V})$ is convex under an appropriate choice of $\rho$. This would yield the convexity of $\mathcal{H}(\mathbf{U},\mathbf{V})$ and therefore would justify a desired representation of the objective function in (4.11). The following result gives us a precise meaning of this statement, which therefore verify the required reduction of (4.11) to DC continuous optimization.

Theorem 4.3

The function

[TABLE]

is convex on $\Delta^{n}\times\mathcal{B}$ provided that

[TABLE]

where $\xi:=r+\underset{1\leq j\leq n}{\mbox{max}}\|a_{j}\|$ and $r:=\sqrt{\sum_{j=1}^{n}\|a_{j}\|^{2}}$.

Proof. Consider the function $\mathcal{G}_{1}(\mathbf{U},\mathbf{V})$ defined in (4.13) for all $(\mathbf{U},\mathbf{V})\in\Delta^{n}\times\mathcal{B}$ and deduce by elementary transformations directly from its construction that

[TABLE]

Next we define the functions $\gamma_{i,j}\colon\mathbb{R}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ for all $i=1,\ldots,k$ and $j=1,\ldots,n$ by

[TABLE]

and show that each of these functions is convex on the set $\{u_{i,j}\in[0,1],\,\|v_{i}\|\leq r\}$, where $r>0$ is taken from Theorem 4.2.

To proceed, consider the Hessian matrix of each function in (4.15) given by

[TABLE]

and calculate its determinant $\det(J_{\gamma_{i,j}}(u_{i,j},v_{i}))$ by

[TABLE]

It follows from the well-known second-order characterization of the convexity that the function $\gamma_{i,j}(u_{i,j},v_{i})$ is convex on $\{u_{i,j}\in[0,1],\,\|v_{i}\|\leq r\}$ if $\det(J_{\gamma_{i,j}}(u_{i,j},v_{i}))\geq 0$. Using [5, Theorem 1] gives us the estimate

[TABLE]

Then we get from the construction of ${\mathcal{B}}$ in Theorem 4.2 that $0<\|v_{i}-a_{j}\|\leq\|v_{i}\|+\|a_{j}\|\leq r+\underset{1\leq j\leq n}{\mbox{max}}\|a_{j}\|=:\xi$, and therefore

[TABLE]

It allows us to deduce from the aforementioned condition for the convexity of $\gamma_{i,j}(u_{i,j},v_{i})$ that we do have this convexity if $\rho$ satisfies the estimate (4.14). $\hfill\square$

5 Design and Implementation of the Solution Algorithm

Based on the developments presented in the previous sections and using the established smooth DC structure of problem (4.11) with the subsequent $\rho-$parameterization of the objective function therein as ${\mathcal{G}}({\mathbf{U}},{\mathbf{V}})-{\mathcal{H}}({\mathbf{U}},{\mathbf{V}})$, we are now ready to propose and implement a new algorithm for solving this problem involving both DCA-2 and Nesterov’s smoothing.

To proceed, let us present the problem under consideration in the equivalent unconstrained format by using the infinite penalty via the indicator function:

[TABLE]

where $\mathcal{B}$, $\Delta$, and $\rho$ are taken from Section 4.

We first explicitly compute the gradient of the convex function $\mathcal{H}(\mathbf{U},\mathbf{V})$ in (5.19). Denoting

[TABLE]

we have $\mathcal{Y}=\nabla\mathcal{H}_{\mathbf{U}}(\mathbf{U},\mathbf{V})$ and $\mathcal{Z}=\nabla\mathcal{H}_{\mathbf{V}}(\mathbf{U},\mathbf{V})$. Thus for each $i=1,\ldots,k$ and $j=1,\ldots,n$ the $(j,i)-$entry of the matrix $\mathcal{Y}$ and the $i$th row of the matrix $\mathcal{Z}$ are

[TABLE]

respectively. Let us now describe the proposed algorithm for solving the DC program (5.19) and hence the original problem (1.3) of multifacility location. The symbols $\mathcal{Y}^{l-1}_{[j,:]}$ and $\mathcal{Z}^{l-1}_{i}$ in this description represents the $j$th row of the matrix $\mathcal{Y}$ and the $i$th row of the matrix $\mathcal{Z}$ at the $l$th iteration, respectively. Accordingly we use the symbols $\mathbf{U}^{l}_{[:,j]}$ and $\mathbf{V}^{l}_{i}$. Recall also that the Frobenius norm of the matrices in this algorithm is defined in (2.7).

Algorithm 3: Solving Multifacility Location Problems.

[TABLE]

Next we employ Algorithm 3 to solve several multifacility location problems of some practical meaning. By trial and error we verify that the values chosen for $\mu$ determine the performance of the algorithm for each data set. It can be seen that very small values of the smoothing parameter $\mu$ may prevent the algorithm from clustering, and thus we gradually decrease these values. This is done via multiplying $\mu$ by some number $0<\beta<1$ and stopping when $\mu<\mu_{f}$. Note also that in the implementation of our algorithm we use the standard approach of choosing $\mathbf{U}^{0}$ by computing the distance between the point in question and each group center $\mathbf{V}^{0}$ and then by classifying this point to be in the group whose center is the closest to it by assigning the value of $1$, while otherwise we assign the value of [math].

Let us now present several numerical examples, where we compute the optimal centers by using Algorithm 3 via MATLAB calculations. Fix in what follows the values of $\mu=0.5$, $\beta=0.85$, $\varepsilon=10^{-6}$, $\mu_{f}=10^{-6}$, $\alpha=30$, and $\rho=30$ unless otherwise stated. The objective function is the total distance from the centers to the assigned data point. Note that this choice of the objective function seems to be natural from practical aspects in, e.g., airline and other transportation industries, where the goal is to reach the destination via the best possible route available. This reflects minimizing the transportation cost.

In the following examples we implement the standard $k$-means algorithm in MATLAB using the in-built function kmeans().

Example 5.1

Let us consider a data set with $14$ entries in $\mathbb{R}^{2}$ given by

[TABLE]

with the initial data defined by

[TABLE]

Employing Algorithm 3, we obtain the optimal centers as depicted in Table 1 and Figure 5.1.

Table 1 shows that the proposed Algorithm 3 is marginally better for the given data in comparison to the classical $k$-means approach in terms of the objective function.**

Example 5.2

In this example we test our algorithm on a dataset $\mathbf{X}$ containing 10 distinct points on each boundary of 4 balls of radius $r=0.3$ centered at (2,2), (4,2), (4,4), and (2,4). We generate the points as follows:

[TABLE]

where $C_{i}$ and $r_{i}$ are the center and radius of each ball respectively; see [19]. Typically the centroids are the centers of the balls. Choosing a random point from the boundary of each ball for the initial centers, this algorithm converges to the optimal solution

[TABLE]

Its visualization is shown in Figure 5.2.

Note that a drawback in employing the random approach to choose the initial cluster $\mathbf{V}^{0}$ in Example 5.2 is the need of having prior knowledge about the data. Typically it may not be plausible to extract such an information from large unpredictable real life datasets.

In the next example we choose the initial cluster by the process of random selection and see its effect on the optimal centers. Then the results obtained in this way by Algorithm 3 are compared with those computed by the $k$-means approach.

Example 5.3

Let $\mathbf{X}$ be $200$ standard normally distributed random datapoints in $\mathbb{R}^{2}$, and let the initial data be given by

[TABLE]

We obtain the optimal centers as outlined in Table 2.

Observe from Table 2 that the proposed Algorithm 3 is better for the given data in comparison to the standard $k$-means approach. In addition, our approach gives a better approximation for the optimal center as shown in Figure 5.3.

Note that a real-life data may not be as efficiently clustered as in Example 5.3. Thus a suitable selection of the initial cluster $\mathbf{V}^{0}$ is vital for the convergence of the DCA based algorithms. In the next Example 5.4 we select $\mathbf{V}^{0}$ in Algorithm 3 by using the standard $k$-means method. The results achieved by our Algorithm 3 are again compared with those obtained by using the $k$-means approach.

Example 5.4

Consider the dataset $\mathbf{X}$ consisting of the latitudes and longitudes of $50$ most populous cities in the USA555Available at https://en.wikipedia.org/wiki/List of United States cities by population with

[TABLE]

By using Algorithm 3 we obtain the following optimal centers as given in Table 3.

We see that Algorithm 3 (combined with $k$-means) in which the initial cluster $\mathbf{V}^{0}$ is selected by using $k$-means method performs better in comparison to the standard $k$-means approach (Table 3). Moreover, it gives us optimal centers as depicted in Figure 5.4.

In the next example we efficiently solve yet another multifacility location problem by using Algorithm 3.

Example 5.5

Consider the dataset $\mathbf{X}$ in $\mathbb{R}^{2}$ that consists of the latitudes and longitudes of $988$ US cities [26] with

[TABLE]

The optimal centers illustrated in Figure 5.5 are given by

[TABLE]

The total transportation cost in this problem is 5089.5150.

In the last example presented in this section we efficiently solve a higher dimensional multifacility location problem by using Algorithm 3 and compare it’s value of the cost function with the standard $k$-means algorithm.

Example 5.6

Let $\mathbf{X}$ in $\mathbb{R}^{13}$ be the wine dataset from the UCI Machine Learning Repositiory [28] consisting of $178$ demand points. We apply Algorithm 3 with

[TABLE]

The total costs using Algorithm 3 and the $k$-means algorithm are obtained in Table 4 showing that former algorithm is better than the latter.

6 Concluding Remarks

In this paper we develop a new algorithm to solve a class of multifacility location problems. Its implementation exhibits a better approximation compared to the classical $k$-means approach. This is demonstrated by a series of examples dealing with two-dimensional problems with nonnegative weights. Thus the verification and implementation of the proposed algorithm for real-life multifacility location problems in higher dimensions with arbitrary weights is a central direction of our future work. In addition, refining the initial cluster selection and the stopping criterion is an important area to explore.