Nesterov's Smoothing Technique and Minimizing Differences of Convex Functions for Hierarchical Clustering

TL;DR

This paper reformulates a complex bilevel hierarchical clustering problem as a continuous optimization task, applying Nesterov's smoothing and DCA algorithms to effectively handle nonsmoothness and nonconvexity, demonstrated through numerical examples.

Contribution

It introduces a novel continuous reformulation of NP-hard hierarchical clustering models and applies advanced smoothing and optimization techniques to solve them.

Findings

Successful application of Nesterov's smoothing to hierarchical clustering

Effective use of DCA for nonsmooth nonconvex optimization

Numerical results demonstrate method's viability

Abstract

A bilevel hierarchical clustering model is commonly used in designing optimal multicast networks. In this paper, we consider two different formulations of the bilevel hierarchical clustering problem, a discrete optimization problem which can be shown to be NP-hard. Our approach is to reformulate the problem as a continuous optimization problem by making some relaxations on the discreteness conditions. Then Nesterov's smoothing technique and a numerical algorithm for minimizing differences of convex functions called the DCA are applied to cope with the nonsmoothness and nonconvexity of the problem. Numerical examples are provided to illustrate our method.

Equations

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

$$g^{\ast }(y):=\sup \{⟨y,x⟩-g(x)|x\in {\mathbb{R}}^n\},y\in {\mathbb{R}}^n.$$

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

$$f(x_k)-f(x_{k+1})\ge \frac{{\gamma }_1+{\gamma }_2}{2}\Vert x_{k+1}-x_k{\Vert }^2\text{mbox}\mathrm{forall}k\in \mathbb{N}.$$

$$\varphi_{\mu}(x):=\frac{1}{2\mu}\|x-a\|^{2}-\frac{\mu}{2}\big{[}d(\frac{x-a}{\mu};\mathbb{B})\big{]}^{2}.$$

$${\phi }_{\mu }(x)\le \phi (x)\le {\phi }_{\mu }(x)+\frac{\mu }{2},$$

$$x^{\ast }:=\frac{1}{k}(x^1+x^2+\cdots +x^k).$$

$$\phi (x^1,\dots ,x^k):=\sum _{i=1}^m\underset{\ell =1,\dots ,k}{\min }\Vert x^{\ell }-a^i\Vert +\sum _{\ell =1}^k\Vert x^{\ell }-x^{\ast }\Vert .$$

$$\varphi (x^1,\dots ,x^k):=\sum _{\ell =1}^k\underset{i=1,\dots ,m}{\min }\Vert x^{\ell }-a^i\Vert ,x^1,\dots ,x^k\in {\mathbb{R}}^n.$$

$$\text{mbox}\mathrm{minimize}\phi (x^1,\dots ,x^k)\text{mbox}\mathrm{subjectto}\varphi (x^1,\dots ,x^k)=0.$$

$$\text{mbox}\mathrm{minimize}{f}_{\lambda }(x^1,\dots ,x^k):=\phi (x^1,\dots ,x^k)+\lambda \varphi (x^1,\dots ,x^k),x^1,\dots ,x^k\in {\mathbb{R}}^n,$$

$$\phi (x^1,\dots ,x^k)$$

$$\displaystyle=\Big{[}\sum_{i=1}^{m}\sum_{\ell=1}^{k}\|x^{\ell}-a^{i}\|+\sum_{\ell=1}^{k}\|x^{\ell}-x^{*}\|\Big{]}-\sum_{i=1}^{m}\max_{r=1,\ldots,k}\sum_{\ell=1,\ell\neq r}^{k}\|x^{\ell}-a^{i}\|,$$

$$\varphi (x^1,\dots ,x^k)$$

$$f_{\lambda }(x^1,\dots ,x^k)$$

$$\displaystyle-\Big{[}\sum_{i=1}^{m}\max_{r=1,\ldots,k}\sum_{\ell=1,\ell\neq r}^{k}\|x^{\ell}-a^{i}\|+\lambda\sum_{\ell=1}^{k}\max_{s=1,\ldots,m}\sum_{i=1,i\neq s}^{m}\|x^{\ell}-a^{i}\|\Big{]}.$$

$$f_{\lambda \mu }(x^1,\dots ,x^k)$$

$$-\frac{(1+\lambda )\mu }{2}\sum _{i=1}^m\sum _{\ell =1}^k{\left[d\left(\frac{x^{\ell }-a^i}{\mu };\mathbb{B}\right)\right]}^2-\frac{\mu }{2}\sum _{\ell =1}^k{\left[d\left(\frac{x^{\ell }-x^{\ast }}{\mu };\mathbb{B}\right)\right]}^2$$

