Learning Strictly Orthogonal p-Order Nonnegative Laplacian Embedding via Smoothed Iterative Reweighted Method

1 Introduction

Data that reside in high-dimensional space are often intractable due to the computational complexity and the lack of intuition. In traditional Laplacian embedding (LE), the intrinsic subspace/manifold in high-dimensional space can be explored in such a way that the inherent data structures are well preserved and made more apparent due to the fact that the features less related to others will be pruned. LE is a powerful nonlinear graph based embedding method, which was first introduced as "quadratic placement" in 1970s [Hall, 1970]. Recently, the real power of LE was revealed as its relation to graph clustering [Hagen and Kahng, 1992; Chan et al., 1994; Shi and Malik, 2000]. The eigenvectors of the Laplacian matrix provide the approximation to the Ratio Cut spectral clustering [Chan et al., 1994] and it has been proved that LE and ratio cut clustering are mathematically identical [Luo et al., 2009].

∗Equal contribution.

However, both positive and negative values in the solution of multi-way clustering tasks make the results hard to interpret directly, because the clustering indicator vectors require nonnegative results. In two-way clustering, this is not a problem, because a linear Ψ-transformation [Ding and He, 2004] of the eigenvectors leads to two genuine indicator vectors (each row has only one nonzero positive entry). Thus, mixed-sign solution is a generic difficulty for multi-way spectral clustering. To solve this problem, a clustering step usually has to be performed after the embedding is learned. That is, in traditional way, the clustering indicator vectors approximated by the eigenvectors of the Laplacian matrix will be grouped by using K-means clustering [Hartigan and Wong, 1979] in the eigenvector space. Thus, the traditional clustering solution provided from this process is neither stable nor intuitive, which may also be very sensitive to data outliers. To tackle this difficulty, Nonnegative Laplacian Embedding (NLE) method [Luo et al., 2009] was proposed by additionally imposing the nonnegative constraints on the embeddings.

Despite the fact that the nonnegativity can be achieved in the NLE method, there are still some difficulties of this model that are not well addressed. It has been noted that the NLE method imposes the nonnegative constraint at the price of relaxing the orthogonality on the learned approximations [Ding et al., 2006], although the orthogonality constraint (XT X = I) is of significant importance to guarantee a good performance. The true meaning of the orthogonality constraint is to prevent degenerate solution (X → 0). For one dimensional problem, the orthogonality can avoid that the embedded data collapse into a point. For multi-dimensional problem, the orthogonality can prevent data points from collapsing into a subspace with dimensions less than desired.

In this paper, we propose a new approach to learn LE with strictly guaranteed orthogonality and nonnegativity in the solution. Unlike using the auxiliary function method [Lee and Seung, 2001] to derive the solution algorithm for NLE in [Luo et al., 2009], the orthogonality of our solution is rigorously achieved by using the Alternating Direction Method of Multipliers (ADMM) [Bertsekas, 1996; Boyd et al., 2011], leading to a more stable solution and a better performance in the problem of spectral clustering. We also keep the nonnegativity in the constraint, such that the clustering membership can be readily read off from the embedded data due to the nonnegative constraint, i.e., we can consider each row of the solution X as the posterior clustering probability. In other words, the values in i-th row of the solution can be viewed as the likelihoods that the i-th data point belongs to different clusters, which gives our new approach the soft clustering capability that is crucial in many real-world applications.

