Visualizing High Dimensional Data

Quantifying dependency between attributes.

2011-11-30T13:29:00.000-08:00

Given two sets of attributes, X = (x₁, x₂, ..., x_p) and Y = (y₁, y₂, ..., y_q), how do we quantify the dependency between X and Y?

If both X and Y are one dimensional numerical attributes, we can use linear correlation coefficient to quantify the dependency between X and Y. For multidimensional X and Y (that is our data are numerical matrices of dimension n×p and n×q respectively, assuming that we have sampled n data points) we can use canonical correlation analysis (CCA) to capture the dependency between X and Y. Linear methods are in general effective and simple to use. VisuMap offers support for CCA and many other linear analysis services through a plugin module MultivariateAnalysis. However, linear methods work only for numerical attributes and can only capture linear relationships.

In this note I'll describe a more general method to quantify dependency between attributes. The method can be applied to non-numerical data and can capture non-linear relationships. The method is based on the concept conditional entropy from information theory. Assuming that X and Y are two random variables, the information gain ratio from X to Y is defined as following:

where H(Y) is the entropy of Y, and H(Y|X) is the conditional entropy of Y under the conditional that X is known. The numerator H(Y) - H(Y|X) is also called the mutual information between X and Y, and it is normally denoted as I(X; Y). So, the above formula can also be defined as R(X;Y):=I(X;Y)/H(Y). The following diagram shows the relationship between various types of entropy relating to two random variables:

In general, R(X; Y) tells us how much information we gain about Y when the value of X becomes known. R(X; Y) is a number between 0 and 1.0. If R(X;Y)=0, then X and Y are independent from each other, knowing X does not provide any information about Y. If R(X;Y)=1.0, then Y is completely dependent on X, i.e. Y is completely determined by the value of X.

Now, let's go back to our original problem. We assume that X and Y are now two tables of dimension n×p and n×q respectively, and X and Y may contain numerical and non-numerical values. In order to calculate the information gain ratio between X and Y we first convert them to random variables. To do so, we pick two clustering algorithms C and D that can be applied on X and Y. The clustering algorithm will assign each of the n data point to a label of label sets X and Y respectively:

C: N → X

D: N → Y

where N:= {1,2,3, ..., n} is the set of data sample indexes. We then define two random variables X' and Y' over X and Y as following:

With above definitions we can then apply the formula for R(X;Y) to the random variables X' and Y'. In particular, we have the following detailed formula to calculate the generalized information gain ratio (GIGR) of X for Y:

We notice that GIGR in above formula is indexed by the two clustering algorithms C and D. This means that the calculation is dependent on what clustering algorithms we choose for X and Y. In other words, the method described here is actually a meta-algorithm that requires, apart from two input data tables, also two clustering algorithms as input.

The clustering algorithms C and D plays an essential role in GIGR. They determine the granularity of the entropy calculation. More particularly, R_C,D will be maximum if C and D are the finest clustering algorithms by creating unique label for each distinguished data point. On the other side, R_C,D will be zero, if C and D, as the coarsest clustering algorithm, simply assign all data points to single label.

We can also consider the clustering algorithm in GIGR as a quantization method that converts data from multidimensional discrete or continues domain to one-dimensional discrete domain, so that we can effectively apply the information theoretical calculation on them.

In practise, C and D can be different clustering algorithms. The main requirement on C and D is that they can be applied on data tables X and Y separately and can assign a label to each data row. In fact, the data X and Y don't even have to be tables. They can be, for instance, a tree structures or DNA sequences. To apply our method, we just need to pick appropriate clustering algorithms for the data types.

It should be pointed out here that GIGR described here does not bias towards any type of clustering algorithms. In fact, I believe that any clustering algorithm would be appropriate for certain type of problems. The attached script code implemented 6 different groups of clustering algorithms to calculate GIGR with VisuMap.

Since the two clustering algorithms in GIGR are kind of input, we can use GIGR to investigate similarity between different clustering algorithms. For instance, if R_C,D(X; X) is close 1.0 for a data set X, then the two clustering algorithms C and D produce similar results with respect to our data. Thus, in the practise, we might use the algorithm that is easier to implement, if the GIGR indicates that two algorithm produce similar results.

Also along the same line, GIGR can be used as an indicator for the stability of clustering algorithms. Many clustering algorithms, like the k-mean algorithm, are non-deterministic in the sense that repeating the algorithm on the same data normally yields different result. We would normally favor an algorithm that maintain certain level of stability so that the clustering results remain more or less the same. When we calculate R_C,C(X; X) by repeating algorithm C twice on X, we would normally get a GIGR value that is smaller than 1.0 for non-deterministic clustering algorithms. In these cases, there closer R_C,C(X; X) is to the value 1.0, the more stable is C with respect to the data set X.

Appendix:
The following JavaScript code helps users to calculate the GIGR with VisuMap software application:

K-Mean clustering under no-optimal conditions

2011-11-29T12:40:00.001-08:00

K-Mean algorithm (KMA) is a widely used clustering algorithm for multidimensional data. In this note, I'll show some intuitive behaviors of KMA with a very simple quantization task: Grouping a series of numeric values into multiple groups. An intuitive requirement for this task is: if some values are concentrated at a location, those values should be grouped together.

I have generated 5 datasets, each with 5000 values between -5.0 and +20.0. Data points in the 5 datasets are successively concentrated at 1 to 5 locations according to Gaussian distribution. I used KMA to cluster these data into 3 clusters. I ran KMA on each dataset twice; the result is shown in the following picture:

In above picture, each spectrum shows the result of one run of KMA on a dataset. Each bar in the spectrum indicates a data point in a dataset. The curve on the spectrum shows the density of the data points. The bar colors indicate the clusters created by the KMA. We notice that all bars are colored by 3 different colors as the number of clusters, K, used KMA is 3.

Let N be the natural clusters present in the dataset (i.e. N=1,2,3,4,5 for these 5 datasets,) it would be helpful to distinguish the behaviors of KMA in to the following three cases:

(1) N < K : In this case the KMA has to group fewer natural clusters into larger number of clusters. As be shown in spectrum 1.a, 1.b, 2.a and 2.b, KMA assigned more than one clusters (colors) to a natural cluster, but didn't assign one to multiple natural clusters. For instance, in spectrum 2.a, the first natural cluster has been split evenly into two clusters.

(2) N = K : In this case the KMA has to group exact N natural clusters in to N clusters. We see in spectrum 3.a and 3.b that KMA did a nice job in this case. KMA acutally found the natural clusters.

(3) N > K : KMA has to group more natural clusters into fewer clusters. As be shown in spectrum 4.a, 4.b, 5.a and 5.b, KMA, in this case, grouped multiple natural clusters together. KMA has not, at least not significantly, split natural clusters. How KMA groups the natural clusters can change from run to run, depending on the initialization of the algorithm. For instance, spectrum 4.a and 4.b show two different grouping created by two different KMA runs.

In above test, our datasets all have clearly separated natural clusters. I would say that KMA did a pretty good job: its results reflect pretty much our intuitive expectation. However, if the natural cluster are not so clearly separated, the result of KMA would not be such nice. The following picture shows how KMA groups less separated natural clusters:

In above picture, the 5 datasets are created similarly as before except that these natural clusters overlap each other somewhat. We see that, only for the case K=N, KMA found the exact natural clusters. For other causes, KMA does not respect the natural clusters. Actually, we can visually do a better job based on the density curves.

VisuMap v3.5 Relased.

2011-05-13T13:25:00.000-07:00

Today we have released VisuMap version 3.5.864. This version is a major upgrade from its previous version 3.2; it is marked by the following major changes:

Multi-core optimization: we have implemented extensive optimization for most algorithms to exploit multi-core computer system. By default, VisuMap 3.5 will automatically make use of all available CPU cores with multi-threading technology. The user can also explicitly limit the number of processing unit through the application configuration property WorkThreads. On our 6-core Windows 7 test machine we have often achieved 3-5 times performance speedup.
Updated DirectX Library: VisuMap used DirectX+MDX (managed directX) to implement hardware accelerated visualization service. Since MDX has been discontinued as a component in new generations of .NET platform, we have to switch to another product. In VisuMap v3.5 we have replaced MDX by SlimDX that is an open source library with similar goal as MDX. With this change, we not only significantly improved the rendering performance, but also reduced the dependency on old obsoleted software components. VisuMap can now be installed on most Windows 7 machines without pre-installing any other software. Some older systems may need to upgrade to the latest DirectX runtime library through the normal DirectX web installer.

Many bug fixes and enhancements: Based on customer feedback VisuMap has undergone many small changes. Most of them have already made available through previous minor releases. Worth mentioning is here the extension of the XY plot that made this simple service a generic tool for large datasets; as well as addition of new glyphsets and utility scripts. We have recently also released some VisuMap plugins which extend VisuMap for special user groups. For instance, the plugin for Multivariate Analysis and Wave Transformations. Also, the Customer Importer plugin has been extended to import FCS (Flow Cytometry Standard) files.

VisuMap v3.5 is backwards compatible with its previous versions. All existing datasets will silently upgraded to the new version when loaded into VisuMap. Existing users can upgrade to the new version, as long as the license type allows it, by simply clicking on the menu Help>Check for Updates in VisuMap's main window.

Self-similarity of high dimensional random walk process

2011-04-06T14:55:00.000-07:00

A high dimensional random walk process (RWP) is the trajectory of a vector variable whose components change independently and randomly step by step, in discrete time space, for a small constant percentage. High dimensional RWP can be used as the starting point to investigate changing complex systems. For instance, as discussed in this blog, the stocks price history of 500 stocks can be considered as a 500 dimensional RWP.

Self-similarity means that an object is, from certain point of view, similar to parts of itself. A high dimensional RWP is self-similar in the sense that each sub-section of it follow the same statistical constraints. Thus, the randomness is a property shared by RWP and its sub-sections. We can generate series of 1000 data points of a 500-dimensional RWP with the following VisuMap's JavaScript code:


var n = 1000;
var dim = 500;
var rwp = New.NumberTable(n, dim);
var m = rwp.Matrix;
for( var col=0; col < dim; col++) {
    m[0][col] = 1.0;
    for(var t=1; t < n; t++) {
        m[t][col] = m[t-1][col] * (1 + 0.01*(Math.random() - 0.5) );
    }
}
rwp.ShowValueDiagram();

How can we visualize the self-similar randomness in high dimensional RWP? A simple way is to use a dimensionality reduction method to map the high dimensional trajectory to low dimensional space, then plot it out on paper or screen. The following picture for instance shows two CCA (curvilinear component analysis) maps of the 1000 data points mapped to the 3 dimensional space:

We can see the randomness of the trajectory in above picture, but the self-similarity is hardly apparent. This is because, as intrinsically to the human perception, we normally only good at recognizing similarity between geometric patterns, but not similarity between random patterns.

Fortunately, when we use principal component analysis (PCA) to project the trajectory to low dimensional space, we get much easier recognizable patterns. For instance, the following picture shows the projections of our sample RWP to some major principal components (ie. eigenvectors):

The PCA projection to the major principal components apparently filtered out the randomness of the data, like a low pass filter suppresses high frequency noises. Now, we can select three principal components, say the second, third and forth components as our projection axes; then plot the PCA projection of the data (or sub-sections of it) to these three axes. As be illustrated in the following picture, we can see that they are all geometrically similar to each other albeit with different densities.

Notice that in above picture, the projection axes are the second, third and forth principal components of the selected data, not those of the complete dataset. The following video clip shows what we have discussed above in a more intuitively way:

One practical use of studying RWP is to find non-randomness ( that is information ) in seemly random data. Using the PCA technique we can show random data as easy recognizable patterns which enable us to detect deviations from those patterns, so that we can quickly find potential information in apparently random data.

Pirated software web sites

2011-04-01T08:21:00.000-07:00

I have notice that more than a dozens so-called freeware web sites now provide pirated version of VisuMap software. All those sites seem to be originated from a single source, as they are all organized in a similar way and all pirated the version 3.2.854 of VisuMap. I have spent 2 or 3 hours to investigate those sites, the followings are what I have found:

These sites have done a pretty good SEO (search engine optimization) job, they pretty much over flooded search results when I searched for "VisuMap download" on Google. Google is clearly loosing the battle against those mis-information.

Most of these sites are not really free. They usually re-direct users through a chain of redirections with dubious ads (they hold the user with popups for popups, you need to shutdown the browser to get away from them), then at the end they usually require users to sign up for paid service for "fast" download.
All product information about VisuMap on those sites are copy and pasted from VisuMap's web site. Some sites used translation engine to translate English to other languages with horrible results.
All information about downloading sites, speed, etc. on those site are faked and likely automatically generated to fool users and search engines.

I have tested the pirated software on an isolated machine. I scanned the software with anti-virus software, no virus or spyware has been detected.

The pirated software seemed to run smoothly on the test machine, although it felt a little sluggish. Since it has no proper license, it won't be able to get online update from VisuMap's web site.

The pirated software has an obsoleted version, so that it won't be able to use some new plugin extensions for VisuMap.

With so many limitations and apparent dubious practice, I don't think those pirate sites will have serious impact on the legitimated software vendors. But, they clearly post a challenge for Internet searching engines.

Visualizing dynamics in time series data.

2010-11-04T15:52:00.000-07:00

We have just released a short video clip that shows how VisuMap visualizes dynamics in time series data. This demo uses a feature of the spectrum view that smoothly changes the display from one configuration to another one. This video is inspired by the some samples from google data explorer. It is amazing how human eyes can quickly capture such characteristics like speed and variations if the data is properly visualized.

This video uses some features in VisuMap available in VisuMap version 3.2.855. The sample datasets used in the video are included in the standard installation of VisuMap.

Principal Components Analysis (PCA) of Random Walk Process

2010-10-04T09:26:00.000-07:00

One sample dataset distributed with the VisuMap software package is the normalized weekly stock price history of 500 S&P stocks for the year 2002. This dataset as shown in the following picture can be considered as 500 time series for about 50 time points (i.e. weeks).

Picture 1