$$-\sum _{i=1}^m\underset{r=1,\dots ,k}{\max }\sum _{\ell =1,\ell \ne r}^k\Vert x^{\ell }-a^i\Vert -\lambda \sum _{\ell =1}^k\underset{s=1,\dots ,m}{\max }\sum _{i=1,i\ne s}^m\Vert x^{\ell }-a^i\Vert .$$

$$\text{mbox}\mathrm{minimize}{f}_{\lambda \mu }(x^1,\dots ,x^k)=g_{\lambda \mu }(x^1,\dots ,x^k)-h_{\lambda \mu }(x^1,\dots ,x^k),x^1,\dots ,x^k\in {\mathbb{R}}^n.$$

$$g_{\lambda \mu }(x^1,\dots ,x^k):=g_{\lambda \mu }^1(x^1,\dots ,x^k)+g_{\lambda \mu }^2(x^1,\dots ,x^k),$$

$$h_{\lambda \mu }(x^1,\dots ,x^k):=h_{\lambda \mu }^1(x^1,\dots ,x^k)+h_{\lambda \mu }^2(x^1,\dots ,x^k)+h_{\lambda \mu }^3(x^1,\dots ,x^k)+h_{\lambda \mu }^4(x^1,\dots ,x^k),$$

$$g_{\lambda \mu }^1(x^1,\dots ,x^k):=\frac{1+\lambda }{2\mu }\sum _{i=1}^m\sum _{\ell =1}^k\Vert x^{\ell }-a^i{\Vert }^2,g_{\lambda \mu }^2(x^1,\dots ,x^k):=\frac{1}{2\mu }\sum _{\ell =1}^k\Vert x^{\ell }-x^{\ast }{\Vert }^2,$$

$$h_{\lambda \mu }^1(x^1,\dots ,x^k):=\frac{(1+\lambda )\mu }{2}\sum _{i=1}^m\sum _{\ell =1}^k{\left[d\left(\frac{x^{\ell }-a^i}{\mu };\mathbb{B}\right)\right]}^2,h_{\lambda \mu }^2(x^1,\dots ,x^k):=\frac{\mu }{2}\sum _{\ell =1}^k{\left[d\left(\frac{x^{\ell }-x^{\ast }}{\mu };\mathbb{B}\right)\right]}^2,$$

$$h_{\lambda \mu }^3(x^1,\dots ,x^k):=\sum _{i=1}^m\underset{r=1,\dots ,k}{\max }\sum _{\ell =1,\ell \ne r}^k\Vert x^{\ell }-a^i\Vert ,h_{\lambda \mu }^4(x^1,\dots ,x^k):=\lambda \sum _{\ell =1}^k\underset{s=1,\dots ,m}{\max }\sum _{i=1,i\ne s}^m\Vert x^{\ell }-a^i\Vert .$$

$$\mathbf{A}=\left(\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \dots & a_{1n}\\ a_{21}& a_{22}& a_{23}& \dots & a_{2n}\\ ⋮& ⋮& ⋮& & ⋮\\ a_{m1}& a_{m2}& a_{m3}& \dots & a_{mn}\end{array}\right).$$

$$\mathbf{X}=\left(\begin{array}{ccccc}{x}_{11}& x_{12}& x_{13}& \dots & x_{1n}\\ x_{21}& x_{22}& x_{23}& \dots & x_{2n}\\ ⋮& ⋮& ⋮& & ⋮\\ x_{k1}& x_{k2}& x_{k3}& \dots & x_{kn}\end{array}\right).$$

$$\big{\|}\mathbf{X}\big{\|}_{F}:=\sqrt{\Big{\langle}\mathbf{X},\mathbf{X}\Big{\rangle}}=\sqrt{\sum_{\ell=1}^{k}\langle x^{\ell},x^{\ell}\rangle}=\sqrt{\sum_{\ell=1}^{k}\|x^{\ell}\|^{2}}.$$

$$g_{\lambda \mu }(\mathbf{X})=g_{\lambda \mu }^1(\mathbf{X})+g_{\lambda \mu }^2(\mathbf{X}).$$

$$g_{\lambda \mu }^1(\mathbf{X})$$

$$=\frac{1+\lambda }{2\mu }\sum _{i=1}^m\sum _{\ell =1}^k\left[\Vert x^{\ell }{\Vert }^2-2⟨x^{\ell },a^i⟩+\Vert a^i{\Vert }^2\right]$$

$$\displaystyle=\;\frac{1+\lambda}{2\mu}\left[m\big{\|}\mathbf{X}\big{\|}_{F}^{2}-2\Big{\langle}\mathbf{X},\mathbf{E}_{km}\mathbf{A}\Big{\rangle}+k\big{\|}\mathbf{A}\big{\|}_{F}^{2}\right],$$

