Using Graph Algebra to Optimize Neighborhood for Isometric Mapping

1 Introduction

Classic dimensionality reduction approaches can not be reliably applied to find the meaningful low-dimensional structures hidden in the high-dimensional observations for exploratory data analysis such as classification and visualization, so that two new approaches have recently been developed. One is locally linear embedding (LLE) that is based on local approximation of the geometry of the manifold[Roweis and Saul, 2000]. LLE has many variants, such as Laplacian eigenmaps [Belkin and Niyogi, 2003]and Hessian eigenmaps [Donoho and Grimes, 2003], incremental LLE [Kouropteva and Pietikainen, 2005], supervised LLE [de Ridder et al., 2003; Kouropteva and Pietikainen, 2005], integrated LLE with classic PCA or LDA [Abusham et al., 2005; Chang et al., 2004], integrated LLE with SOM [Xiao et al., 2005] etc. The other is Isomap that preserves the manifold geometry at all scales and has better ability to deal with nonlinear subspaces [Tenenbaum et al., 2000]. It also has many variants, such as Landmark Isomap [Silva and Tenenbaum, 2003], supervised Isomap [Geng et al., 2005], spatio-temporal extension of Isomap [Jenkins and Mataric, 2004], incremental extension of Isomap [Law and Jain, 2006], integrated Isomap with fuzzy LDA [Weng et al., 2005] etc. The other related studies are also performed such as the selection of the optimal parameter value for LLE and Isomap [Kouropteva et al., 2002; Shao and Huang, 2005], integration of LLE and Isomap [Saxena et al., 2004] etc. LLE and Isomap have their own superiorities so that they have been being developed simultaneously for various context applications.

Despite Isomap performs well in many cases, it often fails to nicely deal with noisy data or sparsely sampled data [Geng et al., 2005]. This is because in these cases the local neighborhood structure on which Isomap largely depends is critically distorted. This makes Isomap vulnerable to short-circuit errors that some neighors of the current point come from other different folds so that these neighbors are not nearest ones on manifold [Silva and Tenenbaum, 2003]. This can in turn lead to drastically incorrect low-dimensional embedding. Accordingly some supervised Isomap are proposed [Vlachos et al., 2002; Geng et al., 2005] which employ the class labels of the input data to guide the manifold learning. They improve the Isomap in classification and visualization, but they can not work when the class labels of data are not available. Due to that Isomap can not works on the disconnected neighborhood graph, some approaches to constructing connected neighborhood graph are proposed such as the k-connected or k-edge-connected spanning subgraph of the complete Euclidean graph of all data points[Yang, 2004; 2005]. This paper deals with the neighborhood graph from another perspective. It considers the implicit correlation among data points using the path algebra on the neighborhood graph, which is further applied to improve Isomap.

2 Approach to optimize the neighborhood

In many data analysis and pattern recognition tasks, similarity between two objects can be established by direct comparison and induced by mediating objects. Namely, two objects might be considered similar when they are connected by a chain of intermediate objects where all dissimilarities or distances between neighboring objects in the chain are small [Fischer and Buhmann, 2003b]. This concept can be generalized by assuming that object similarity behaves transitive in many applications. Based on this intuition, the path algebra-based clustering approach is proposed that can extract elongated structures from the data in a robust way which is particularly useful in perceptual organization [Fischer and Buhmann, 2003a; 2003b]. This paper applies it to build the neighborhood graph for Isomap. The intuitive picture that two objects should be considered as neighbors if they can be connected by a mediating path of intermediate objects.

Let G = (V, V × V ) be the complete Euclidean graph of all data points. It weights the edges using Euclidean distance de . A path l from v0 to vn in G is defined as a sequence of vertex (v0, · · · , vi, · · · , vn) such that (vi, vi+1) ∈ V × V , 0 ≤ i < n.

Definition Let P(v0, vn)={l|l be a path from v0 to vn in G}, we define path algebra as a set P with two binary operations ∨ and · which have the following properties [B.Carre, 1979]:

1. For all x, y, z ∈ P, the ∨ is idempotent, commutative, and associative: x ∨ x = x, x ∨ y = y ∨ x,(x ∨ y) ∨ z = x ∨ (y ∨ z)
2. For all x, y, z ∈ P, the · is associative, and distributive over ∨: (x · y) · z = x · (y · z), x · (y ∨ z) = (x · y) ∨ (x · z), (y ∨ z) · x = (y · x) ∨ (z · x)
3. The set P contains a zero element φ and a unit element e such that: φ · x = φ, φ∨x = x = x∨φ, e · x = x = x · e