Finally, we recognize that the squared ℓ2-norm distance used in the traditional LE and NLE objectives does not guarantee the optimal embedding [Wang et al., 2015] and is also notoriously known to be sensitive to the outliers [Wang et al., 2012; Nie et al., 2013; Wang et al., 2013c; Nie et al., 2016; Liu et al., 2017]. With strict orthogonality and nonnegativity guaranteed simultaneously in the solution, we are also interested in promoting the robustness of our new NLE model by using the p-th order (0 < p ≤ 2) of the ℓ2-norm distance in the objective. As a result, the proposed optimization objective is a quadratic function with both orthonormal and nonnegative constraints, which is highly nonlinear and nonconvex in its feasible domain. The p-th order term further makes the objective nonsmooth and difficult to efficiently optimize in general. To solve this challenging optimization problem, we propose a novel smoothed iterative reweighted method. Compared to the iterative reweighted method proposed in [Candes et al., 2008; Nie et al., 2010] to solve the ℓ1-norm or ℓ2,1-norm minimization problems, our new optimization framework explicitly adds a smoothness term which can improve numerical stability. Most importantly, as an important theoretical contribution, we rigorously prove the convergence of our new iterative algorithm with the smoothness term, which, though, was not present in [Candes et al., 2008; Nie et al., 2010] and their following works.

To evaluate the proposed robust NLE objective that uses the p-th order of the ℓ2-norm distances and our new smoothed iterative reweighted method, we have performed extensive empirical studies. The promising experimental results have validated the effectiveness of our new methods.

2 Strictly Orthogonal p-Order Nonnegative Laplacian Embedding

Given a set of n data points, we can represent the pairwise similarities between these data points by a graph G = {V, E}, where the data points are represented by the vertices V and |V| = n. Suppose that W ∈ ℜn×n denotes the affinity matrix of the graph G where wij measures the similarity between the i-th and the j-th vertices, quadratic placement [Hall, 1970] aims to embed the vertices of the graph into the onedimensional space with coordinates (x1, . . . , xn), such that if the i-th and the j-th vertices are similar (i.e., wij is large), they should be adjacent in embedded space, i.e., (xi -xj)2 should be small. This can be achieved by the following objective [Hall, 1970]:

where x = [x1, . . . , xn]T , and D = diag (d1, d2, . . . , dn) is the degree matrix of the graph with di = ∑j wij .

The one-dimensional quadratic placement in Eq. (1) can be generalized to r-dimensional LE by minimizing the following objective [Luo et al., 2009]:

where X = [x1, x2, . . . , xn]T ∈ ℜn×r . Obviously, the i-th row of X, i.e., xTi ∈ ℜr , is the embedding of the i-th data point in the r-dimensional space. Here, the orthonormal constraint of XT X = I is imposed in Eq. (2) to avoid degenerate solutions, which is critical as analyzed in [Luo et al., 2009; Ding et al., 2006].

To decode the clustering membership from X in an easier way, Luo et al. [Luo et al., 2009] further developed LE by additionally imposing the nonnegative constraint onto the embedding matrix X by minimizing the following objective:

The squared ℓ2-norm distances used in the both objectives in Eq. (2) and Eq. (3) do not tolerate large value of distance, thus making the distances in the embedded space tend to be even, i.e., not too large but also not too small. Therefore, solving the objective in Eq. (2) or Eq. (3) may not find the optimal embedding such that most of the distances of local data pairs are minimized but a few of them are large [Wang et al., 2015]. Motivated by recent papers that use not-squared ℓ2-norm distances [Wang et al., 2012; Nie et al., 2013; Wang et al., 2013c; Wang et al., 2014; Nie et al., 2016; Liu et al., 2017] or the p-th order of the ℓ2-norm distances [Wang et al., 2011; Wang et al., 2013a; Wang et al., 2015] to promote the robustness of learning models, in this paper we propose to solve the following problem to find the optimal spectral embedding from an input graph:

where 0 < p ≤ 2. Obviously, the NLE method in Eq. (3) [Luo et al., 2009] is a special case of our new method when p = 2. More importantly, by setting p ≤ 1, the method will focus on minimizing most of the distances of local data pairs. Thus, we call Eq. (4) as the proposed strictly orthogonal pOrder Nonnegative Laplacian Embedding (PO-NLE) method.

3 Smoothed Iterative ReweightedMethod and its Convergence

Although the motivation of our new objective in Eq. (4) is clear, it is nonsmooth and difficult to efficiently solve in general. To solve this challenging optimization problem, in this section we will first introduce a novel smoothed iterative reweighted method.

