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

0.57 MB·

by  Nguyen Mau Nam

NESTEROV’S SMOOTHING TECHNIQUE AND MINIMIZING
DIFFERENCES OF CONVEX FUNCTIONS FOR HIERARCHICAL
CLUSTERING
N. M. NAM1, W. GEREMEW2, S. REYNOLDS3 and T. TRAN4.
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
7
problem, a discrete optimization problem which can be shown to be NP-hard. Our approach is to
1
reformulate the problem as a continuous optimization problem by making some relaxations on the
0
2 discreteness conditions. Then Nesterov’s smoothing technique and a numerical algorithm for mini-
r mizing differences of convex functions called the DCA are applied to cope with the nonsmoothness
a
and nonconvexity of the problem. Numerical examples are provided to illustrate our method.
M
Key words. DC programming, Nesterov’s smoothing technique, hierarchical clustering, subgradi-

6 ent, Fenchel conjugate.
AMS subject classifications. 49J52, 49J53, 90C31
]
C
1 Introduction
O
.
h
Although convex optimization techniques and numerical algorithms have been the topics
t
a of extensive research for more than 50 years, solving large-scale optimization problems
m
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

2
sets where convexity is not assumed. One of the most successful approaches to go beyond
v
4 convexity is to consider the class of functions representable as differences of two convex
6
functions. Functions of this type are called DC functions, where DC stands for difference
4
4 of convex. It was recognized early by P. Hartman [10] that the class of DC functions has
0
many nice algebraic properties. For instance, this class of functions is closed under many
.
1
operations usually considered in optimization such as taking linear combination, maximum,
0
7 or product of a finite number of DC functions.
1
: Given a linear space X, a DC program is an optimization problem in which the objective
v
i functionf: X → Rcanberepresentedasf = g−h,whereg,h: X → Rareconvexfunctions.
X
This extension of convex programming enables us to take advantage of the available tools
r
a
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
1Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR
97207, USA ([email protected]). Research of this author was partly supported by the National
Science Foundation under grant #1411817.
2School of General Studies, Stockton University, Galloway, NJ 08205, USA
([email protected])
3Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR
97207, USA ([email protected] ).
4Fariborz Maseeh Department of Mathematics and Statistics, Portland State University, Portland, OR
97207, USA ([email protected]).
1
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 cluster-
ing 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 mathe-
matical 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
hierarchicalclusteringprobleminwhichthedistancemeasurementisdefinedbythesquared
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 hierar-
chical crusting problem: the DCA introduced by Pham Dinh Tao and Nesterov’s smoothing
technique.
We consider throughout the paper DC programming:
minimizef(x) := g(x)−h(x),x ∈ Rn, (2.1)
where g: Rn → R and h: Rn → 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: Rn → R, the Fenchel conjugate of g is defined by
g∗(y) := sup{(cid:104)y,x(cid:105)−g(x) | x ∈ Rn},y ∈ Rn.
Note that g∗: Rn → (−∞,∞] is also a convex function. In addition, x ∈ ∂g∗(y) if and only
if y ∈ ∂g(x), where ∂ denotes the subdifferential operator in the sense of convex analysis;
see, e.g., [19-21]
2
LetuspresentbelowtheDCAintroducedbyTaoandAn[3,9]asappliedto(2.1). Although
the algorithm is used for nonconvex optimization problems, the convexity of the functions
involved still plays a crucial role.
The DCA
INPUT: x ∈ Rn, N ∈ N.
0
for k = 1,...,N do
Find y ∈ ∂h(x ).
k k−1
Find x ∈ ∂g∗(y ).
k k
end for
OUTPUT: x .
N
Let us discuss below a convergence result of DC programming. A function h: Rn → R is
called γ-convex (γ ≥ 0) if the function defined by k(x) := h(x)− γ(cid:107)x(cid:107)2, x ∈ Rn, is convex.
2
If there exists γ > 0 such that h is γ−convex, then h is called strongly convex. We say that
an element x¯ ∈ Rn is a critical point of the function f defined by (2.1) if
∂g(x¯)∩∂h(x¯) (cid:54)= ∅.
Obviously, in the case where both g and h are differentiable, x¯ is a critical point of f if and
only if x¯ satisfies the Fermat rule ∇f(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 } generated
k
by the DCA. The following properties are valid:
(i) If g is γ -convex and h is γ -convex, then
1 2
γ +γ
f(x )−f(x ) ≥ 1 2(cid:107)x −x (cid:107)2 for all k ∈ N.
k k+1 k+1 k
2
(ii) The sequence {f(x )} is monotone decreasing.
k
(iii) If f is bounded from below, g is γ -convex and h is γ -convex with γ +γ > 0, and
1 2 1 2
{x } is bounded, then every subsequential limit of the sequence {x } is a critical point of f.
k k
Now we present a direct consequence of Nesterov’s smoothing technique given in [16]. In
the proposition below, d(x;Ω) denotes the Euclidean distance and P(x;Ω) denotes the
Euclidean projection from a point x to a nonempty closed convex set Ω in Rn.
Proposition 2.2 Given any a ∈ Rn and µ > 0, a Nesterov smoothing approximation of
ϕ(x) := (cid:107)x−a(cid:107) defined in Rn has the representation
ϕ (x) := 1 (cid:107)x−a(cid:107)2− µ(cid:2)d(x−a;B)(cid:3)2. (2.2)
µ
2µ 2 µ
Moreover, ∇ϕ (x) = P(x−a;B) and
µ µ
µ
ϕ (x) ≤ ϕ(x) ≤ ϕ (x)+ ,
µ µ
2
where B is the closed unit ball of Rn.
3
3 The Bilevel Hierarchical Clustering Problem
Given a set of m points (nodes) a1,a2,...,am in Rn, our goal is to decompose this set into
k clusters. In each cluster, we would like to find a point xi 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 a1,a2,...,am, 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 x1,x2,...,xk that could be anywhere in Rn. Then we define the total
center as the centroid of x1,x2,...,xk given by
1
x∗ := (x1+x2+···+xk).
k
The total cost of the tree is given by
m k
(cid:88) (cid:88)
ϕ(x1,...,xk) := min (cid:107)x(cid:96)−ai(cid:107)+ (cid:107)x(cid:96)−x∗(cid:107).
(cid:96)=1,...,k
i=1 (cid:96)=1
Note that here each ai 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
k
(cid:88)
φ(x1,...,xk) := min (cid:107)x(cid:96)−ai(cid:107),x1,...,xk ∈ Rn.
i=1,...,m
(cid:96)=1
Observe that φ(x1,...,xk) = 0 if and only if for every (cid:96) = 1,...,k, there exists i ∈
{1,...,m} such that x(cid:96) = ai, which means that x(cid:96) is a real node. Therefore, we consider
the constrained minimization problem:
minimize ϕ(x1,...,xk) subject to φ(x1,...,xk) = 0.
This problem can be converted to an unconstrained minimization problem:
minimize f (x1,...,xk) := ϕ(x1,...,xk)+λφ(x1,...,xk),x1,...,xk ∈ Rn, (3.1)
λ
where λ > 0 is a penalty parameter. Similar to the situation with the clustering prob-
lem, 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
4
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 ad-
dressed 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 ϕ and φ in (3.1) belong to the class of DC functions with the
following DC decompositions:
m k k k
(cid:88)(cid:104)(cid:88) (cid:88) (cid:105) (cid:88)
ϕ(x1,...,xk) = (cid:107)x(cid:96)−ai(cid:107)− max (cid:107)x(cid:96)−ai(cid:107) + (cid:107)x(cid:96)−x∗(cid:107)
r=1,...,k
i=1 (cid:96)=1 (cid:96)=1,(cid:96)(cid:54)=r (cid:96)=1
m k k m k
(cid:104)(cid:88)(cid:88) (cid:88) (cid:105) (cid:88) (cid:88)
= (cid:107)x(cid:96)−ai(cid:107)+ (cid:107)x(cid:96)−x∗(cid:107) − max (cid:107)x(cid:96)−ai(cid:107),
r=1,...,k
i=1 (cid:96)=1 (cid:96)=1 i=1 (cid:96)=1,(cid:96)(cid:54)=r
k m k m
(cid:88)(cid:88) (cid:88) (cid:88)
φ(x1,...,xk) = (cid:107)x(cid:96)−ai(cid:107)− max (cid:107)x(cid:96)−ai(cid:107).
s=1,...,m
(cid:96)=1 i=1 (cid:96)=1 i=1,i(cid:54)=s
It follows that the objective function f in (3.1) has the DC decomposition:
λ
k m k
(cid:104) (cid:88)(cid:88) (cid:88) (cid:105)
f (x1,...,xk) = (1+λ) (cid:107)x(cid:96)−ai(cid:107)+ (cid:107)x(cid:96)−x∗(cid:107)
λ
(cid:96)=1 i=1 (cid:96)=1
m k k m
(cid:104)(cid:88) (cid:88) (cid:88) (cid:88) (cid:105)
− max (cid:107)x(cid:96)−ai(cid:107)+λ max (cid:107)x(cid:96)−ai(cid:107) .
r=1,...,k s=1,...,m
i=1 (cid:96)=1,(cid:96)(cid:54)=r (cid:96)=1 i=1,i(cid:54)=s
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 approx-
imate the objective function f by a new DC function favorable for applying the DCA. To
λ
accomplish this goal, we simply replace each term of the form (cid:107)x−a(cid:107) from the first part
of f (x1,...,xk) (the positive part) by the smooth approximation (2.2), while keeping the
λ
second part (the negative part) the same. As a result, we obtain
(1+λ)µ (cid:88)m (cid:88)k (cid:13)(cid:13)x(cid:96)−ai(cid:13)(cid:13)2 µ (cid:88)k (cid:13)(cid:13)x(cid:96)−x∗(cid:13)(cid:13)2
f (x1,...,xk) := (cid:13) (cid:13) + (cid:13) (cid:13)
λµ 2 (cid:13) µ (cid:13) 2 (cid:13) µ (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
i=1 (cid:96)=1 (cid:96)=1
(1+λ)µ (cid:88)m (cid:88)k (cid:20) (cid:18)x(cid:96)−ai (cid:19)(cid:21)2 µ (cid:88)k (cid:20) (cid:18)x(cid:96)−x∗ (cid:19)(cid:21)2
− d ;B − d ;B
2 µ 2 µ
i=1 (cid:96)=1 (cid:96)=1
m k k m
(cid:88) (cid:88) (cid:88) (cid:88)
− max (cid:107)x(cid:96)−ai(cid:107) − λ max (cid:107)x(cid:96)−ai(cid:107).
r=1,...,k s=1,...,m
i=1 (cid:96)=1,(cid:96)(cid:54)=r (cid:96)=1 i=1,i(cid:54)=s
The original bilevel hierarchical clustering problem now can be solved using a DC program:
minimize f (x1,...,xk) = g (x1,...,xk)−h (x1,...,xk), x1,...,xk ∈ Rn.
λµ λµ λµ
5
In this formulation, g and h are convex functions on (Rn)k defined by
λµ λµ
g (x1,...,xk) := g1 (x1,...,xk)+g2 (x1,...,xk),
λµ λµ λµ
h (x1,...,xk) := h1 (x1,...,xk)+h2 (x1,...,xk)+h3 (x1,...,xk)+h4 (x1,...,xk),
λµ λµ λµ λµ λµ
with their respective components defined as
m k k
1+λ (cid:88)(cid:88) 1 (cid:88)
g1 (x1,...,xk) := (cid:107)x(cid:96)−ai(cid:107)2, g2 (x1,...,xk) := (cid:107)x(cid:96)−x∗(cid:107)2,
λµ 2µ λµ 2µ
i=1 (cid:96)=1 (cid:96)=1
(1+λ)µ (cid:88)m (cid:88)k (cid:20) (cid:18)x(cid:96)−ai (cid:19)(cid:21)2 µ (cid:88)k (cid:20) (cid:18)x(cid:96)−x∗ (cid:19)(cid:21)2
h1 (x1,...,xk) := d ;B , h2 (x1,...,xk) := d ;B ,
λµ 2 µ λµ 2 µ
i=1 (cid:96)=1 (cid:96)=1
m k k m
(cid:88) (cid:88) (cid:88) (cid:88)
h3 (x1,...,xk) := max (cid:107)x(cid:96)−ai(cid:107), h4 (x1,...,xk) := λ max (cid:107)x(cid:96)−ai(cid:107).
λµ λµ
r=1,...,k s=1,...,m
i=1 (cid:96)=1,(cid:96)(cid:54)=r (cid:96)=1 i=1,i(cid:54)=s
To facilitate the gradient and subgradient calculations for the DCA, we introduce a data
matrix A and a variable matrix X. The data A is formed by putting each ai, i = 1,...,m,
in the ith row, i.e.,
 
a a a ... a
11 12 13 1n
a a a ... a 
 21 22 23 2n
A =  . . . . .
 .. .. .. .. 
 
a a a ... a
m1 m2 m3 mn
Similarly, if x1,...,xk are the k cluster centers, then the variable X is formed by putting
each x(cid:96), (cid:96) = 1,...,k, in the (cid:96)th row, i.e.,
 
x x x ... x
11 12 13 1n
x x x ... x 
 21 22 23 2n
X =  . . . . .
 .. .. .. .. 
 
x x x ... x
k1 k2 k3 kn
Then the variable matrix X of the optimization problem belongs to Rk×n, the linear space
of k by n real matrices equipped with the inner product (cid:104)X,Y(cid:105) := trace(XTY). The
Frobenius norm on Rk×n is defined by
(cid:118) (cid:118)
(cid:114) (cid:117) k (cid:117) k
(cid:13) (cid:13) (cid:68) (cid:69) (cid:117)(cid:88) (cid:117)(cid:88)
(cid:13)X(cid:13) := X,X = (cid:116) (cid:104)x(cid:96),x(cid:96)(cid:105) = (cid:116) (cid:107)x(cid:96)(cid:107)2.
F
(cid:96)=1 (cid:96)=1
Finally, we represent the average of the k cluster centers by x∗, i.e., x∗ := 1 (cid:80)k xj.
k j=1
Gradient and Subgradient Calculations for the DCA
Let us start by computing the gradient of
g (X) = g1 (X)+g2 (X).
λµ λµ λµ
6
Using the Frobenius norm, the function g1 can equivalently be written as
λµ
m k
1+λ (cid:88)(cid:88)
g1 (X) = (cid:107)x(cid:96)−ai(cid:107)2
λµ 2µ
i=1 (cid:96)=1
m k
1+λ (cid:88)(cid:88)(cid:104) (cid:105)
= (cid:107)x(cid:96)(cid:107)2−2(cid:104)x(cid:96),ai(cid:105)+(cid:107)ai(cid:107)2
2µ
i=1 (cid:96)=1
1+λ (cid:104) (cid:13) (cid:13)2 (cid:68) (cid:69) (cid:13) (cid:13)2(cid:105)
= 2µ m(cid:13)X(cid:13)F −2 X,EkmA +k(cid:13)A(cid:13)F ,
where E is a k ×m matrix whose entries are all ones. Hence, one can see that g1 is
km λµ
differentiable and its gradient is given by
1+λ
∇g1 (X) = [mX−E A].
λµ µ km
Similarly, g2 can equivalently be written as
λµ
k
1 (cid:88)
g2 (X) = (cid:107)x(cid:96)−x∗(cid:107)2
λµ 2µ
(cid:96)=1
k
1 (cid:88)(cid:104) (cid:105)
= (cid:107)x(cid:96)(cid:107)2−2(cid:104)x(cid:96),x∗(cid:105)+(cid:107)x∗(cid:107)2
2µ
(cid:96)=1
(cid:20) (cid:21)
1 (cid:13) (cid:13)2 2(cid:68) (cid:69) 1(cid:68) (cid:69)
= 2µ (cid:13)X(cid:13)F − k X,EkkX + k X,EkkX
(cid:20) (cid:21)
1 (cid:13) (cid:13)2 1(cid:68) (cid:69)
= 2µ (cid:13)X(cid:13)F − k X,EkkX ,
where E is a k×k matrix whose entries are all ones. Hence, g2 is differentiable and its
kk λµ
gradient is given by
(cid:20) (cid:21)
1 1
∇g2 (X) = X− E X .
λµ µ k kk
Since g (X) = g1 (X)+g2 (X), its gradient can be computed by
λµ λµ λµ
∇g (X) = ∇g1 (X)+∇g2 (X)
λµ λµ λµ
(cid:20) (cid:21)
1+λ 1 1
= [mX−E A]+ X− E X
km kk
µ µ k
(cid:20) (cid:21)
1 1
= (1+λ)mX−(1+λ)E A+X− E X
km kk
µ k
(cid:20)(cid:20) (cid:21) (cid:21)
1 1
= [(1+λ)m+1]I − E X−(1+λ)E A .
kk kk km
µ k
Therefore,
1(cid:104)(cid:0) (cid:1) (cid:105) 1
∇g (X) = ((1+λ)m+1)I −J X−(1+λ)S , where J = E , and S = E A.
λµ kk kk km
µ k
7
Our goal now is to compute ∇g∗(Y), which can be accomplished by the relation
X = ∇g∗(Y) if and only if Y = ∇g(X).
The latter can equivalently be written as
[[1+(1+λ)m]I −J]X = [(1+λ)S+µY].
kk
Then with some algebraic manipulation we can show that
(cid:20) (cid:21)
1 1
∇g∗(Y) = X = I + J [(1+λ)S+µY].
kk
1+(1+λ)m [1+(1+λ)m](1+λ)m
Next, we will demonstrate in more detail the techniques we used in finding a subgradient
for the convex function h . Recall that h is defined by
λµ λµ
4
(cid:88)
h (X) = hi (X).
λµ λµ
i=1
We will start with the function h1 given by
λµ
(1+λ)µ (cid:88)m (cid:88)k (cid:20) (cid:18)x(cid:96)−ai (cid:19)(cid:21)2
h1 (X) = d ;B .
λµ 2 µ
i=1 (cid:96)=1
From its representation, one can see that h1 is differentiable, and hence its subgradient
λµ
coincides with its gradient, that can be computed by the partial derivatives with respect to
x1,··· ,xk, i.e.,
∂h1λµ(X) = (1+λ)(cid:88)m (cid:20)x(cid:96)−ai −P (cid:18)x(cid:96)−ai;B(cid:19)(cid:21).
∂x(cid:96) µ µ
i=1
∂h1
Thus, for (cid:96) = 1,2,...,k, ∇h1 (X) is the k×n matrix U whose (cid:96)th row is λµ(X).
λµ ∂x(cid:96)
Similarly, one can see that the function h2 given by
λµ
µ (cid:88)k (cid:20) (cid:18)x(cid:96)−x∗ (cid:19)(cid:21)2
h2 (X) = d ;B
λµ 2 µ
(cid:96)=1
is differentiable with its partial derivatives computed by
∂h2µ(X) = (cid:20)x(cid:96)−x∗ −P (cid:18)x(cid:96)−x∗;B(cid:19)(cid:21) − 1 (cid:88)k (cid:20)xj −x∗ −P (cid:18)xj −x∗;B(cid:19)(cid:21).
∂x(cid:96) µ µ k µ µ
j=1
Hence, for (cid:96) = 1,2,...,k, ∇h2(X) is the k×n matrix V whose (cid:96)th row is ∂h2µ(X).
µ ∂x(cid:96)
Unlike h1 and h2 , the convex functions h3 and h4 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 h3 as an example. We have
λµ
m k m
(cid:88) (cid:88) (cid:88)
h3 (X) = max (cid:107)x(cid:96)−ai(cid:107) = γ (X),
λµ i
r=1,...,k
i=1 (cid:96)=1,(cid:96)(cid:54)=r i=1
8
where, for i = 1,...,m,
 
k
 (cid:88) 
γ (X) := max γ (X) = (cid:107)x(cid:96)−ai(cid:107), r = 1,...,k .
i ir
 
(cid:96)=1,(cid:96)(cid:54)=r
Then, for each i = 1,...,m, we find W ∈ ∂γ (X) according to the subdifferential rule for
i i
the maximum of convex functions. Then define W := (cid:80)m W to get a subgradient of the
i=1 i
function h3 at X by the subdifferential sum rule. To accomplish this goal, we first choose
λµ
an index r∗ from the index set {1,...,k} such that
k
(cid:88)
γ (X) = γ (X) = (cid:107)x(cid:96)−ai(cid:107).
i ir∗
(cid:96)=1,(cid:96)(cid:54)=r∗
Using the familiar subdifferential formula of the Euclidean norm function, the (cid:96)th row w(cid:96)
i
for (cid:96) (cid:54)= r∗ of the matrix W is determined as follows
i

 x(cid:96)−ai if x(cid:96) (cid:54)= ai,
w(cid:96) := (cid:107)x(cid:96)−ai(cid:107)2
i
 u ∈ B if x(cid:96) = ai.
The r∗th row of the matrix W is wr∗ := 0.
i i
The procedure for calculating a subgradient of the function h4 given by
λµ
k m
(cid:88) (cid:88)
h4 (x1,...,xk) = λ max (cid:107)x(cid:96)−ai(cid:107),
λµ
s=1,...,m
(cid:96)=1 i=1,i(cid:54)=s
is very similar to what we just demonstrated for h3 .
λµ
At this point, we are ready to give a new DCA based algorithm for Model I.
Algorithm 1.
INPUT: X ,λ ,µ ,N ∈ N.
0 0 0
while stopping criteria (λ, µ) = false do
for k = 1,2,3,··· ,N do
Find Y ∈ ∂h (X ).
k 1λµ k−1
Find X ∈ ∂g∗ (Y ).
k 1λµ k
end for
update λ and µ.
end while
OUTPUT: X .
N
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 xk+1 to denote the total center. At
9
first we allow the total center xk+1 to be a free point in Rn, the same as the k cluster
centers. Then the total cost of the tree is given by
m k
(cid:88) (cid:88)
ϕ(x1,...,xk+1) := min (cid:107)x(cid:96)−ai(cid:107)+ (cid:107)x(cid:96)−xk+1(cid:107),x1,...,xk+1 ∈ Rn.
(cid:96)=1,...,k
i=1 (cid:96)=1
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
k+1
(cid:88)
φ(x1,...,xk+1) := min (cid:107)x(cid:96)−ai(cid:107) = 0.
i=1,...,m
(cid:96)=1
Our goal is to solve the optimization problem
minimize ϕ(x1,...,xk+1)
subject to φ(x1,...,xk+1),x1,...,xk+1 ∈ Rn.
Similar to the first model, this problem formulation can be converted to an unconstrained
minimization problem involving a penalty parameter λ > 0:
minimize f (x1,...,xk+1) := ϕ(x1,...,xk+1)+λφ(x1,...,xk+1),x1,...,xk+1 ∈ Rn. (3.2)
λ
Next, we apply Nesterov’s smoothing technique to get an approximation of the objective
function f given in (3.2) which involves two parameter λ > 0 and µ > 0:
m k+1 k m
1+λ (cid:88)(cid:88) 1 (cid:88) 1 (cid:88)
f (X) := (cid:107)x(cid:96)−ai(cid:107)2 + (cid:107)x(cid:96)−xk+1(cid:107)2 − (cid:107)xk+1−ai(cid:107)2
λµ
2µ 2µ 2µ
i=1 (cid:96)=1 (cid:96)=1 i=1
λµ (cid:88)m (cid:20) (cid:18)xk+1−ai (cid:19)(cid:21)2 (1+λ)µ (cid:88)m (cid:88)k (cid:20) (cid:18)x(cid:96)−ai (cid:19)(cid:21)2
− d ;B − d ;B
2 µ 2 µ
i=1 i=1 (cid:96)=1
µ (cid:88)k (cid:20) (cid:18)x(cid:96)−xk+1 (cid:19)(cid:21)2 (cid:88)m (cid:88)k (cid:88)k+1 (cid:88)m
− d ;B − max (cid:107)x(cid:96)−ai(cid:107) − λ max (cid:107)x(cid:96)−ai(cid:107).
2 µ r=1,...,k s=1,...,m
(cid:96)=1 i=1 (cid:96)=1,(cid:96)(cid:54)=r (cid:96)=1 i=1,i(cid:54)=s
As we will show in what follows, it is convenient to apply the DCA to minimize the function
f . This function can be represented as the differences of two convex functions defined on
λµ
R(k+1)×n using a variable X whose ith row is xi for i = 1,...,k+1:
f (X) = g (X)−h (X), X ∈ R(k+1)×n.
λµ λµ λµ
In this formulation, g and h are convex functions defined on R(k+1)×n by
λµ λµ
g (X) := g1 (X)+g2 (X)
λµ λµ λµ
and
h (X) := h1 (X)+h2 (X)+h3 (X)+h4 (X)+h5 (X)+h6 (X),
λµ λµ λµ λµ λµ λµ λµ
10