Obviously, there are many ways to define the path algebra based on different intuitions. This paper focuses on a simple but important way to define the path algebra as follows, where for all x, y ∈ R, φ = 0, and e = ∞

In this way a new distance dm can be defined as follows [Fischer and Buhmann, 2003b]:

The difference between dm and de can be illustrated in Figure 1. This difference can significantly impact the construc-

)} , ( )}, , ( ), , ( ), , ( min{max{ ) , ( b a d b d d d c d c a d b a d e e e e m =

Figure 1: The difference between dm and de

tion of neigborhood graph on data manifold. For example, the categories of 5-nearest neighbors of x using de distance,circled as in Figure 2, is

Generally points belonging to the same class are often closer to each other than those belonging to different classes. Obviously the short circuiting problem happens here because points marked by square are not the most similar points to x in terms of categories. By contrast, if dm distance is utilized to determine the neighbors for x, the result is

This is a correct result that indicates all neighbors of x lie on the manifold. It seems that dm can pulls points belonging to the same class closer and propels those belonging to different classes further away. It ensures to preserve the intrinsic structure of the data set and therefore can be applied to improve the Isomap.

3 Proposed Isomap with Optimized Neighborhood

Isomap can perform dimensionality reduction well when the input data are well sampled and have little noise. In the presence of the noise or when the data are sparsely sampled, short-circuit edges pose a threat to Isomap by further deteriorating the local neighborhood structure of the data. Subsequently Isomap generates drastically incorrect low-dimensional embedding. To relieve this negative effect of the noise, dm can be applied to determine the true neighborhoods as opposed to using Euclidean distance de. This can get better estimation of geodesic distances and in turn give lower residual variance and robustness. It also reduces the time complexity because the better neighborhood structure can speed up the subsequent optimization process.

To avoid computational expense to calculate the dm between any two points, we find an approximate way to build the neighborhood graph using dm that can be computed efficiently. It first applies de to build the neighborhood graph, and then utilize the idea of dm to optimize this neighborhood graph. The algorithm is given as follows.

Algorithm 1: OptimizeNeighborhood(X, k, m, d)

Input: X = {xi} be the high dimensional data set, k be the neighborhood size, m be the scope for optimization of neighborhood, and d be the dimension of the projected space.

Output: The optimized neighborhood set N = {N(xi)} for all points.

1. Calculate the neighborhood N(xi) for any point xi using de, where N(xi) is sorted ascendingly.
2. Compute n, which is the number of points in X.
3. For i=1 to n
4. For j= 1 to k
5. Select j-th point fromN(xi), denoted as xij
6. For p= 1 to m
7. Select p-th point fromN(xij), denoted as xijp

1. If de(xij , xijp) < de(xi, xik) and xijp /∈ N(xi) and parent(xijp) ∈ N(xi)
2. Return N = {N(xi)} that is the optimized neighborhoods for X

9. Delete xik fromN(xi) 10. Insert xijp to N(xi) ascendingly 11. Let j=1 12. Break 13. EndIf 14. EndFor 15. EndFor 16. EndFor

In algorithm 1, step 1 calculates the neighborhood graph using Euclidean distance, as the same as most dimensionality reduction approaches such as Isomap do. Its time complexity < O(n3 ). From step 3 to step 16, the neighborhood of the given point xi determined in step 1 is optimized using dm. It includes two cycles. Although step 11 reset j, k and m are neighborhood sizes so that they can be regarded as small constants. It means that optimization of neighborhood of a point can be finished in O(1). Therefore the neighborhood graph can be efficiently optimized with additional time complexity O(n). In Step 8, condition parent(xijp) ∈ N(xi) means that xijp must be reachable from xi in the neighborhood graph determined in step 1.

In order to illustrate the idea behind the algorithm 1, we give an example to optimize the neighborhood of point x, shown as Figure 3, where k=3 and m=2. Figure 3(1) shows the neighborhood graph of x determined using Euclidean distance, where only parts are described. It can be observed that x has neighbors N(x) = {x1, x2, x3}. In algorithm 1, step 5 chooses its first neighbor x1 and then step 7 selects the neighbor x11 for x1. In step 8, x11 will be compared with the largest neighbor of x1, namely x3. Because the condition is satisfied, All steps between step 9 and step 12 will be executed. This results in that x11 will replace x3. Consequently the neighorhood of x is optimized as N(x) = {x1, x11, x2}. It is shown as Figure 3(2). Because the neighborhood of point x has changed, the optimization restarts from the smallest neighbor x1 of the point x and enter the next cycle. Due to x11 has been in neighborhood of x, x12 will be explored, shown as Figure 3(3). Because the condition in step 8 is satisfied, all steps between step 9 and step 12 will be executed. This results in that x12 will replace x2.Consequently the neighorhood of x is optimized as N(x) = {x1, x11, x12},shown as Figure 3(4). Due to m=2, x13 will never be exploited. So far all neighbors of x1 have been tested. Next cycle will explore the second neighbor of x in step 5, namely x11. It goes through all steps as exploring the first neighbor x1 of x, which is not explained more here.

Now OptimizeNeighborhood can be applied to improve the basic Isomap. For hereafter comparison, we denote this improved Isomap as m2-Isomap, which is described as algorithm 2.

Algorithm 2: m2-Isomap(X, k, m, d)

• Input: X = {xi} be the high dimensional data set, k be the neighborhood size, m be the scope for optimization of neighborhood, and d be the dimension of the projected space.

Output: The dimensionally reduced dataset Y = {yi}.

1. Utilize OptimizeNeighborhood(X, k,m, d) to calculate the optimized neighborhoods for each point in input data set X, which will be applied to construct the neighbor graph in next step.
2. Construct the weighted neighbor graph Ge = (V,E) for the dataset X, by connecting each point to all its knearest neighbors, where (vi, vi+1) ∈ E,if xj is a member of k neighbors of xi determined in step 1. This edge has weight de(xi, xj ).
3. Employ the Ge to approximately estimate the geodesic distance dg(xi, xj) between any pair of points as the shortest path through graph that connects them.
4. Construct d-dimensional embedding from Ge using MDS

Compared with the basic Isomap, m2-Isomap adds step 1 and then modifies the step 2. It should be noted that the edges of neighborhood graph are weighted by de instead of dm. The dm is only applied to optimize the neighborhood rather than applied to estimate the geodesic distance. The step 3 and step 4 remains the same as the basic Isomap. Therefore the time complexity increases from step 1 where the additional complexity is O(n). However the optimized neighborhood will speed up the step 4 so as to decrease the running time of the whole algorithm. This will be proved empiricially in later experiment.

4 Experimental Results

Several experiments are conducted to compare Isomap with m2-Isomap in visualization and time complexity. Isomap uses the published matlab code [Tenenbaum et al., 2000]. The m2-Isomap is also implemented in MATLAB 7.0. In experiments, there are two parameters involved. Isomap has a parameter k to determine the size of local neighborhoodwhereas m2-Isomap introduces an additional parameter m that controls the scope of the neighborhood to be modified, where k=5-15 and m =1-10 will be explored.

A. Swiss Roll Data set

Swiss Roll dataset is widely applied to compare many non-linear dimensionality reduction approaches [Roweis and Saul, 2000; Tenenbaum et al., 2000; Balasubramanian and Schwartz, 2002; Geng et al., 2005]. We take many random samples from the Swiss Roll surface shown as Figure 4, and do visualization experiments on them to compare Isomap with m2-Isomap in the following cases: (1)well sampled data without noise (2) well sampled data with noise (3)sparsely sampled data.

The quality of visualization can be evaluated subjectively, while it is often quantitatively evaluated by the residual variance [Roweis and Saul, 2000; Tenenbaum et al., 2000; Balasubramanian and Schwartz, 2002; Geng et al., 2005]. The lower the residual variance is, the better high-dimensional data are represented in the embedded space.

Experiment 1 On well sampled noiseless data sets

To make comparison on well sampled noiseless data sets, five sample data sets with 1000 points are sampled randomly from Swiss Roll surface. Firstly on the first sampled data set, two approaches are tested to choose the optimal parameters in the given scope, and then these parameters are applied to run the remaining four data sets. The mapped results from left to right corresponding to five data sets are shown as Figure 5, where the parameters and residual variance for each data set are added as caption. For example on the first data set,in Figure 5(a)Isomap, the caption "6/0.0015475" indicates that the parameter k takes 6 whereas the residual variance is 0.0015475. In Figure 5(b)m2-Isomap, the caption

"12/6/0.00059146" indicates that the parameter k takes 12,m takes 6, and the residual variance is 0.0.00059146. It can be observed that on well sampled data sets without noise, m2Isomap outperforms Isomap. It gets better mapping results and lower residual variances on any data set.

Experiment 2 On well sampled but noise contaminated data sets

To compare the topologically stability of Isomap and m2Isomap on data sets with noise, we do experiments on the first data set with 1000 points used in experiment 1, but here it is contaminated by adding random Gaussian noise. The mean of the noise is 0 and the variance is 0.64. It can be seen from Figure 6 that Isomap gets the best result with residual variance 0.001286 when k=6, while m2-Isomap gets the best result with the lower residual variance 0.0010773 with parameters k=12 and m=6. Obviously m2-Isomap outperforms Isomap in terms of visualization and residual variance. It can be also observed from Figure 6 that as k increases, the results of Isomap change drastically, but those of m2-Isomap do not change much. This indicate that m2-Isomap is also less sensitive to k than Isomap.

Experiment3 On sparsely sampled data sets

Generally, in sparsely sampled data sets, the Euclidean distance between points in neighborhood becomes larger as compared to the distance between different folds of the manifold. This easily makes two approaches faces the problem of short-circuiting. This experiment makes comparison between two approaches about robust to samples with different sampling density. Five data sets with 300,400,500,600 and 700 points respectively are randomly sampled from Swiss Roll surface. Two approaches are tuned to select the best visualization results on each data set. It can be observed from Figure 7 that m2-Isomap outperforms Isomap on each datset in terms of visualization performance and residual variance.

B. Face image data set

This data set consists of 698 images (each contains 64×64 pixels) of a human face rendered with different poses and lighting directions [Tenenbaum et al., 2000]. To compare two approaches on sparsity data set, we choose 350 (half of the data set)images from this data set to form a new data set. The mapped results are shown in Figure 8 and Figure 9. In fig-

100

80

60

40

20

0

-20

-20

0

20

40

60

80

100

120

140

100

80

60

40

20

0

-20

-20

0

20

40

60

80

100

120

140

160

Figure 9: m2-Isomap (k=5 m=1)

ures, the x-axis represents the left-right poses of the faces, and the y-axis represents the up-down poses of the faces. The corresponding points of the successive images from left to right in the middle are marked by circles and linked by lines. The nine critical face samples are marked by plus at the leftbottom corner of each image indicates the point representing the image. It can be observed from these figures that Isomap can hardly reveal the different face poses. The middle leftright line is heavily curved, and the arrangement of the nine face samples is tangle some. By contrast, m2-Isomap puts the middle left-right line better. The nine face samples are also mapped to the approximately right positions corresponding to the face poses. These indicate that m2-Isomap performs better than Isomap.

This also can be supported from their residual variances. It can be observed from Figure 10 that m2-Isomap gets lower residual variances and the correct intrinsic dimensionality of data set, while Isomap gives the higher estimate to intrinsic dimensionality than that of the data set. The intrinsic dimensionality of the data can be estimated by looking for the "elbow" at which this curve ceases to decrease significantly with added dimensions.

C. Time complexity comparison

We generate 10 data sets with the same size randomly sampled from Swiss Roll surface. Two approaches run on these ten data sets respectively and then average running time for two approaches are calculated. And then the experiments are performed on different sizes: 500 900 1300 1700 2100 2500 2900 3300. It can be observed from Table 1 that m2-Isomap also exceeds the Isomap in time complexity. The larger the data set size is, the superiorities of m2-Isomap is more obvious. This is because good neighborhood structure of data set is beneficial to speed up the subsequent optimization process.

5 Conclusion

This paper applies the graph algebra to build the neighborhood graph for Isomap. The improved Isomap outperforms the classic Isomap in visualization and time complexity. It also has stronger topological stability and less sensitive to parameters. This indicates that the more complicated or even time-consuming approaches can be applied to construct the better neighborhood structure whilst the whole time complexity will not raise. Because the transitivity of similarity is not absolutely correct for any data set, in the future, we will apply the fuzzy theory and probability theory to define the graph algebra and then to build the neighborhood graph. And besides, the proposed approach should be combined with kedge-connected spanning subgraph so as to guarantee to build the connected neighborhood graph for any data set.

Acknowledgments

This work is supported by Natural Science Foundation of China(60003019),Hubei Science and Technology Project (2005AA101C17),and UK Engineering and Physical Sciences Research Council under grant number GR/N15764/01.