$$\nabla g_{\lambda \mu }^1(\mathbf{X})=\frac{1+\lambda }{\mu }\left[m\mathbf{X}-{\mathbf{E}}_{km}\mathbf{A}\right].$$

$$g_{\lambda \mu }^2(\mathbf{X})$$

$$=\frac{1}{2\mu }\sum _{\ell =1}^k\left[\Vert x^{\ell }{\Vert }^2-2⟨x^{\ell },x^{\ast }⟩+\Vert x^{\ast }{\Vert }^2\right]$$

$$\displaystyle=\;\frac{1}{2\mu}\left[\big{\|}\mathbf{X}\big{\|}_{F}^{2}-\frac{2}{k}\Big{\langle}\mathbf{X},\mathbf{E}_{kk}\mathbf{X}\Big{\rangle}+\frac{1}{k}\Big{\langle}\mathbf{X},\mathbf{E}_{kk}\mathbf{X}\Big{\rangle}\right]$$

Full text

**NESTEROV’S SMOOTHING TECHNIQUE AND MINIMIZING DIFFERENCES OF CONVEX FUNCTIONS FOR HIERARCHICAL CLUSTERING

** N. M. NAM111Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR 97207, USA (mau.nam.nguyen@pdx.edu). Research of this author was partly supported by the National Science Foundation under grant #1411817., W. GEREMEW222School of General Studies, Stockton University, Galloway, NJ 08205, USA (wondi.geremew@stockton.edu), S. REYNOLDS333Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR 97207, USA (ser6@pdx.edu ). and T. TRAN444Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR 97207, USA (tuyen2@pdx.edu)..

**Abstract. A bilevel hierarchical clustering model is commonly used in designing optimal multicast networks. In this paper, we consider two different formulations of the bilevel hierarchical clustering problem, a discrete optimization problem which can be shown to be NP-hard. Our approach is to reformulate the problem as a continuous optimization problem by making some relaxations on the discreteness conditions. Then Nesterov’s smoothing technique and a numerical algorithm for minimizing differences of convex functions called the DCA are applied to cope with the nonsmoothness and nonconvexity of the problem. Numerical examples are provided to illustrate our method.

Key words. DC programming, Nesterov’s smoothing technique, hierarchical clustering, subgradient, Fenchel conjugate.

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

1 Introduction

Although convex optimization techniques and numerical algorithms have been the topics of extensive research for more than 50 years, solving large-scale optimization problems without the presence of convexity remains a challenge. This is a motivation to search for new optimization methods that are capable of handling broader classes of functions and sets where convexity is not assumed. One of the most successful approaches to go beyond convexity is to consider the class of functions representable as differences of two convex functions. Functions of this type are called DC functions, where DC stands for difference of convex. It was recognized early by P. Hartman [10] that the class of DC functions has many nice algebraic properties. For instance, this class of functions is closed under many operations usually considered in optimization such as taking linear combination, maximum, or product of a finite number of DC functions.

Given a linear space $X$, a DC program is an optimization problem in which the objective function $f\colon X\to\mathbb{R}$ can be represented as $f=g-h$, where $g,h\colon X\to\mathbb{R}$ are convex functions. This extension of convex programming enables us to take advantage of the available tools from convex analysis and optimization. At the same time, DC programming is sufficiently broad to use in solving many nonconvex optimization problems faced in recent applications. Another feature of DC programming is that it possesses a very nice duality theory; see [3] and the references therein. Although DC programming had been known to be important for many applications much earlier, the first algorithm for minimizing differences of convex functions called the DCA was introduced by Tao and An in [3, 9]. The DCA is a simple but effective optimization scheme used extensively in DC programming and its applications.

Cluster analysis or clustering is one of the most important problems in many fields such as machine learning, pattern recognition, image analysis, data compression, and computer graphics. Given a finite number of data points in a metric space, a centroid-based clustering problem seeks a finite number of cluster centers with each data point assigned to the nearest cluster center in a way that a certain distance, a measure of dissimilarity among data points, is minimized. Since many kinds of data encountered in practical applications have nested structures, they are required to use multilevel hierarchical clustering which involves grouping a data set into a hierarchy of clusters. In this paper, we apply the mathematical optimization approach to the bilevel hierarchical clustering problem. In fact, using mathematical optimization in clustering is a very promising approach to overcome many disadvantages of the $k-$mean algorithm commonly used in clustering; see [1, 2, 5, 15] and the references therein. In particular, the DCA was successfully applied in [2] to a bilevel hierarchical clustering problem in which the distance measurement is defined by the squared Euclidean distance. Although the DCA in [2] provides an effective way to solve the bilevel hierarchical clustering in high dimensions, it has not been used to solve the original model defined by the Euclidean distance measurement proposed in [5] which is not suitable for the resulting DCA according to the authors of [2]. By applying Nesterov’s smoothing technique and the DCA, we are able to solve the original model proposed in [5] in high dimensions.