While above picture shows some common trend (like the large downturn in middle of the year), the price development is predominately random. How can we characterize this randomness? Well, one way to model the randomness is considering the stock price developments as independent random walk processes: Assuming that you have invested one dollar in each of the 500 stocks, their values will change more or less randomly from week to week for a small percentages. As a reference model we can generate 500 such random walk process with the following JavaScript code (as supported by VisuMap):

var randomWalk = New.NumberTable(500,50);
for(var row=0; row < randomWalk.Rows; row++) {
  var v = 1.0;
  for( var col=0; col < randomWalk.Columns; col++) {
    randomWalk.Matrix[row][col] = v;
    v *= 1 + 0.01*(Math.random() - 0.5);
  }
}
randomWalk.ShowValueDiagram();

The last line in above code block creates the following diagram (the diagram has been colored with k-mean algorithm):

Picture 2

Comparing Figure 1 and 2, we can say that Figure 2, more or less, resembles Figure 1 if we remove some common trends from Figure 1 (i.e. de-trending the down-turn in the middle and the overall slightly down-wards trend). The resemblance is however a collective similarity. It
does not make sense to find similarity between individual curves in the two figures.

Another way to compare two sets of time series is using principal components analysis (PCA). Among many other services, PCA provides a systematic way to decompose high dimensional variances in a dataset into many 1-dimensional directions. These directions (called principal components) are also ordered in the way, so that the beginning components have larger variance and therefore more information. Thus, PCA might be interesting way to characterize the variance (i.e. some kind of randomness) in a dataset.

It is very simple to visualize the PCA components with VisuMap. To do so, we just need to open the PCA view and open the Projection Analyzer window. Then select some PCs and choose the context menu "Show PCAs". The following two picture shows the first 24 PCs of the S&P500 and the RandomWalk datasets:

Picture 3

Picture 4

We can clearly see some similarity between the first few PCs in above two pictures. This indicates that random walk process models pretty well the major variance directions of the S&P dataset. The discrepancy between the PCs in above two picture may be used to characterize how the S&P price history differs from random walk process. In this way the random walk process is used as a null-space reference model.

Looking at Picture 4 we can also notice an interesting thing: the PCs clearly resembles the curves of sine functions. This is actually not too surprising since the PCA algorithm is basically a sequence of high dimensional rotation operations, which, as we know, lead to a lot of sine/cosine functions. Nevertheless, it would be interesting to determine the exact mathematical formula for the random walk process. By doing so, we can have a quick statical approximation for PCAs of many random walk alike datasets.

Visualizing linkage disequilibrium clusters of genotype SNPs

2010-07-18T11:08:00.000-07:00

SNP and Linkage Disequilibrium

I have just released a sample dataset Haplotype Analysis that shows how to use VisuMap to visualize clusters among SNPs (single nucleotide polymorphism). SNPs are about 0.1% of base-pairs locations in DNA sequences that vary from population to population; and from individual to individual. It has been said that the difference of many phonetic traits, like height, eye color, etc, of human beings can be attributed to variations of SNPs.

Haplotype analysis aims to find correlation between SNPs and phonetic traits. One of the often used method in haplotype analysis is the concept linkage disequilibrium (LD) that measures the correlation between SNP pairs with respect to a population. In abstract terms, LD induces an information structure over SNPs. Such a relational structure may offer a helpful framework to study the correlation between SNPs and phonetic variations.

The numerical calculation of LD is actually pretty straightforward. But all descriptions about LD I can find are too vague for step-by-step calculation. For the sack of reference, I will briefly describe the way how LD is calculated here. Assuming that we have obtained the genotype data about two SNPs, SNP_a and SNP_b, for a group of individuals I₁, I₂, I₃, ..., I_k as in the following table:

	I₁	I₂	I₃	...	I_k
SNP_a	T/T	G/T	T/T	...	G/G
SNP_b	C/T	C/C	C/T	...	T/T