First, let us consider a general problem as follows:

When gi(x) is a vector output function, tr ((gT i (x)gi(x))p2 ) becomes the following term:

Equation (5) is nonsmooth, thus we turn to solve the following smooth problem [Liu et al., 2017]:

where δ > 0 is a small constant. When δ → 0, Eq. (7) is reduced to Eq. (5) since the following equation holds:

Before deriving the algorithm for optimizing the problem in Eq. (7), we need the following lemmas. First, according to the chain rule in calculus, we have:

Lemma 1 Suppose g(x) is a matrix output function, h(x) is a scalar output function, x is a scalar, vector or matrix variable, then we have:

According to the chain rule in Lemma 1, we can easily derive the following two lemmas:

Lemma 2 Suppose g(x) is a scalar, vector or matrix output function, x is a scalar, vector or matrix variable, then we have:

Lemma 3 Suppose g(x) is a scalar, vector or matrix output function, x is a scalar, vector or matrix variable, D is a constant and D is symmetrical if D is a matrix, then we have:

Now we derive the algorithm to optimize the problem in Eq. (7). The Lagrangian function of the problem in Eq. (7) is:

where r(x, λ) is a Lagrangian term for the constraint x ∈ C. By setting the derivative of Eq.(12) w.r.t. x to zero, we have

According to Lemma 2, Eq.(13) can be rewritten as

Algorithm 1 The algorithm to solve the problem (7)

Initialize x ∈ C.

while not converge do

1. For each i, calculate

1. Update x by solving the problem:

end while

If we can find a solution x that satisfies the Eq.(14), we usually find a local or global optimal solution to the problem in Eq. (7) according to the Karush-Kuhn-Tucker (KKT) conditions [Boyd and Vandenberghe, 2004]. However, directly finding a solution x that satisfies Eq.(14) is not a easy task. In this paper, following [Candes et al., 2008; Nie et al., 2010] we propose an iterative algorithm to find it using the following observation: if

is a constant, Eq.(14) is reduced to:

which is equivalent to solving the following problem:

Based on the observation above, we can first guess a solution x. Then we calculate Di using the current solution of x and update x by the optimal solution of the problem in Eq. (17) by calculating Di. We iteratively perform this procedure until it converges, which is summarized in Algorithm 1.

The convergence of Algorithm 1 is guaranteed by the following theorem. The proof of Theorem 1 will be supplied in the extended journal version of this paper due to space limit.

Theorem 1 The Algorithm 1 will monotonically decrease the objective of the problem (7) in each iteration until the algorithm converges.

In the convergence, the equality in Eq. (14) will hold, thus the KKT condition of problem (7) is satisfied. Therefore, the Algorithm 1 will converge to a local optimum solution to the problem (7). If the problem (7) is convex, the Algorithm 1 will converge to a global optimum solution.

Here we note that the iterative reweighted method introduced in [Candes et al., 2008; Nie et al., 2010] solves the nonsmooth ℓ1-norm or ℓ2,1-norm minimization problems. However, the method described in [Candes et al., 2008; Nie et al., 2010] does not explicitly use the smoothness constant (i.e., δI in Eq. (7)). Without this smoothness term, the algorithm is heavily impacted by the singularity problem due to inverted matrices that divide 0s, which routinely leads to inferior learning performances. To improve the numerical stability, in [Nie et al., 2010] and the following works by the same group of authors [Wang et al., 2013b; Nie et al., 2013], a smoothness term was informally added for empirical purpose. But they only theoretically proved the convergence of the algorithm that does not use the smoothness term and did not provide any theoretical analysis on the objectives that use the smoothness term. As an important theoretical contribution of this paper, we formally introduce the smoothness term (i.e., δI in Eq. (7)) into our algorithm and theoretically prove the convergence of our algorithm in which the smoothness term leads to much more stable solutions. We call Algorithm 1 as the proposed Smoothed Iterative Reweighted Method, which can be broadly used to solve a variety of difficult machine learning problems that minimize the objectives using the p-th order of ℓ2-norm distances, the p-th order of ℓp-norm distances, or the p-th order of the Schatten p-norm.