The paper is organized as follows. In Section 2, we present basic definitions and tools of optimization that are used throughout the paper. In section 3, we study two models of bilevel hierarchical clustering problems, along with two new algorithms based on Nesterov’s smoothing technique and the DCA. Numerical examples and conclusions are presented in Section 4 and Section 5, respectively.

2 Basic Definitions and Tools of Optimization

In this section, we present two main tools of optimization used to solve the bilevel hierarchical crusting problem: the DCA introduced by Pham Dinh Tao and Nesterov’s smoothing technique.

We consider throughout the paper DC programming:

[TABLE]

where $g\colon\mathbb{R}^{n}\to\mathbb{R}$ and $h\colon\mathbb{R}^{n}\to\mathbb{R}$ are convex functions. The function $f$ in (2.1) is called a DC function and $g-h$ is called a DC decomposition of $f$.

Given a convex function $g\colon\mathbb{R}^{n}\to\mathbb{R}$, the Fenchel conjugate of $g$ is defined by

[TABLE]

Note that $g^{*}\colon\mathbb{R}^{n}\to(-\infty,\infty]$ is also a convex function. In addition, $x\in\partial g^{*}(y)$ if and only if $y\in\partial g(x)$, where $\partial$ denotes the subdifferential operator in the sense of convex analysis; see, e.g., [19-21]

Let us present below the DCA introduced by Tao and An [3, 9] as applied to (2.1). Although the algorithm is used for nonconvex optimization problems, the convexity of the functions involved still plays a crucial role.

The DCA

[TABLE]

Let us discuss below a convergence result of DC programming. A function $h\colon\mathbb{R}^{n}\to\mathbb{R}$ is called $\gamma$-convex ($\gamma\geq 0$) if the function defined by $k(x):=h(x)-\frac{\gamma}{2}\|x\|^{2}$, $x\in\mathbb{R}^{n}$, is convex. If there exists $\gamma>0$ such that $h$ is $\gamma-$convex, then $h$ is called strongly convex. We say that an element $\bar{x}\in\mathbb{R}^{n}$ is a critical point of the function $f$ defined by (2.1) if

[TABLE]

Obviously, in the case where both $g$ and $h$ are differentiable, $\bar{x}$ is a critical point of $f$ if and only if $\bar{x}$ satisfies the Fermat rule $\nabla f(\bar{x})=0$. The theorem below provides a convergence result for the DCA. It can be derived directly from [9, Theorem 3.7].

Theorem 2.1

*Consider the function $f$ defined by (2.1) and the sequence $\{x_{k}\}$ generated by the DCA. The following properties are valid:

(i) If $g$ is $\gamma_{1}$-convex and $h$ is $\gamma_{2}$-convex, then*

[TABLE]

(ii)* The sequence $\{f(x_{k})\}$ is monotone decreasing.

(iii) If $f$ is bounded from below, $g$ is $\gamma_{1}$-convex and $h$ is $\gamma_{2}$-convex with $\gamma_{1}+\gamma_{2}>0$, and $\{x_{k}\}$ is bounded, then every subsequential limit of the sequence $\{x_{k}\}$ is a critical point of $f$.*

Now we present a direct consequence of Nesterov’s smoothing technique given in [16]. In the proposition below, $d(x;\Omega)$ denotes the Euclidean distance and $P(x;\Omega)$ denotes the Euclidean projection from a point $x$ to a nonempty closed convex set $\Omega$ in $\mathbb{R}^{n}$.

Proposition 2.2

Given any $a\in\mathbb{R}^{n}$ and $\mu>0$, a Nesterov smoothing approximation of $\varphi(x):=\|x-a\|$ defined in $\mathbb{R}^{n}$ has the representation

[TABLE]

Moreover, $\nabla\varphi_{\mu}(x)=P(\frac{x-a}{\mu};\mathbb{B})$ and

[TABLE]

where $\mathbb{B}$ is the closed unit ball of $\mathbb{R}^{n}$.

3 The Bilevel Hierarchical Clustering Problem

Given a set of $m$ points (nodes) $a^{1},a^{2},\ldots,a^{m}$ in $\mathbb{R}^{n}$, our goal is to decompose this set into $k$ clusters. In each cluster, we would like to find a point $x^{i}$ among the nodes and assign it as the center for this cluster with all points in the cluster connected to this center. Then we will find a total center $x^{*}$ among the given points $a^{1},a^{2},\ldots,a^{m}$, and all centers are connected to this total center. The goal is to minimize the total transportation cost in this tree computed by the sum of the distances from the total center to each center and from each center to the nodes in each cluster. This is a discrete optimization problem which can be shown to be NP-hard. We will solve this problem based on continuous optimization techniques.