Notice that each element in above table is a pair of nucleotides A, C, G or T. For certain reason (which I don't know), all SNPs are bi-allelic, that means each row in above table will only have maximal two different nucleotides. So, in above sample, SNP_a has G and T, whereas SNP_b has C and T. In order to calculate LD between the two SNPs, we first select a nucleotide arbitrarily from each row; and then count how many times that nucleotide appears in that row in each column. For instance, if we have selected T and C for SNP_a and SNP_b respectively, we will get a frequency table as following:

	I₁	I₂	I₃	...	I_k
SNP_a	2	1	2	...	0
SNP_b	1	2	1	...	0

The LD between the SNP_a and SNP_b (more accurately the R-Squared LD) is then defined as the squared value of the linear correlation coefficient between the two frequency rows. Notice that since the SNPs are bi-allelic, the final LD value won't depend on which nucleotide we have selected from each row for the calculation. We can verify this easily based on the formula below. More formerly, let {a₁, a₂, a₃, ..., a_k} and {b₁, b₂, b₃, ..., b_k} be the two frequency row vectors in above table; let a and b be their mean values; then the R-Squared LD between SNP_a and SNP_b is defined as following:

The Sample Dataset

The uploaded sample dataset is obtained from the HapMap Project website that provides accesses to a large collection of genotype datasets. More interestingly, from that website we can select genotype data for special SNPs and populations. To generate the sample dataset I have select all SNPs from the region of the first 2 Million base pairs in the 9-th chromosome for the CEU population that consists of 169 individuals. There are altogether 1900 SNPs in the selected region. Our sample dataset is thus basically a table with 1900x169 nucleotide pairs.

In order to import the downloaded data into VisuMap we have implemented a plugin module HaploExplorer that enables users to import data downloaded from the HapMap website (files must have the extension .hmap.) The HaploExplorer plugin also implements a special metric named "LD R-Sqaured", so that users can simply select this metric to generate maps to visualize LD values. The following picture shows a sphere map of the 1900 SNPs:

In above map, each glyph represents a SNP from the dataset. The size of the glyph indicates the chromosome location of the SNP, ie. smaller glyphs represent SNPs located at the beginning of the chromosome. Most importantly, the distances between the glyphs indicate the LD between the SNPs; that means closely located glyphs correspond to SNPs with high LD. We can see clearly various type of clusters in above picture.

We can also visualize the chromosome locations more directly with the help a spectrum view. The following picture shows, for instance, a 2D-RPM map of the the dataset together with a spectrum view of their chromosome locations:

When interactively exploring about views, the user can select a group of SNPs in the lower window, the upper spectrum window will automatically show the chromosome locations of the selected SNPs.

For comparison purpose the HapMap project recommends a program called Haploview that allows user to directly visualize LD values with a triangle matrix, where the LD values are represented by different levels of gray colors. The following picture shows such a map generated by Haploview:

In above picture, all 1900 SNPs are sequentially lined up to the upper edge, so that the matrix is basically a 1900x1900 triangle matrix. In order support exploration of such large matrix Haploview implements an overview window (displayed on the lower left corner) that allows users to select a small section of the data for close investigation.

Comparing above maps, we can see that VisuMap provides a more direct and intuitive way to visualize patterns among SNPs. More importantly, as an integrated software package, VisuMap offers simple framework to investigate different type of LD relationships between SNPs. We can, for instance, easily experiment with any other comparable distance metric available in VisuMap; and we can use any of clustering algorithms in VisuMap to cluster the SNPs.

At last, but not at least, after we have imported the data in to VisuMap, we can also visualize patterns among populations by transposing the dataset table. In order to study relationships between individuals or populations a special distance metric, called the IBS (identity-by-state) distance, has been suggested by some researchers. The IBS distance metric is also implemented in the HaploExplorer plugin. After we have transposed a SNP dataset, we can select the IBS distance metric to produce population maps.

Spherical Multidimensional Scaling in RPM way

2010-03-18T09:55:00.000-07:00

We have just released VisuMap version 3.0.844. In addition to many enhancements (like the map gallery view), this release includes the implementation of a new kind of MDS algorithm, called the manifold RPM algorithm. Manifold RPM works similarly as the original RPM (relational perspective map), except that it maps data to different 2-dimensional surfaces.

The manifold RPM service in this release supports following surfaces as image space:

Flat real projective plane.
Flat klein-bottle
Flat sphere;
3D-sphere.

Whereas the first 3 image spaces so far didn't seem to produce significantly new results, the

3D-sphere surface (also termed as S2) has produced surprisingly good results.

The spherical RPM often produced better results than the original toroidal RPM in the sense

that spherical RPM maps often have less ad-hoc fragmentations. The main reason for this improvement is probably because that 3D-sphere is more symmetrical than torus. Flat torus is symmetrical in shifting (isometric), but not symmetrical in rotation (isotropic), so that some directions are treated differently as I have once blogged. Obviously, 3D-sphere is both isometric and isotropic.

From the point of view of MDS (multidimensional scaling), the spherical RPM algorithm

basically replaces the distance metric by angle metric in image space. Although many trigonometrical calculation are invoked in the algorithm, the implementation of the algorithm has turned out to be significantly faster than the original RPM algorithm.

The major disadvantage of the sphere map is that 3D sphere maps is a little difficult

to explore on 2D computer screens or printed papers. In this release we have implemented a 3D viewer, the sphere view, to help people to explore 3D sphere maps. With the growing popularity of tools like Google Earth, I hope people will find a easy and useful tool in sphere viewer. The following link is a short video demonstrating the 3D sphere viewer implemented in VisuMap:

"Data Earth" - Exploring Data with Sphere Map.

New productivity features in VisuMap 3.0

2009-10-24T11:15:00.000-07:00

I have just posted a short video that shows some new productivity features implemented in VisuMap version 3.0.

A layman's introduction to principal component analysis.

2009-09-23T16:11:00.001-07:00

I have just uploaded a very short introduction for PCA to youtube. The video is probably the shortest introduction to PCA (1.5min). This introduction emphasizes the geometrical aspects, instead of the usual statistical nature. I especially like the statement that the covariance matrix is just a method to measure the average extend of an object (a point cloud) along any axis.

VisuMap 3.0 Released

2009-08-31T10:23:00.000-07:00

Today we have released VisuMap 3.0. This is a major upgrade from the previous version 2.7.

The main change in this release is the implementation of the Atlas service. An atlas in VisuMap 3.0 is frame-less window that enables users to organize and compose data views of different types. The following picture shows the snapshot of a sample atlas:

Atlas service in VisuMap is thus similar as dashboard service in some other visualization centric software applications. It allows users to group visual information according to the subject instead of information type itself. And more importantly for VisuMap users, atlas service enables users to combine multiple algorithms to generate and compose heterogeneously structured information visualization. For instance, this release includes a script (PartitionAtlas.js) that automatically partitions a dataset in several clusters, and generates various maps for each cluster.

Together with this release we have extensive updated the scripting interface. Unfortunately, some old scripts may need minor adjustments to convert to this new release. We have also extended the data format extensively to accommodate the new service. However, the data format has kept compatible in the sense that old dataset files can be loaded into new version, but datasets of new format can not be read by previous versions.

Click the following image to watch a short video that shows how to create an atlas interactively:

Here is another video that shows how a script uses kmean algorithm to partition a dataset and generate different data views for the partitions.

This video shows using pre-configured SOG cluster to create mutiple data views for a dataset.

Data Clustering with Self-Organizing Graph

2009-07-01T09:00:00.000-07:00

We have just released VisuMap version 2.7.832. In this release we have added a very interesting clustering service called the self-organizing graph (SOG). SOG is basically an extension to the self-ogranizing map (SOM) that simulates a homogeneous artificial neural network. Unlike SOM, SOG simulates a network of arbitrary structure. The network can be in fact any weighted undirected graph. The network is a kind of parameter for the algorithm: the user defines the network depending on his/her knowledge or interest about the dataset.

The following picture illustrates an application of SOG. On the left side is a map of a dataset that shows roughly two data point clouds. The right side shows the network we used to classify the dataset. During the learning process each node of the network will become associated with some data points. After the learning process has finished each data point will be classified as their corresponding nodes in the SOG network.

We see that the SOG algorithm not only correctly classified the two data point clouds, but also captured those data points located a the peripheries of the clouds (data points shown as blue triangles).

The following is a short video tutorial for the SOG clustering service in VisuMap:

Comparing with traditional clustering algorithms, like k-Mean and hierarchical clustering algorithms, SOG provides more flexible clustering service. With SOG the user can not only specify the number of clusters, but also the size and other geometric or topological relationships between the clusters. In fact, SOG enables the user to specify a model to reflect existing knowledge and interests about the data.

There are many variants of the self-organizing map that use irregular networks. Those algorithms often employ adaptive strategies to generate the irregular networks depending on the data. Those irregular network are results of those algorithms. Contrary to those algorithms SOG uses the network as input for the algorithm, the user can express his/her interest and knowledge about the data through the network.

SOG is our first attempt to provide a new kind of clustering service. While its framework resembles to provide general pattern matching capability. The learning algorithm it employed has only limited capability to capture structural information. In view of this there are still many challenges and potentials for SOG algorithm in the future.

Visualizing the real projecitve plane.

2009-03-13T15:14:00.001-07:00

I have just uploaded a video clip to YouTube that shows some ways to explorer the real projective plane with VisuMap. The link to this video is included here:

This dataset contains 14,358 data points (4-dimensional vectors) sampled from the surface of RP2 embedded into R⁴. The analytical formula is given in the video clip. A smaller sample dataset in VisuMap format can be downloaded at KleinBottle4D that also contains the dataset for the Klein-Bottle surface.

It is interesting to notice that RPM algorithm flattens the 2D sphere by cutting it into two equal sized half; but it flattens the real projective plane in the way that the whole plane is kept in a single connected region.

Projective and Injective Exploration Methods

2009-03-09T10:26:00.000-07:00

In a recent blog post Terence Tao has put mathematical approaches to study a space X (which can be a space, a manifold etc.) into two fundamental camps:

1. By looking into a map f: Y->X that maps certain (better known) space Y into X.

2. By looking into a map f: X->Y that maps X into certain (better known) space Y.

I think that these statements not only apply to mathematical studies, they also provide useful framework for the exploration of high dimensional data. In our case, the unknown space X becomes a high dimensional dataset that we wish to understand. The space Y is usually a lower dimensional space that we know better and is normally easy to visualize.

For the sake of simplicity I would like to refer to methods from the above two camps as injective and projective methods respectively. That means, injective methods inject known space into the target space, whereas the projective methods project the target space into a better known space.

Following this point of view, algorithms to visualize high dimensional data can be grouped as injective or projective. Most MDS (multidiemsnional scaling) based visualization algorithms are projective. Algorithm like Sammon map and PCA (principal component analysis) map the high dimensional dataset into a low dimensional Euclidean space. The RPM algorithm (relational perspective map) maps high dimensional to flat torus.

On the other side, most clustering algorithms can be considered as injective algorithm. For instance, the method of vector quantization maps a finite discrete set into the target space, whereas the Kohonen net (or Self-organizing map) maps a discrete set with a mesh topology into the target space. The hierarchical clustering algorithms (like the agglomerative clustering algorithm) basically map binary tree structures into the target space.

One major benefit of this classification of algorithms is that it can guide us to apply existing algorithms to new data, in the sense that algorithms from the same camp might be able to share some methods/concepts/tricks. For instance, we might be apply an optimization method used by one clustering algorithm to another injective algorithm.

Furthermore, as be pointed out in Tao's blog post, there is certain "duality" between injective and projective methods. That means, for a given method in one camp you can sometimes find its corresponding counter part in the other camp. This duality can then help us to transfer a method from camp to anther one. As an example for this "duality", the RPM algorithm is basically a dual algorithm of the self-organizing map: they are in different camp, but are dual to each there as they share the flat torus as the image space (i.e. the space Y mentioned at the beginning of this post).

Similar to the duality concept I have blogged before about the "symmetry" between these two camps in more restricted situations, namely when the datasets are multidimensional vector spaces (e.g. data tables). In these cases, we can often transfer a method from one camp to the another one in a straightforward way (e.g. by transposing the data matrix).

The injective and projective algorithm are realized in VisuMap by clustering and mapping services, which form the two main groups of services in the VisuMap application.

Dimension analysis on mulitdimensional scaling

2009-03-05T09:09:00.000-08:00

This post brings forward an inconsistency problem in the study of multidimensional scaling from the perspective of dimensional analysis. This post is rather of philosophical nature. Nevertheless, I find it is worthwhile being aware of this issue
for further theoretical investigations.

Dimensional analysis is a simple and useful rule to check the consistence of a mathematical formula or equation. Roughly speaking, this rule demands that both sides of an equation (that has any physical or chemical bearings) must have the same dimension, or unit of measurement. So for instance the equation E = mC² is consistent with the dimensional analysis, since both sides of the equation have Energy as dimension. The equation E=mC is not consistent, since the right side would have momentum as dimension.

Similarly, if a formula is composed of multiple terms by addition, then all terms must have the same dimension. So for instance, the formula: S:=π·r² + r²/2 makes sense, but not S:=π·r² + r/2.

Now, let's turn to a problem of multidimensional scaling (MDS) from the perspective of the dimensional analysis. Roughly speaking, MDS maps data points from an abstract space with a distance function to a low dimensional, normally, Euclidean space like R². The goal of a MDS algorithm is to construct the map in a way, so that the distance in the low dimensional space approximates the distance in the original abstract space.

More formally, let D_ij be the distance between two data points in the original space, d_ij be the distance low dimensional space after MDS mapping; Then a MDS algorithm tries to minimize the following quantity (often called stress):

S := Σ ( D_ij - d_ij)²

Now, from the view point of dimensional analysis the stress function in its above form is not consistent, since D_ij and d_ij have arguably different dimensions as they measure distance in different spaces.

One work-around to solve this inconsistency issue is to add a special constant or apply a monotonous function on d_ij like:

S := Σ ( D_ij - C·d_ij)² or S := Σ ( D_ij - f(d_ij))²

The constant C or the function f should magically convert the dimension of d_ij to that of D_ij. This kind of ad-hoc techniques look pretty artificial. Most scientists probably won't like this method.

It should be noticed that the RPM (relational perspecitve map) algorithm avoided this inconsistency by using the following stress function:

S := Σ D_ij/d_ij

The dimensional analysis rule does not prohibit composing terms of different dimensions by multiplication. In above formula, both sides of the equation have a completely new dimension, that is the dimension of D_ij divided by the dimension of d_ij.

Experiements with kernel PCA

2009-01-15T15:04:00.000-08:00

The newly released VisuMap software version 2.7 has added the kernel PCA (kPCA) service. From the user's perspective kPCA basically provides the same service as MDS methods, namely maping an abstract set of data points with a similarity distance to a multidimensional vector space, so that the Euclidean distance in the dimensional space somehow reflect the original similarity distance. Thus, the input for kPCA is a dataset with a similarity distance (or matrix), the output is a table of numerical values with one row vector for each data point.

As there are many effective analysis methods operating on dimensional vector space (like k-mean clustering, support vector machine), kPCA enables us, like MDS, to apply effective methods to a broader range of data.

In order see how kPCA works, I created a dataset with about 3000 data points which form together a sphere in the 3D space. With the Gaussian kernel (a way to measure similarity between data points) kPCA of VisuMap created a 50 dimensional datasets within about 100 seconds. That means, kPCA mapped a 3 D dataset into the 50 D space. The first 3 dimensions in the 50 D space basically mirror the original 3 dimensions. In order to see how other dimensions look like, I fixed the x- and y-axises of a 3D view window to the first two dimensions, then assign the z-axis to other dimensions one after the other while rotating the 3D view window. Then following small video clip shows how the spherical dataset looks like with those extra dimensions:

The first few seconds of the video shows how the first 3 dimension re-constructs the original sphere. Then, each time when the z-axis switched to another dimension, the sphere turned in to another geometrical shape. I am not sure how those geometrical shape get formed and how to take advantage of those new dimensions. But, those dimensions surely look interesting and a little mysterious.

K-Center and Affinity Propagation Clustering

2008-10-30T10:19:00.000-07:00

Recent release of VisuMap implements a new clustering algorithm called Affinity Propagation (AP). AP provides similar service as the metric sampling (MS) algorithm that has been available in VisuMap for a while. In this blog I will compare AP with MS from practical point of view with concrete samples.

First, a few words about the implementations of these algorithms. MS is basically a variation of the k-mean algorithm that uses the expectation maximization (EM) method to solve the k-medoid problem. Such algorithms are sometimes called k-center or k-medoid algorithm. Simple direct implementation of k-medoid algorithm often just find bad solutions (i.e. local minimums). People normally has to repeat the algorithm many many times to find an acceptable solution. The MS algorithm implemented in VisuMap uses a stochastic mechanism to avoid or reduce such local minimum problem.

AP uses message passing mechanism to find representative exemplars within dataset with a similarity structure. AP is basically a deterministic algorithm, it does not require random initialization or randomness during the optimization process. However, AP algorithm often suffers from the oscillation problem in the way that the algorithm wanders around several local minimums for quite long time. To reduce the oscillation, we can use a large dumping factor as suggested in the original paper of AP. But a more efficient way seems to be using small randomness when assigning preferences to the data points. By default, the implementation of AP in VisuMap enables this randomness.

Thus, both MS and AP in VisuMap are of stochasitc nature. I have applied the two algorithms each 500 times on the dataset SP500 (each data point contains the weekly price history of a stock in a year) and got all together 1000 different clustering solutions (each solution is a selection of 25 exemplar points). Using the t-SNE algorithm I created the map of these 1000 solutions as shown in the following picture:

In above map each red/green dot represents a solution found by AP/MS algorithm respectively. As a special characteristic of t-SNE algorithm, identical solutions will be arranged as one more circles. That means, dots in a rectangle in above map represent acutally identical solutions.

We notice first that the two algorithms found very different solutions as red and green dots are clearly separated from each other. There are only two small regions (marked by ovals) where the two algorithms found similar solutions. Secondly, we notice that MS found a large number of different solutions. AP found just few different solutions. This means that the AP optimization process often converged to common solutions.

The following table shows the statistics about the solutions found by AP and MS:

	Metric Sampling	Affinity Propagation
Min Error	0.4972	0.4901
Max Error	0.5058	0.5035
Average Error	0.5004	0.4960
Process Time in Seconds	5	4

In the above table the error of a solution is the average distance of the data points to their cluster center (i.e. exemplars). The average error is calculated over the 500 solutions found by each of the two algorithms. We can see that AP has found, on the average, better solutions than MS, but the improvement is only about 0.4%. So, it is fair to say that these two algorithms performed similarly with respect to this dataset.

The situation changed a lot when I applied the two algorithms to a different dataset, namely the dataset PharMine171 discussed in previous blog entry. This dataset contains about 1000 binary 171-dimensional vectors , whereas the SP500 dataset contains 500 40-dimensional vectors with real values. The following is the map of the clustering solutions and the corresponding statistics.

	Metric Sampling	Affinity Propagation
Min Error	2.1775	2.1101
Max Error	2.2321	2.1686
Average Error	2.1972	2.1603
Process Time in Seconds	6	30

From above map we can see that AP and MS still produced different set of solutions. AP also produced more different solutions than with the previous dataset. This means that, with respect to this dataset, AP is more sensitive to small variations/noise in the input data. This phenomenon might be explained as following: Because this dataset contains binary vectors, many similarities between data points might have exact equal values. Those equivalence of similarity will be broken different by small variations in the input data. Since AP indeed uses operation like max() to distinguish even infinitesimal differences, small variations may lead to significantly different solutions.

We also notice that AP is substantially slower than MS. For larger datasets, the performance difference will be even larger. For instance, for a dataset with about 4000 data points, MS is more than 100 times faster than AP. It should be pointed that both implementations of AP and MS uses the full similarity matrix. It remains to see how much speedup can be achieved by using sparse matrix.

More interestingly we see that the solutions found by AP is about 3.7% better than those found by MS; and even the worst solution found by AP is better than the best solution found by MS.

Now, one would ask the question: Does AP always find better solutions than MS? Unfortunately, this is not the case. The following picture shows the 3D presentation of the spheroid dataset that contains 1825 data points forming together a spheroid:

As be shown in the following test statistics, solutions found by MS is about 3.7% better than those found by AP:

	Metric Sampling	Affinity Propagation
Min Error	5.9029	5.9029
Max Error	6.4034	6.6483
Average Error	6.0948	6.3199
Process Time in Seconds	9	284

Conclusions:
I believe AP is a useful and promising addition to the arsenal of clustering algorithms. Not only that AP can find better solutions, it also finds different solutions as existing algorithms. As we know that there are many different clustering problems. AP doesn't seem to be designed to directly attack the k-medoid problem, but nevertheless it is competitive to those algorithms directly targeting the k-medoid problem. This indicates that AP may have potential to attack a much larger class of clustering problems.

Clustering clusters

2008-10-21T12:39:00.000-07:00

Many clustering algorithms, like k-mean algorithm, are of stochastic nature in the sense that they produce different partitions (of a data set) depending on the random initialization. In this blog I am going to use the mapping service provided by VisuMap to investigate the partitions generated by a particular clustering algorithm. In other word, I want to visualize the clusters among partitions generated by a clustering algorithm.

The clustering algorithm I am going to look at is the metrical sampling (MS) clustering algorithm which is basically a k-mean algorithm with some randomization extension. For a given set of data points with a distance matrix, MS algorithm tries to find a set of sample points, so that minimizes the average error, that is average distance of the data points to their nearest sample point (also called exemplar). Thus, MS algorithm is one of many available algorithms to solve the so called k-medoid clustering problem (which is an NP complete problem according to certain literature).

For a dataset of N data points, MS algorithm produces a N dimensional binary vectors. Each component is a binary value that decides whether a data point is selected as a sample point. So, when we run the MS algorithm many times with different random initializations, we will get a list of N dimensional binary vectors. I am interested to see clusters or patterns among those binary vectors. For this study, I picked up the dataset PharMine171 that contains 1025 fingerprints of 171 pharmacologically interesting compounds. I have run MS algorithm on the dataset for 1000 times. The following is a map of the 1000 partitions ( produced with the t-SNE mapping algorithm):

In above map, each spot represents a partition generated by one run of MS algorithm. Closely located partitions means they share more sample points. The coloring is done manually to highlight the clusters. We can easily see that there are 8 or 9 clusters among the 1000 partitions.

Now, one question that raises naturally is: Do the partitions in one cluster have similar qualities (i.e. average error)? In order to answer this questions I opened a bar view window that displays the average error of each partition in an ordered bar view as follows:

In above picture, each bar (in yellow or red color) represents a partition. The height of the bar indicates their average error. The red bars represent those partitions selected manually on the partition maps. So, when we select a cluster in the partition map, we can easily see the corresponding distribution of average error of that clusters. The following two pictures show two such distributions:

We notice that the average error in a cluster is spread all over the whole error range (from very bad to very good). This means, similar partitions can have very different qualities. Thus, the answer to above questions is negative.

More in general, I suppose this study indicates that any algorithm trying to find good solutions for the k-medoid problem need a mechanism to allow global variations in the searching space. The MS algorithm implemented in VisuMap does this through controlled randomness. Another interesting algorithm, named affinity propagation seems to do that with help of 2N² addition system states (called messages). Affinity propagation is also supported in VisuMap. I will do some comparison between MS and AP in another blog.

Feature selection by clustering

2008-09-13T08:09:00.000-07:00

I have blogged before that there is a "symmetry" between data clustering and dimensionality reduction: whereas the former reduces the number of rows; the latter reduces the number of columns in a data table. Thus, we can theoretically convert any clustering algorithm to a dimensionality reduction method by transposing the data table, and vise verse.

In this blog I am going to demonstrate how to use the metric sampling algorithm to reduce dimensionality with a concrete dataset. My sample dataset (PharMine171) contains data for 171 pharmacologically interesting compounds. Each data point contains 1025 binary vectors which indicate whether the compound contains particular fragments (ie. molecular sub-structurs). Thus, the dataset is a 171 x 1025 table with binary values and its dimension is 1025.

In order the reduce the dimensionality we first transpose the table to get a 171 dimensional dataset with 1025 data points (each data point now represents a fragment). In order to get feeling about the transposed dataset I have created 2 dimensional map for it (with VisuMap) as be shown in the following picture:

Map of 1025 binary features

Above map shows clearly some clusters among the 1025 data points. In order to find representative columns in our original dataset, we need to find representative rows in the new data table. In order to do so, we apply the metric sampling algorithm (available in VisuMap) that selects representative subset of data points. I have configured the algorithm to select 25 data points. The following picture shows how the 25 data points are distributed among the 1025 data points:

We can then export the 25 data points as a 25 x 171 table, and transpose the table to get an 171 x 25 data table. This table is then our new dataset with reduced dimension.

Now, how can we tell that the reduced dataset approaches the original dataset? We can do this pretty easily with VisuMap by creating maps for both datasets. If the two datasets contain about the same information, their maps should be visually similar to each other. The following picture shows the map of the original dataset (on the left side) and of the reduced dataset (on the right side). They are both created with the SMACOF mapping method. We see that the map of the reduced dataset reveals clearly the three clusters as the original dataset, but it shows less details within the clusters.

Notice that the metric sampling algorithm implemented in VisuMap is a variation of the widely known k-medoid algorithm. The recent release of VisuMap (version 2.6.807) utilized a randomization strategy in the optimization process that made the algorithm much more robust and less prone to local minimum problem; which has been a big issue for the k-medoid algorithm.

On similarity metrics for chemical compounds

2008-09-04T11:05:00.000-07:00

Recently, Yap Chun Wei has posted a dataset on the pharmine blog. The dataset consists of fingerprints of 171 pharmacologically interesting compounds. Just to recapture, the fingerprint of a compound is here a vector of 1025 binary flags, each flag in the vector indicates whether a particular molecular fragment is present in the compound. There are many ways to calculate fingerprints. Depending on the nature of the problem, you can use different algorithms or different collection of fragments. The mentioned dataset, for instance, used OpenBabel to calculate those fingerprints.

The dataset constains 3 different groups of compounds: penicillins, cephalosporins, fluoroquinolones. Using VisuMap Yap created different maps of these 171 compounds which showed, more or less, the cluster structure of the dataset. I personally find that the PCA map provides the best visualization. The following picture is a PCA map I created with VisuMap for this dataset:

Compound Map of 171 compounds

In above map, the 3 compound groups are displayed as glyphs in 3 different colors. The coloring are done manually with VisuMap based on known information about these groups. Although, for this dataset, you can get almost exactly these 3 clusters using the k-mean algorithm provided by VisuMap. The bar diagram in the picture shows the presence frequency of the 1025 fragments among the 171 selected compounds. That means, a higher bar indicates that a particular fragment is present in more compounds.

The above map visualizes the similarity information between the 171 compounds. That means, closely located compounds will have similar fragment collections and therefore similar pharmacological properties. The similarity information are basically encoded in the fingerprints. Thus, the method to calculate of those fingerprints is naturally critical for this kind of data analysis.

In order to better understand those fingerprints we can created a map of those 1025 features with VisuMap. In order to do so, we simply transpose the binary data table (via the menu Edit>Filter Data>Transpose Table in VisuMap), so that each binary feature becomes a new data point; and each compound become a feature in the transposed dataset. We can then pick a mapping algorithm and metrics to create a feature map. The following picture shows such a map created with the t-SNE algorithm and the tanimoto dissimilarity metric:

Feature map of 1025 binary features

Above picture shows 4 or 5 clear clusters on the left side represented by colored glyphs. The rest are more or less randomly distributed. It turned out that those yellow-square features are those fragments which are NOT present in any of those 171 compounds (all bits are zero). Therefore, they carry no direct information about our compounds. Interestingly, these zero vectors form together a homogeneous cluster in the map.

Other clusters in above map represent groups of fragments which have high frequency and are informative to distinguish the three compound groups in the original dataset. We can verify this with the help of the bar diagram in VisuMap as follows:

We first open the feature map in a separate window (via the menu Tools>Map Snapsot). Then open the compound map, and then select all compounds and open the bar view through the context menu "Bar View". The bar view by default displays the frequency in the order as given in the transposed data table. We then sort the bars through the context menu "Sort Values" so that bars are displayed in the order from low frequency to high frequency as depicted in the following two picture:

The sorted bar view shows 3 plateaus which correspond to clusters in the feature map and in the compound maps. In order to see the correspondence we select a plateau in the bar view with the mouse, the snapshot window of the feature map will automatically high light those selected features. As we can see in the following picture, the selected plateau of features clearly correspond to a particular cluster in the feature map (the marked cluster at lower left corner).

Correspondence between high frequency features and feature clusters

We notice that some of the clusters in the feature map show some fine sequential structure that may lead to more hints about the internal structure of the fragment collections.

With the knowledge about informative feature clusters we can, for instance, reduce the number of features significantly without significant loss of information about the clusters in the original dataset. The VisuMap dataset folder PharMine171.xvmz (zipped XML file) includes a reduced dataset with 298 features that characterizes similar similarity information as the original dataset with 1025 features.

The above VisuMap dataset folder also contains feature maps created with other mapping algorithms. It is interesting to notice that for the feature map dataset, the t-SNE algorithm provides the best visualization, whereas the result of the PCA algorithm is rather disappointing.

New mapping algorithm t-SNE

2008-08-13T08:56:00.000-07:00

Recently I spent some time to experiment with the t-SNE mapping algorithm proposed by L. Maaten and G. Hinton in Visualizing Data Using t-SNE. The algorithm is a variation of the SNE algorithm developed few years ago.

The t-SNE home site provides source code, technical papers and many samples, so that it is pretty easy to test the algorithm. The matlab source code is provided as a component in a matlab toolbox that also implements many other popular methods for dimensionality reduction. The web site also provides a binary executable that is optimized for the Intel/Windows platform. It took me a while to figure out its binary data format, but the time spent is worthwhile, as it runs much faster than the matlab version.

As demonstrated in its technical paper, t-SNE has the nice capability to show the cluster structure of high dimensional datasets in low dimensional maps. More interestingly, I noticed that t-SNE is able to segment dataset with complex structure in a natural way to preserve relevant information from certain perspective. For instance, the left map in the following picture shows a dataset that forms two crossing squares in 3D space, the map on the right side is a 2D map created with t-SNE. We notice that t-SNE segmented the connected shape into 5 clusters. The segmentation occurred roughly along the center line where the two squares crossing each other.

In general, I think there are two major challenges for mapping algorithms intended to visualize high dimensional complex data. The first one, I call it the singularity problem, is present in the above sample. The dataset is mostly a 2D surface except the center region where the two square cross each other. In order to map this dataset into 2D map in non-overlapping manner, an algorithm has to find a non-trivial way to segment the dataset along the singularity region. t-SNE did a good job here compared to other algorithms like CCA (curvilinear component analysis) which segments the dataset rather randomly.

The second challenge, I call it the closedness problem, occurs when the dataset forms certain kind of closed manifold in high dimensional space, like sphere, torus (e.g. surfaces of 3D objects). In this case, a mapping algorithm needs to have a mechanism to cut-off the closed structure underlying the dataset. The t-SNE algorithm behaviors also pretty well in addressing this kind of problems. The following video clip (produced with VisuMap 2.6.800) shows how t-SNE gradually deforms two overlapping spheres to a 2D map.

We notice that t-SNE first captures the global structure of dataset (like the Sammon algorithm or PCA), then deforms the map locally to achieve non-overlapping mappings. This kind of dynamical behavior is quite similar as the CCA algorithm.

The implementation of t-SNE is pretty straightforward, and it can be easily tuned for a wide range of data. Even with its current simple implementation, some of its features already surpassed many other popular mapping algorithms. So, I believe it will become a valuable tool for analyzing high dimensional data.

Villarceau Circles and Diagonal Rotation of RPM Map

2008-04-12T20:24:00.000-07:00

It has been known for long time that there are 4 perfect circles passing every points on the a torus surface. We can get two of these circles buy cutting the torus vertically and horizontally through that points. These two circles are the normal circles most people can easily visualize. The other two, not so apparent, circles are the so-called Villarceau circles. They are made by cutting the torus diagonally as depicted in the following animation (IE7 users need to enable the flag Tools>Internet Options>Advanced>"Play animations in web pages" to see the animation):

Figure 1: Creating the Villarceau Circles on a torus.

Above animation is very helpful to visualize the Villarceau Circles. Another simple way to show the Villarceau circles is using the fundamental polygon representation where a torus is represented as a rectangle whose opposite edges are glued together. In the polygon represetnation, the horizontal and vertical edges are the two normal circles (horizontal and vertical circles). The two diagonals are the Villarceau circles as depicted in the following diagram:

Figure 2: Four circles on a torus passing through a point.

We know that RPM algorithm simulates a multi-particle dynamic system confined on the flat torus surface. For some unknown reasons, RPM maps often form symmetry patterns alone the diagonal lines, e.g. along the Villarcea circles. For instance, when we map the sphere to the tours, RPM algorithm cuts the sphere into two equal sized half and map them as two discs to the torus surface. The two discs are positioned alone a a diagonal line as depicted in the following picture:

Figure 3: RPM Map of a sphere

We see in above picture that only one disc is completely displayed as a whole, the other one has been cut by the edges into 3 or 4 pices (which will form a complete disc when we glue the opposite edge together). The whole picture is more or less symmetric with respect to the left-up-to-right-bottom diagonal line. We can get a better visualization of above image when we tile multiple copies of above map together as be show in the following picture:

Figure 4: Tiling duplicates of Figure 3.

Now, it would be helpful for the purpose of visualization if we can show the two discs as whole discs in a single map without using the tiling techniques. We can achieve this by transforming the coordinator from the two vertical and horizontal edges to the two diagonals (i.e. the two Villarceau circles). The transformation is schematically depicted in the following diagram:

Figure 5: Rotation to the diagonal circles.

We call the above transformation the diagonal rotation as it basically rotate the coordinator to the diagonals. Formally, the diagonal rotation is done by the following forms:

In above forms w and h are the width and height of the fundamental rectangle respectively. d is the length of its diagonal. The above forms transforms points from the region D to D'. For points from other regions, we have to apply some simple shifting after the transformation as indicated in Figure 5.

As indicated in Figure 5, the two diagonals will become the edges of the transformed map. Thus, in order to avoid the resulting map be cut off by the edges we should shift a RPM map away from the diagonals through torus cyclical shifting. For our sphere example we can shift the map, for instance, to the following position:

Figure 6: Sphere map shifted off from the diagonal lines.

When we now apply the diagonal rotation on the map in Figure 6 we will get the following rotated map:

Figure 7: The sphere map after diagonal rotation.
Comparing with the original RPM map in Figure 3, the rotated map in figure 7 is obviously easier to investigate for human eyes.

It should be pointed out that a RPM map after diagonal rotation is no long a torus map in the sense that you cannot tile the map to fill 2D space the same way as be shown in Figure 4. Since this will break the edge-connectivity as indicated in Figure 5. But you can fill the 2D space with the rotated RPM map like building a bridge wall: you tile the map site by site to build a continuous lay, then build next lay the same while shift the lay to the left (or right) for half block. The following picture depicts this characteristics:

Figure 8: Tiling with rotated map.

The Equivalence of Cooling and Expanding Dynamics

2008-04-04T17:00:00.001-07:00

The RPM (relational perspective map) algorithm has a seemly add-hoc measure that gradually decreases the learning rate to force the convergence of the algorithm. I have tried to avoid such measure with various methods (e.g. conjugated algorithm exploring second order gradient), but was not able to find anything working.

I have notice a significant problem here: In the original Newton-Raphson algorithm, the algorithm should converge to a local minimum if the learning rate is sufficiently small. But, this not so for the optimization problem of RPM. The algorithm never converges, no matter how small the learning rate is (in a way that is doable with double precision arithmetic). If I keep the learning rate constant at certain small value, the simulated dynamic system goes in to an equilibrium state with constant variance.

Thus, it looks like NR method, as it is, does not work for massively complex dynamic system, it is not stable enough. Letting the learning rate gradually vanish seems to be the simplest method to enforce the convergence.

Mathematically, a very fine granulated gradient descent algorithm should be able to find a local minimum. But, I think this is practically not possible. To understand this, we notice that the learning rate practically controls the average variance of the positions of the simulated particles. That means, the learning rate is a kind of "temperature" of the system. In physics, we know that if we keep the temperature constant, the system will never converge to a frozen state. We can only reach that state when we gradually reduce the temperature to zero.

It is interesting to notice that the cooling dynamic system simulated by RPM is theoretically equivalent to a n-body dynamic system that is ever expanding. To see this, we notice that the learning rate is directly linked to the average variance of the particle's positions. Thus, instead of gradually lowering learning rate, we can just reduce the scale of unit length of the space compared to total size of the space (e.g. the size of the torus). The latter is again equivalent to a system in which we keep the scale of unit length constant, but gradually expand the size of the whole space.

We see here certain similarity of RPM model with the inflationary big-bang model from cosmology. Whereas RPM keeps the total size of the system stationary and gradually reduces the energy, the big-bang model keeps the total energy constant but let the universe expands. It is probably not a too crazy idea to consider a dual big-bang model that keeps the apparent universe stationary, but magically reduces the total energy/mass/temperature.

The "Big Three" Data Categories

2008-03-10T08:28:00.000-07:00

Data analysis is an essential part of science and engineering. As a guide to deal with various data analysis tasks, I find it helpful to group data into three basic categories:

Multidimensional tabular data (e.g . tables of numbers and text).
Time series data (e.g. image, audio/video data).
Graphics oriented data (e.g. tree and network structured data).

All data analysis methods can be grouped in to three categories depending on which data category they target. So, for instance, wavelet transformation is a method for time series; whereas MDS are methods for multidimensional data. Context free grammar and text parsers are mainly tools for graphics oriented data.

Furthermore, all software applications can also be grouped in to three corresponding categories. A software package that supports multiple data categories (like java system, Matlab) is either a generic tool that requires further software for end-users; or it is hopelessly difficult to design, implement, understand and to use.

When I encounter a new data analysis problem I first identify its category, then I start to search for methods and software packages targeting that category. I would not try, for instance, to use MDS or SQL to analyze images or video/audio data.