4 Algorithm to Solve the Problem in Eq. (4)

Equipped with Algorithm 1, we can derive the solution algorithm to our new PO-NLE objective in Eq. (4) now. According to Step 2 of Algorithm 1, the key step to solve Eq. (4) is to solve the following problem:

Denote ˜Wij = wijdij and let ˜D be the diagonal matrix with the i-th diagonal entry as ∑j ˜wij = wijdij. The problem in Eq. (20) can be written as following:

Obviously, Eq. (21) is identical to the NLE objective in Eq. (3), which was proposed in [Luo et al., 2009]. In [Luo et al., 2009], a solution algorithm was derived using the auxiliary function method [Lee and Seung, 2001]. However, as analyzed in [Ding et al., 2005; Ding et al., 2006] the orthogonal constraint indeed are not guaranteed, which, though, is very important to avoid degenerate solutions [Ding et al., 2006]. Thus, instead of using the solution algorithm provided in [Luo et al., 2009], we derive the solution algorithm to solve Eq. (21) using the ADMMmethod.

Denoting L = ˜D - ˜W for brevity 1 , we can write the objective in Eq. (21) as following:

1 In practice, due to the zero mode of the Laplacian matrix of a graph [Wang et al., 2010], we compute L = ˜D- ˜W+ ˜W++ n2 eT e to ensure that L is positive definite, where ˜W++ = ∑i,j ˜W and e is the vector with all entries to be 1.

We can solve Eq. (22) by solving the following equivalent optimization problem:

where the constraint of XT X = I in Eq. (22) is implicitly enforced by the constraints of Y = X and YT Y = I.

Following the ADMM optimization framework, we need to solve the following optimization problem:

in which we introduced the Lagrangian multiplier Λ for the constraint of Y = X. The detailed procedures to solve Eq. (24) using the ADMM method will be supplied in the extended journal version of this paper due to space limit.

5 Experiment and Results

In this section we empirically evaluate our new PO-NLE method on one synthetic data set, four data sets from the UCI Machine Learning Data Repository, and three image data sets. We will compare our new method against its counterparts: NLE, Normalized Cut (NCut) [Shi and Malik, 2000] and Laplacian Embedding (LE).

In our evaluations, we use clustering accuracy and clustering purity to measure the performance of the compared methods. We also study the robustness of our method on the real world data sets when they are contaminated with noise. The performance variations when we increase the value of p will be shown to validate our hypothesis that the optimal solution is usually obtained when p is less than 2 and close to 1 (it depends on data sets), given that noises are present in the data. Orthogonality of the solution will be illustrated and compared against the NLE method in [Luo et al., 2009].

5.1 Experiments on a Synthetic Data Set

To illustrate the effectiveness of our new PO-NLE method, we create a synthetic data set as follows. We first randomly generate 3 data points as centroids in the 30-dimensional space. Then we generate 3 groups of data points and each group consists of 39 data points which are randomly distributed around one of the three centroids. A threshold is set to make the distance of groups large enough. As shown in Figure 1, different colors (red, black and blue) and shapes are used to represent different groups of data points. We randomly initialize X (0 ≤ X ≤ 1) and set ρ = 1.02, µ = 0.1 and p = 0.8 in our algorithm. we use K-Nearest Neighbors with heat kernel to construct our adjacency matrix W. The variation of the objective value when our algorithm iterates are shown as the red curve in Figure 1. For visualization purpose, we set r = 3, i.e., we embed the original data into the 3-dimensional space using our new PO-NLE algorithm. The x, y and z axes of the 3D plots in the figure correspond to the first, second and third row in matrix X, respectively.