3.1 The Bilevel Hierarchical Clustering: Model I

The difficulty in solving this hierarchical clustering problem lies in the fact that the centers and total center have to be among the nodes. We first relax this condition with the use of artificial centers $x^{1},x^{2},\ldots,x^{k}$ that could be anywhere in $\mathbb{R}^{n}$. Then we define the total center as the centroid of $x^{1},x^{2},\ldots,x^{k}$ given by

[TABLE]

The total cost of the tree is given by

[TABLE]

Note that here each $a^{i}$ is assigned to its closest center. However, what we expect are the real centers, which can be approximated by trying to minimize the difference between the artificial centers and the real centers. To achieve this goal, define the function

[TABLE]

Observe that $\phi(x^{1},\ldots,x^{k})=0$ if and only if for every $\ell=1,\ldots,k$, there exists $i\in\{1,\ldots,m\}$ such that $x^{\ell}=a^{i}$, which means that $x^{\ell}$ is a real node. Therefore, we consider the constrained minimization problem:

[TABLE]

This problem can be converted to an unconstrained minimization problem:

[TABLE]

where $\lambda>0$ is a penalty parameter. Similar to the situation with the clustering problem, this new problem is nonsmooth and nonconvex, which can be solved by smoothing techniques and the DCA. Note that a particular case of this model in two dimensions was considered in [5] where the problem was solved using the derivative-free discrete gradient method established in [4], but this method is not suitable for large-scale settings in high dimensions. The DCA was used in [2] to solve a similar model with the squared Euclidean distance function used as the distance measurement. In this paper, the authors also addressed the difficulty of dealing with (3.1) using the DCA as not suitable for the resulting DCA. Nevertheless, we will show in what follows that the DCA is applicable to this model when combined with Nesterov’s smoothing technique.

Note that the functions $\varphi$ and $\phi$ in (3.1) belong to the class of DC functions with the following DC decompositions:

[TABLE]

It follows that the objective function $f_{\lambda}$ in (3.1) has the DC decomposition:

[TABLE]

This DC decomposition is not suitable for applying the DCA because there is no closed form for a subgradient of the function $g^{*}$ involved.

In the next step, we apply Nesterov’s smoothing technique from Proposition 2.2 to approximate the objective function $f_{\lambda}$ by a new DC function favorable for applying the DCA. To accomplish this goal, we simply replace each term of the form $\|x-a\|$ from the first part of $f_{\lambda}(x^{1},\ldots,x^{k})$ (the positive part) by the smooth approximation (2.2), while keeping the second part (the negative part) the same. As a result, we obtain

[TABLE]

The original bilevel hierarchical clustering problem now can be solved using a DC program:

[TABLE]

In this formulation, $g_{\lambda\mu}$ and $h_{\lambda\mu}$ are convex functions on $(\mathbb{R}^{n})^{k}$ defined by

[TABLE]

with their respective components defined as

[TABLE]

To facilitate the gradient and subgradient calculations for the DCA, we introduce a data matrix $\mathbf{A}$ and a variable matrix $\mathbf{X}$. The data $\mathbf{A}$ is formed by putting each $a^{i}$, $i=1,\ldots,m$, in the $i^{th}$ row, i.e.,

[TABLE]

Similarly, if $x^{1},\ldots,x^{k}$ are the $k$ cluster centers, then the variable $\mathbf{X}$ is formed by putting each $x^{\ell}$, $\ell=1,\ldots,k$, in the $\ell^{th}$ row, i.e.,

[TABLE]

Then the variable matrix $\mathbf{X}$ of the optimization problem belongs to $\mathbb{R}^{k\times n}$, the linear space of $k$ by $n$ real matrices equipped with the inner product $\langle\mathbf{X},\mathbf{Y}\rangle:=\mbox{\rm trace}(\mathbf{X}^{T}\mathbf{Y})$. The Frobenius norm on $\mathbb{R}^{k\times n}$ is defined by

[TABLE]

Finally, we represent the average of the $k$ cluster centers by $x^{*}$, i.e., $x^{*}:=\frac{1}{k}\sum_{j=1}^{k}x^{j}$.

Gradient and Subgradient Calculations for the DCA

Let us start by computing the gradient of

[TABLE]

Using the Frobenius norm, the function $g^{1}_{\lambda\mu}$ can equivalently be written as

[TABLE]

where $\mathbf{E}_{km}$ is a $k\times m$ matrix whose entries are all ones. Hence, one can see that $g^{1}_{\lambda\mu}$ is differentiable and its gradient is given by

[TABLE]

Similarly, $g^{2}_{\lambda\mu}$ can equivalently be written as

[TABLE]

where $\mathbf{E}_{kk}$ is a $k\times k$ matrix whose entries are all ones. Hence, $g^{2}_{\lambda\mu}$ is differentiable and its gradient is given by

[TABLE]

Since $g_{\lambda\mu}(\mathbf{X})=g^{1}_{\lambda\mu}(\mathbf{X})+g^{2}_{\lambda\mu}(\mathbf{X})$, its gradient can be computed by

[TABLE]

Therefore,

[TABLE]

Our goal now is to compute $\nabla g^{*}(Y)$, which can be accomplished by the relation

[TABLE]

The latter can equivalently be written as

[TABLE]

Then with some algebraic manipulation we can show that

[TABLE]

Next, we will demonstrate in more detail the techniques we used in finding a subgradient for the convex function $h_{\lambda\mu}$. Recall that $h_{\lambda\mu}$ is defined by

[TABLE]

We will start with the function $h^{1}_{\lambda\mu}$ given by

[TABLE]

From its representation, one can see that $h^{1}_{\lambda\mu}$ is differentiable, and hence its subgradient coincides with its gradient, that can be computed by the partial derivatives with respect to $x^{1},\cdots,x^{k}$, i.e.,

[TABLE]

Thus, for $\ell=1,2,\ldots,k$, $\nabla h^{1}_{\lambda\mu}(\mathbf{X})$ is the $k\times n$ matrix $\mathbf{U}$ whose $\ell^{th}$ row is $\frac{\partial h^{1}_{\lambda\mu}}{\partial x^{\ell}}(\mathbf{X})$.

Similarly, one can see that the function $h^{2}_{\lambda\mu}$ given by

[TABLE]

is differentiable with its partial derivatives computed by

[TABLE]

Hence, for $\ell=1,2,\ldots,k$, $\nabla h^{2}_{\mu}(\mathbf{X})$ is the $k\times n$ matrix $\mathbf{V}$ whose $\ell^{th}$ row is $\frac{\partial h^{2}_{\mu}}{\partial x^{\ell}}(\mathbf{X})$.

Unlike $h^{1}_{\lambda\mu}$ and $h^{2}_{\lambda\mu}$, the convex functions $h^{3}_{\lambda\mu}$ and $h^{4}_{\lambda\mu}$ are not differentiable, but both can be written as a finite sum of the maximum of a finite number of convex functions. Let us compute a subgradient of $h^{3}_{\lambda\mu}$ as an example. We have

[TABLE]

where, for $i=1,\ldots,m$,

[TABLE]

Then, for each $i=1,\ldots,m$, we find $\mathbf{W}_{i}\in\partial\gamma_{i}(\mathbf{X})$ according to the subdifferential rule for the maximum of convex functions. Then define $\mathbf{W}:=\sum_{i=1}^{m}\mathbf{W}_{i}$ to get a subgradient of the function $h^{3}_{\lambda\mu}$ at $\mathbf{X}$ by the subdifferential sum rule. To accomplish this goal, we first choose an index $r^{*}$ from the index set $\{1,\ldots,k\}$ such that

[TABLE]

Using the familiar subdifferential formula of the Euclidean norm function, the $\ell^{th}$ row $w_{i}^{\ell}$ for $\ell\neq r^{*}$ of the matrix $\mathbf{W}_{i}$ is determined as follows

[TABLE]

The $r^{*th}$ row of the matrix $\mathbf{W}_{i}$ is $w_{i}^{r^{*}}:=0$.

The procedure for calculating a subgradient of the function $h^{4}_{\lambda\mu}$ given by

[TABLE]

is very similar to what we just demonstrated for $h^{3}_{\lambda\mu}$.

At this point, we are ready to give a new DCA based algorithm for Model I.

Algorithm 1.

[TABLE]

3.2 The Bilevel Hierarchical Clustering: Model II

In this section, we introduce the second model to solve the bilevel hierarchical clustering problem. In this model, we use an additional variable $x^{k+1}$ to denote the total center. At first we allow the total center $x^{k+1}$ to be a free point in $\mathbb{R}^{n}$, the same as the $k$ cluster centers. Then the total cost of the tree is given by

[TABLE]

To force the $k+1$ centers to be chosen from the given nodes (or to make them as close to the nodes as possible), we set the constraint

[TABLE]

Our goal is to solve the optimization problem

[TABLE]

Similar to the first model, this problem formulation can be converted to an unconstrained minimization problem involving a penalty parameter $\lambda>0$:

[TABLE]