From Figure 1, we observe that the objective function monotonically decreases in each iteration, which empirically confirms the convergence of the solution algorithm to solve PO-NLE derived by our new smoothed iterative reweighted method. Moreover, for each checkpoint shown by the black circles on the objective curve, the clustering structure of the experimental data becomes more and more clear in the 3D plots when the algorithm iterates. The three clusters of data points gradually find a solution to separate themselves apart and fall on different axes. Note that, due to the nonnegative constraints on X, data points will finally converge on the positive part of each axis. This observation clearly demonstrate the effectiveness of the proposed new method.

5.2 Studies of the Orthogonality of the Solutions of Our New Method

An important improvement of our new method over the NLE method [Luo et al., 2009] is that the orthogonality of our solution is rigorously guaranteed, which, as analyzed in [Ding et al., 2005; Ding et al., 2006] is very important to avoid degenerate solutions. Thus, in this subsection we empirically study the orthogonality of the solutions of our new method and compare them against the solutions obtained from the NLE method. Figure 2 compare the visualizations of XT X learned from the two compared methods on the Glass data set by the heatmaps. The heatmap of our new method is on the top and that of the NLE method is at the bottom. From the results we can see that the learned embeddings from our method are strictly orthogonal as shown in Figure 2(a), which will in return lead to better clustering performances and robustness after embedding. In contrast, the NLE method failed to guarantee the orthogonality, as can be seen in Figure 2(b).

5.3 Experiments on Noiseless Real Data Sets

Now we compare our new method, NLE, NCut and LE on the seven standard data sets as summarized in Table 1. Each data set will be tested by different algorithms independently for 200 times. For NCut and LE algorithms, we run K-means

(b) Visualization of XT X learned by the NLE method.

Figure 2: The comparison of orthogonality between our method and NLE on glass data set.

clustering with random initialization for 50 times and report the best results.

The performances of the compared methods are reported in the top half of Table 2, from which we can see that our method clearly outperforms all other competing methods, especially on those comparatively noisier data sets. Due to the nonnegative solutions of our new method, we do not need any additional clustering step. Instead, the clustering membership can be readily read off directly from the learned embeddings. The strictly guaranteed orthogonality constraint avoids degenerate solution and helps improve the performance compared with loosely constrained NLE method which does not have such desirable property. To illustrate the convergence of the objective function of our new method, Figure 3(a) and Figure 3(b) show a typical run of our algorithm on two UCI benchmark data sets. As can be seen from the figures, when the algorithm iterates and the objective value decreases, the accuracy shows a relatively smoothly increasing line.

5.4 Experiments on Noisy Real Data Sets

To study the impacts of the value of p in our new embedding model, we randomly contaminate 20% of the data points in all 7 data sets and we run our method with increasing p on those data sets. For each p, we run 200 times for the same contaminated data and original data respectively. Other algorithms are also tested for 200 times on each data set for comparison. The performances of the clustering methods on contaminated noisy data sets are reported in the bottom half of Table 2. Among all the best and the average values of clustering accuracy and clustering purity, our method is consistently better than its counterparts. The results of our approach generally decreases less than other methods on the contaminated data sets, especially for those noisier data sets.

6 Conclusion

In this paper, we proposed a new robust Laplacian embedding approach that uses the p-th order of the ℓ2-norm distances in the objective and strictly satisfies both orthogonality and nonnegativity constraints at the same time. This results in an objective that is neither convex nor smooth, which is difficult to efficiently solve in general. We thereby proposed a novel smoothed iterative reweighted method to solve this challenging optimization problem, in which a smoothness term is formally and explicitly introduced for improved numerical stability. As an important theoretical contribution of this paper, we rigorously we proved the convergence of our new algorithm with the smoothness term. Using this new and improved optimization framework, our objective can be elegantly solved. We have performed extensive experiments, in which the superior performance of our new method has demonstrated its effectiveness and the potential to give a new perspective for nonlinear graph based clustering tasks.

Acknowledgments

All correspondence should be addressed to: Hua Wang (huawangcs@gmail.com). This work was partially supported by National Science Foundation under Grants IIS-1652943 and IIS-1849359.