Next, we apply Nesterov’s smoothing technique to get an approximation of the objective function $f$ given in (3.2) which involves two parameter $\lambda>0$ and $\mu>0$:

[TABLE]

As we will show in what follows, it is convenient to apply the DCA to minimize the function $f_{\lambda\mu}$. This function can be represented as the differences of two convex functions defined on $\mathbb{R}^{(k+1)\times n}$ using a variable $X$ whose $i^{th}$ row is $x^{i}$ for $i=1,\ldots,k+1$:

[TABLE]

In this formulation, $g_{\lambda\mu}$ and $h_{\lambda\mu}$ are convex functions defined on $\mathbb{R}^{(k+1)\times n}$ by

[TABLE]

and

[TABLE]

with their respective components given by

[TABLE]

Gradient and Subgradient Calculations for the DCA

Let us start by computing the gradient of the first part of the DC decomposition, i.e.,

[TABLE]

By applying similar techniques used in computing gradients/subgradients for Model I, $\nabla g^{1}_{\lambda\mu}(\mathbf{X})$ can be written as

[TABLE]

where $\mathbf{E}$ is the $(k+1)\times m$ matrix whose entries are all ones. Similarly, $\nabla g^{2}_{\lambda\mu}(\mathbf{X})$ can also be written as

[TABLE]

where $\mathbf{I}$ is the $(k+1)\times(k+1)$ identity matrix, and $\mathbf{T}$ is the $(k+1)\times(k+1)$ matrix whose entries are all zeros except its $(k+1)^{th}$ row and $(k+1)^{th}$ column, which both are filled by the vector $(-1,-1,-1,\ldots,-1,k-1)$.

It follows that

[TABLE]

where $c_{1}=1+(1+\lambda)m$. Our goal now is to compute $\nabla g_{\lambda\mu}^{*}(\mathbf{Y})$, which can be accomplished by the relation

[TABLE]

The latter can be equivalently written as

[TABLE]

whose solutions can be explicitly computed by its $\ell^{th}$ row for $\ell=1,\ldots,k+1$:

[TABLE]

In the representation

[TABLE]

the convex functions $h^{1}_{\lambda\mu},\,h^{2}_{\lambda\mu},\,h^{3}_{\lambda\mu},\,\text{and}\,h^{4}_{\lambda\mu}$ are differentiable. The partial derivatives of $h^{1}_{\lambda\mu}$ are given by

[TABLE]

Then $\nabla h^{1}_{\lambda\mu}(\mathbf{X})$ is the $(k+1)\times n$ matrix $\mathbf{L}$ whose $\ell^{th}$ row is $\frac{\partial h^{1}_{\lambda\mu}}{\partial x^{\ell}}(\mathbf{X})$ for $\ell=1,\ldots,k+1$. Similarly, the partial derivatives of $h^{2}_{\lambda\mu}$ are given by

[TABLE]

Then $\nabla h^{2}_{\lambda\mu}(\mathbf{X})$ is the $(k+1)\times n$ matrix $\mathbf{M}$ whose $\ell^{th}$ row is $\frac{\partial h^{2}_{\lambda\mu}}{\partial x^{\ell}}(\mathbf{X})$, for $\ell=1,\ldots,k+1$. Now we compute the gradient of $h^{3}_{\lambda\mu}$. We have

[TABLE]

Thus, $\nabla h^{3}_{\lambda\mu}(\mathbf{X})$ is the $(k+1)\times n$ matrix $\mathbf{U}$ whose $l^{th}$ row is $\frac{\partial h^{3}_{\lambda\mu}}{\partial x^{\ell}}(\mathbf{X})$ for $\ell=1,\ldots,,k+1$.

Let us compute the gradient of $h^{4}_{\lambda\mu}$. We have

[TABLE]

Similarly, $\nabla h^{4}_{\lambda\mu}(\mathbf{X})$ is the $(k+1)\times n$ matrix $V$ whose $\ell^{th}$ row is filled with $\frac{\partial h^{4}_{\lambda\mu}}{\partial x^{\ell}}(\mathbf{X})$, for $\ell=1,\ldots,k+1$.

The procedure for computing subgradients of the last two nondifferentiable components, $h^{5}_{\lambda\mu}$ and $h^{6}_{\lambda\mu}$, is similar to the procedure we described for computing subgradients of $h^{3}_{\lambda\mu}$ in Model I.

The following is a DCA based algorithm for Model II.

Algorithm 2.

[TABLE]

4 Numerical Experiments

We use MATLAB to code our algorithms and perform numerical experiments on a MacBook Pro with 2.2 GHz Intel Core i7 Processor, and 16 GB 1600 MHz DDR3 Memory. For our numerical experiments, we use three data sets: one artificial data set with 18 data points in $\mathbb{R}^{2}$ (see Figure 1(a)), the EIL76 and the PR1002 from [17] (see Figure 1(b) and Figure 1(c)), respectively.

The two MATLAB codes used to implement our two algorithms have two major parts: an outer loop for updating the penalty and the smoothing parameters and an inner loop for updating the cluster centers. The penalty parameter $\lambda$ and the smoothing parameter $\mu$ are updated as follows. Choosing $\mu_{0}>0$, $\lambda_{0}>0$ and $\sigma_{1}>1$, $\sigma_{2}\in(0,1)$, we update $\lambda_{i+1}=\sigma_{1}\lambda_{i}$ and $\mu_{i+1}=\sigma_{2}\mu_{i}$ for $i\geq 0$ after each outer loop. To choose $\sigma_{1}$ and $\sigma_{2}$, we first let $N$ be the number of outer iterations and choose $\lambda_{0},\lambda_{\textup{max}}$ as the initial and final values of $\lambda$, respectively, and similarly for $\mu_{0},\mu_{\textup{min}}$. Then choose growth/decay parameters according to $\sigma_{1}=\left(\lambda_{\textup{max}}/\lambda_{0}\right)^{1/N}$ and $\sigma_{2}=\left(\mu_{\textup{min}}/\mu_{0}\right)^{1/N}$.

By trial and error we find that the values chosen for $\lambda_{0}$, $\lambda_{\textup{max}}$, $\mu_{0}$, and $\mu_{\textup{min}}$ in large part determine the performance of the two algorithms for each data set. Intuitively, we see that very large values of $\lambda$ will over-penalize the distance between an artificial center and its nearest data node and may prevent the algorithm from clustering properly. We therefore use $\lambda_{0}\leq 1\ll\lambda_{\textup{max}}$ so that the algorithm has a chance to cluster the data before the penalty parameter takes effect. Similarly, we choose $\mu_{\textup{min}}\ll 1\leq\mu_{0}$.

We select the starting center $\mathbf{X}_{0}$ at a certain radius, $\gamma\;\text{rad}(\mathbf{A})$, from the median point, $\text{median}(\mathbf{A})$, of the entire data set, i.e.,

[TABLE]

where $\text{rad}(\mathbf{A}):=\max\{\|a^{i}-\text{median}(\mathbf{A})\|\;|\;a^{i}\in\mathbf{A}\}$, $\gamma$ is a randomly chosen real number from the standard uniform distribution on the open interval (0,1), $\mathbf{U}$ is a $k\times n$ matrix whose $k$ rows are randomly generated unit vectors in $\mathbb{R}^{n}$, and the sum is in the sense of adding a vector to each row of a matrix.

As showed in Tables 1, 2, and 3, both algorithms identify the optimal solutions with reasonable amount of time for both DS18 and EIL76 with two and three cluster centers, respectively. To get a good starting point which yields a better estimate of the optimal value for bigger data sets such as PR1002, we use a method called radial search described as follows. Given initial radius $r_{0}>0$ and $m\in\mathbb{N}$, set $\gamma=ir_{0}$ for $i=1,\ldots,m$. Then we test the algorithm with different starting points given by $\mathbf{X}_{0}(i)=\text{median}(\mathbf{A})+ir_{0}(\text{rad}(\mathbf{A})\mathbf{U})$ for $i=1,\ldots,m$. Figure 2 shows the result of the method applied to PR1002 with six cluster centers, where the $y$-axis represents the optimal value returned by ALG1 with different starting points $\mathbf{X}_{0}(i)$, as represented on the $x$-axis.

For comparison purposes, both Cost1 and Cost2 are computed by the same way. First, we systematically reassign the $k$ cluster centers returned by the respective algorithms by $k$ real nodes that are close to them, i.e., for $\ell=1,\ldots,k$

[TABLE]

Then the total center $x^{*}$ will be a real node, from the remaining nodes, whose sum of distances from the $k$ reassigned centers is the minimal, i.e.,

[TABLE]

The total cost is computed by adding the distance of each real node to its closest center (including the total center), and the distances of the total center from the $k$ cluster centers, i.e.,

[TABLE]

5 Conclusions

In our numerical experiments, one of the major challenge we face is how to choose optimal parameters for our algorithms. We observe that parameter selection is the decisive factor in terms of accuracy and speed of convergence of our proposed algorithms. The performance of the proposed algorithms highly depends on the initial values set to the penalty and smoothing parameters, $\lambda_{0}$ and $\mu_{0}$; and their respective growth/decay factors $\sigma_{1}$ and $\sigma_{2}$. A better guidance of selecting the optimal parameter than the methods we suggest in this paper will be addressed in or future research.