Data Visualization & Clustering: Generative Topographic Mapping

Downloaded by UNIV OF CALIFORNIA SAN DIEGO on December 12, 2016 | http://pubs.acs.org Publication Date (Web): October 5, 2016 | doi: 10.1021/bk-2016-1222.ch010

Chapter 10

Data Visualization & Clustering: Generative Topographic Mapping Similarity Assessment Allied to Graph Theory Clustering Matheus de Souza Escobar, Hiromasa Kaneko, and Kimito Funatsu* Department of Chemical System Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8654 *E-mail: [email protected]

Chemical systems can be discriminated in several ways. If one considers industrial data, process monitoring of chemical processes can be achieved, with the following applications: anomaly discrimination (fault detection) and characterization (fault diagnosis & identification). This chapter presents an unsupervised methodology for data visualization and clustering combining Generative Topographic Mapping (GTM) and Graph Theory (GT). GTM and its probabilistic nature highlights system features, reducing variable dimensionality and calculating similarity between samples. GT, then, generates a network, clustering samples, normal and anomalous, according to their similarity. Such assessment can be applied, however, to other data sets, such as the ones involved in drug design and discovery, focusing on clustering of molecules with similar characteristics. Two case studies are presented: a simulation data set and Tennessee Eastman process. Principal Component Analysis (PCA), Dynamic PCA and kernel PCA indexes Q and T2, along GTM independent monitoring methodologies are used for comparison, considering supervised and unsupervised approaches. The proposed method performed well for both scenarios, revealing the potential of GTM and network based visualization and clustering.

© 2016 American Chemical Society Bienstock et al.; Frontiers in Molecular Design and Chemical Information Science - Herman Skolnik Award Symposium 2015: ... ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


Introduction It is fair to say, at least in the realm of chemical systems, that machine learning techniques are used with the following key element in mind: highlighting important features and characteristics of a process by analyzing the relationship of variables and samples belonging to it. This notion, as general as it may be, leads to distinct applications, such as pattern recognition, predictive modeling, classification, clustering, etc. On a more particular note, however, the work explored in this chapter relies on evaluating similarity between samples in a data pool through data visualization and clustering techniques. The spark originating the topics discussed in this chapter lies on industrial data, particularly time-series measurements of chemical plants, where process monitoring judges whether for any system its current conditions are normal or abnormal. This can be divided in two main tasks to be accomplished: fault detection and diagnosis. The former is concerned about finding anomalies, discriminating them from the normal states in the system (1). The latter tries then to pinpoint which variables are more likely to be responsible for the anomaly (2). Such fault assessment is fundamental for process optimization, operation and safety. Broadly speaking, however, this assessment reflects topics on data visualization and clustering, whose applications are far wider. One might imagine, for example, that molecular similarity and design applications are too distant from such anomaly strategies. The techniques explored here, however, are broad enough to motivate their use in different fields, which will be hinted throughout the text. Instead of anomaly assessment, for example, clustering of similar molecules is possible using the same techniques displayed here. Another important aspect to have in mind is whether such assessment is performed relying on preexistent labels (3, 4), which might be reliable or not, or if it should rely only on data and, consequently, on the correlations between variables and samples. From an industrial standpoint, labeling, also called a supervised approach, is time-consuming, expensive and whose reliability depends on current plant operational conditions. From a chemical space perspective, find such labels can be also overwhelmingly complex. Once one considers the size of the chemical space and how many of the compounds in it have never been synthetized, labeling seems farfetched. In addition, depending on which chemical property is being targeted, technological limitations might hinder such evaluation. Relying on what data has to offer, an unsupervised strategy therefore, may lead to results that are more reliable and long-lasting. As for applications in the industry, database maintenance, for example, can exclude irrelevant data by creating, from a given data set, a detection model that can identify normal and abnormal samples for future online evaluation. Such technique can, then, improve soft sensors’ (5) accuracy, by using only reliable samples for model generation. The work presented here, despite its application for anomalous scenarios, reveals an aspect of data visualization and clustering that goes beyond such assessment. Instead of anomalous and normal samples, one can isolate any data set into clusters of samples whose similarity is high. When it comes to molecular 176 Bienstock et al.; Frontiers in Molecular Design and Chemical Information Science - Herman Skolnik Award Symposium 2015: ... ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


design, one could, for instance, group structures with similar characteristics and, relying on the web of connections given by the network representation, evaluate the relationship between distinct clusters that share unique traits. Works in the literature play with some of these concepts, particularly molecular similarity. Dimensionality reduction techniques, for instance, have been used several times for such evaluation (6). Approaches such as Principal Component Analysis (PCA) (7), Self-Organizing Maps (SOM) (8) and Generative Topographic Mapping (GTM) (9) are but a few of the ones available. The main goal for such approach is to evaluate and visualize how variables and samples interact with the system and with each other. From this premise, a combined data visualization and clustering approach is presented, which for industrial data characterizes a nonlinear process monitoring tool. By doing so, we hope to achieve a more complete understanding of any system, where data discrimination is more objective and meaningful. When developing such methodologies, the quality of the information available is fundamental for the development of trustworthy models. Real data sets struggle with redundant information and noise, which might hide the true relation between different features and therefore, different samples. Dimensionality reduction, presented through techniques whose data visualization is equally important, identifies regions with similar characteristics and filters irrelevant information from data. One of the most widespread methods for process monitoring is Principal Component Analysis (PCA) (10), which assesses linear correlation between different process variables, so to reduce the dimensionality of highly correlated variables. Its use is so widespread that numerous PCA-based extensions were devised, such as dynamic PCA (DPCA) (4), recursive PCA (11), distributed PCA (12) and maximum-likelihood PCA (13). Such extensions, nonetheless, do not deal directly with process nonlinearity, as kernel PCA (14) do. Other methods also tackle nonlinearity from scratch, such as Support Vector Machines (SVM) (15), Gaussian Mixture Models (GMM) (16), Generative Topographic Mapping (GTM) (17) and even the use of inferential models (18). The main element explored in this and previous works (19) is GTM. Its nonlinear and probabilistic nature leads to a better handling of complex and realistic scenarios. When it comes to similarity assessment, each sample plotted in GTM’s latent space has a unique probability distribution (PD), a fingerprint, associated to each latent grid point. By assuming that samples with correlated PD profiles represent data with similar characteristics, GTM can be used for fault detection and dimensionality reduction simultaneously, promoting discrimination of normal and anomalous data. The clustering element relies on GT (20), where similarity information is used for the establishment of a network. Then, its density and number of connections unravel clusters with different characteristics. The developments presented here reveal different features of this combined approach, allowing better refinement of normal clusters and revealing a myriad of different interpretations for the networks established. Two case studies are defined for performance comparison. Initially, a simulation data set with multiple anomaly scenarios is presented. Secondly, the Tennessee Eastman Process (TEP) (21) is considered for validation of 177

Bienstock et al.; Frontiers in Molecular Design and Chemical Information Science - Herman Skolnik Award Symposium 2015: ... ACS Symposium Series; American Chemical Society: Washington, DC, 2016.

the methodology. The proposed method (GTM+GT) is compared against unsupervised and supervised PCA, DPCA, KPCA and GTM independent approaches (3, 4).


Data Visualization & Dimensionality Reduction The assessment of normal and abnormal states relies greatly on calculating similarity between all elements in the data pool. Data visualization, which also implies dimensionality reduction, can process data in a way to highlight samples’ characteristics while reducing redundant and unnecessary information that might cloud the true relationship between samples. Several techniques are commonly used for such assessment, where PCA is the most straightforward linear approach known, relying on variables being converted into linearly uncorrelated variables called principal components (PC), through an orthogonal transformation (10, 22). The main strategy selects only PCs whose accumulated component contribution is just below 99%. Such selection is also important for assessing future samples, considering that Q and T2, the most relevant PCA monitoring indexes (23), derive from it. Since one of its main limitations is its inherent linear nature, however, its application for more complex, nonlinear systems are limited. From this premise, pseudo-nonlinear and nonlinear extensions, such as DPCA and KPCA respectively, were developed as well. The main approach considered here, however, relies on GTM. GTM is a widely used technique applied for visualization of high dimensional data. It consists of a probabilistic nonlinear approach, where a low-dimensional latent variable z is represented in a 2D space, aiming to approximate original data x as a high-dimensional manifold on the original data space. This manifold is modeled by a Gaussian function. Acting as a bridge between spaces, an intermediary layer of radial basis functions (RBFs), also Gaussian, is created (17). RBFs are embedded in a mapping function y(z;W), which defines the non-Euclidean manifold and connects both spaces. Figure 1 shows the schematic representation behind GTM.

Figure 1. GTM overall concept representation (19). 178 Bienstock et al.; Frontiers in Molecular Design and Chemical Information Science - Herman Skolnik Award Symposium 2015: ... ACS Symposium Series; American Chemical Society: Washington, DC, 2016.

GTM Structure


The main goal behind GTM is to find a representation for the distribution p(x) of data in a D-dimensional space x = (x1, ..., xD) associated to a number L of latent variables z = (z1, ..., zL). A function y(z;W) is devised, mapping points z in the latent space into the equivalent y(z;W). The transformation y(z;W) maps the latent variable space into and L non-Euclidean manifold S embedded within the data space. W is a parameter matrix that governs the mapping from z to x. The distribution of x is chosen, for a given z and W, to be a radially symmetric Gaussian centered on y(z;W) having variance β−1, as shown in Equation 1.

The distribution in x-space is then obtained by integration over the z-distribution, assuming a known value for W, according to Equation 2.

where p(z) is the prior distribution of z. Once a data set of N data points X = (x1, ..., xN) is given, the unknown parameters W and β can be optimized, using maximum likelihood. It is more convenient, though, to maximize log likelihood, as presented in Equation 3.

One problem with this representation, however, is that despite specifying p(z) and the functional form of y(z;W), the integral specified in Equation 2 is usually analytically intractable. To circumvent this issue, y(z;W) is chosen to be a linear function of W and p(z) has to be defined accordingly. One option is to define p(z) as Gaussian, then the integral becomes a convolution of two Gaussians. In this case, however, the model is closely related to PCA, where the maximum likelihood solution for W columns leads to scaled principal eigenvectors. In order to expand this formalism to nonlinear y(z;W) functions, p(z) has to be defined in a specific form, as shown is Equation 4.

where G is the number of nodes in latent space assuming a regular grid. p(z) is given by a sum of delta functions centered on nodes in a latent space grid. This implies that probability distribution is local in each point of the lattice and not continuously distributed along the latent space. The x-distribution function now takes a different form from Equation 2, as presented in Equation 5.

and the log likelihood function is now given by Equation 6

179 Bienstock et al.; Frontiers in Molecular Design and Chemical Information Science - Herman Skolnik Award Symposium 2015: ... ACS Symposium Series; American Chemical Society: Washington, DC, 2016.

Expectation-Maximization Algorithm


This structure can now be optimized for W and β, once y(z;W) is defined. Knowing that the model developed consists of a mixture distribution, the Expectation-Maximization (EM) Algorithm might be the most suited for optimization (24). This algorithm relies on a suitable choice of y(z;W), such as a generalized linear regression model as described in Equation 7.

where ɸ(z) consists of B fixed basis functions ɸi(z), and W is a parameter matrix D x B relating these functions with the non-Euclidean manifold S. For a large class of basis functions, Radial Basis Functions (RBF) are universal approximators (25). These structures, particularly the Gaussian RBFs, are interesting, due to their fast training. GTM training using multi-level programming (MLP), for example, is prohibitive (26). Once the basis function structure is defined, the optimization can be executed. In the expectation step Wold and βold, two parameters assumed to be the current values of W and β, are used to evaluate the posterior probabilities, also called responsibilities, of each Gaussian component i for every data point xi using Bayes’ theorem, as shown in Equation 8.

This leads to the expectation of the log likelihood data presented in Equation 9.

Wnew and βnew can then be obtained on the maximization step, by maximizing Equation 9 with respect to both parameters independently, as shown in Equations 10 and 11.

This cycle of expectation and maximization is repeated until the objective function reaches a maximum, according to a satisfactory convergence. Data Visualization and Latent Probability Distribution Once the map is trained, it is possible to determine for each sample the likelihood of it belonging to each node in the latent grid, establishing a PD profile. The profile comes from the responsibility matrix obtained from the optimization procedure aforementioned, as suggested by equation 12. 180 Bienstock et al.; Frontiers in Molecular Design and Chemical Information Science - Herman Skolnik Award Symposium 2015: ... ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


Such profiles are represented as individual heat maps or as one cumulative plot for all data, as Figure 2 suggests. If we consider that the combination of variables in hyperspace are generally different for each sample, it is possible to infer that PD profiles are unique. Such assumption allows the similarity assessment of any sample against other ones. The overall structure of the proposed methodology presented in this chapter relies on calculating similarity, a requirement that GTM elegantly fulfills.

Figure 2. Representative GTM PD heat map for a) one sample and b) an entire corresponding data set containing samples with distinct feature values (19).

Data can also be visualized by collapsing PD profiles into mean and mode plots for all samples. Since each datum can be represented as a dot, visualization of distinct clusters in the data set might be more apparent. It is important to notice that since data is collapsed, though, information is lost and such assessment is, therefore, not entirely reliable. Figure 3 shows a comparison between PD heat map and mean/mode GTM plots. Complementary to that, one should be aware of how data is dealt within the optimization algorithm. GTM considers all samples to be independent, identically distributed vectors (i.i.d.), which implies that dynamic information is not being considered. If time-series data is not being used, such as structural data sets used in molecular design, such premise has no impact in the final map. For industrial monitoring applications, however, not using dynamic information is, at most, a waste of valuable information in one’s data set. One way to circumvent such issue is to consider time-delayed variables, analogous to DPCA data treatment, so to incorporate existent dynamic information. 181



Figure 3. a) GTM PD heat map and b) mean/mode plot. Different colors in the mean/mode plot represent different process states, where dots are collapsed PD means and circles are collapsed PD modes.

On GTM Hyperparameters One aspect of GTM that deserves special attention is how to set the hyperparameters, which are structural parameters defined previously to optimize the parameters W and β mentioned earlier. GTM relies on the following set of hyperparameters for its utilization: latent grid size, number of RBFs, width of RBFs and regularization parameter λ. The optimal value for each parameter is usually determined via exhaustive search, using cross-validation to look for the minimization of reconstruction error, i.e., distance from the manifold, once data is recreated into the original hyperdimensional space. Root Mean Squared Error (RMSE) is usually used as an index for such assessment, as described in Equation 13.

where N is the number of samples and M is the number of variables. xi,j is the original ith sample value for the jth variable and is the respective remapped value. Regular RMSE, however, does not take into account certain factors, such as the smoothness of the map, which allied to poor choosing of hyperparameters might lead to overfitting, a serious problem in GTM. As pointed out by several works in the literature (27–29), GTM overfitting is often overlooked. Similarity assessment is greatly impacted if the map is overfitting. Firstly, overfitting leads to samples with concentrated PD in the latent grid, reducing the likelihood of high similarity between samples. Secondly, new samples cannot be reliably incorporated to the new map without full re-training, which is time-consuming and inefficient. Root Mean Squared Error of Midpoint (RMSEM) tackles this issue (30), where midpoints to those existent in training data are used for accuracy assessment. If those samples can be predicted accurately, then not only training 182


data has high prediction accuracy, but also the regions in between, preventing overfitting and concentrating sample’s PD. RMSEM is calculated according to Equation 14.


where L is the number of midpoints selected,

is the ith midpoint value for

is the respective remapped point. Midpoints are the jth variable and sampled randomly from all possible combinations of training data, usually in a greater number (five-fold) than the original data set. When it comes to monitoring, the discrepancy between original samples and remapped ones can be used for monitoring, assuming that different samples have well defined Remapping Errors (RE) and that clusters have a RE threshold, according to 99% of the maximum RE in that cluster. Besides RE, similarity between probability distributions can be used as an index. The same matrix used for the network generation can be used for assessment. From an unsupervised perspective, little knowledge can be obtained about the process, since there is no reference. If normal data is known, however, it is possible to see whether external samples are at least 99% similar to at least one of the samples in the normal data pool. Assuming that all normal samples are part of the same state, if any external sample is similar to any sample in the pool, it can be considered normal as well. This idea leads to a unique threshold for each sample, since maximum similarity for each query sample against the normal data pool is different.

Graph Structure & Clustering Graphs are symbolic representations of networks that model pairwise relations between objects (31). For practical purposes, all graph-related structures presented will be referred as graphs. There are two basic elements for every graph: nodes and edges. The former represents observations (samples) and the latter indicates connections between those observations. For a given data set, adjacency matrix (AM) formalizes this web of connections, by representing all connections via a square matrix whose size is directly related to the number of observations available. Figure 4 shows an example of such representation. All null values show that there is no connection between respective pair nodes. Values different from zero, on the other hand, reveal links between nodes, where the strength of the connection is correlated to the respective adjacency value. AM is the core element of any graph, from where graph analysis, visualization and clustering is possible. This matrix is replaced by the similarity matrix obtained through GTM’s similarity assessment for network construction. 183



Figure 4. Schematic representation of a a) weighted AM and b) its respective undirected graph.

Louvain Community Finding Algorithm When it comes to graph clustering, various strategies are available. Spectral Partitioning (SP) (32) and Girvan-Newman Algorithm (GNA) (33) are among the most straightforward ones, for example. SP and GNA, however, have no termination criterion for optimal clustering and GNA relies on betweenness (20), a graph centrality measure for finding important hubs in the graph, which may not be available for a given graph or it is computationally too expensive to calculate. In order to cope with these limitations, Louvain Community Finding (LCF) (34) presents itself as an algorithm with intriguing features, based on, generally speaking, evaluating the density of edges within a group via an index called modularity (35, 36). LCF algorithm has two steps: local modularity optimization and graph update. Initially, a weighted graph of N nodes is created, where different clusters are assigned to each node, i.e., there are as many clusters as nodes. From this framework, a maximization of modularity is trailed, following the pseudo-algorithm below. 1. 2.

3.

For each node i, consider all neighboring communities j of i. Compute modularity gain (ΔQi,j) when i moves to each community j. i moves to the cluster with maximum gain, only if the gain is positive. Otherwise, i stays in its original community. Figure 5 shows the schematic representation of step 2 for one node, when tested against three other communities. Test modularity gain for all nodes in sequence, till no further improvement is encountered. Modularity gain is described in Equation 15.



where m is the total sum of edge weights in the graph, ki,j is the sum of edge weights from i to j, kI is the sum of edge weights incident to i, Σtot is the sum of edge weights incident to nodes in j and λ is called the resolution limit, regulating both terms of Equation 15. Lower λ results in fewer clusters, where higher λ results in more clusters. Once the algorithm stabilizes, the graph is updated by condensing all nodes belonging to a single community into a single node, keeping in mind that edges between nodes of the same community lead to self-loops. After the update, all steps above are repeated until no more modularity gain is achieved. Figure 6 and Figure 7 show graph and modularity evolution during LCF cycles for a trivial example.

Figure 5. Modularity gain test, where different background patterns indicate different communities (19).

Figure 6. Graph evolution according to LCF algorithm.



Figure 7. Modularity evolution during LCF cycles (19). For this simple example, the graph reaches a modularity peak after the second cycle, indicating the optimal configuration. By observing the original graph on cycle 0, one can easily see that it corresponds indeed to the best clustering scenario. Any further clustering beyond that results in a modularity drop. LCF uses modularity in a similar way that Girvan-Newman algorithm, but it provides an intuitive algorithm with a clear termination criterion.

GTM+GT Combined Approach The main methodology explored in this work involves two key elements: extraction of data essential information and effective data clustering. This is achieved by combining GTM and GT (37). Initially, GTM reduces data to a 2D latent plot, removing redundant and irrelevant information from the original data set. Every sample in the latent space has a unique PD profile, which is used for similarity assessment, as represented schematically in Figure 8 for two responsibility vectors r1 and r2.

Figure 8. Correlation assessment between two samples using the same GTM grid. Each sample’s PD can be expanded in a vector, which then is used for squared Pearson product-moment correlation coefficient (r2) calculation. Each assessment between samples fills one element of the AM. Once all samples are cross-evaluated, AM construction is finished. 186 Bienstock et al.; Frontiers in Molecular Design and Chemical Information Science - Herman Skolnik Award Symposium 2015: ... ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


With the AM built, LCF can cluster data into groups with similar characteristics. For unsupervised fault identification, it is assumed that faults are a minority of the system and due to their faulty nature, their behavior is usually more erratic, i.e., less stable. Normal operational data, on the other hand, represents generally a majority of the samples available, where data itself is stable. From a graph theory perspective, this means that normal data has a far higher number of connected nodes combined with higher connection density, which is used as reference for identifying the optimal NC. It is also important to notice that anomalous data might be detected as not one cluster, but several, representing different fault characteristics or different states within one fault development over time.

AM Construction The merit of GTM is to extract only relevant information for assessing similarity between data. How to calculate it properly, however, becomes the challenge. It was presented earlier that r2 was used for calculation. The reason behind this assumption, however, was not clarified. PD profiles are analogous to images, where comparing two images relies on evaluating the discrepancies in pixels. From image processing, similarity can be defined according to structural discrepancies, also called Structural Similarity (SSIM) (38). It relies on three distinct features: luminance, contrast and structure, described by the Equations 16-18 (39), respectively, for two PD (responsibility) vectors r1 and r2.

where σ is standard deviation, r̅1̅ and r̅2̅ are the average values of r1 and r2 and Ci are arbitrary constants to avoid instability when the denominator is very close to zero. is shown in Equation 19.

Remembering that G is the number of nodes in latent space assuming a regular grid. Luminance considers differences in the average PD value, contrast compares variance changes in PDs and structure calculates the correlation between PDs. When two images are compared to assess degradation, for example, all those three elements are important and equally relevant. For PD evaluation, however, luminance and contrast are far less important than structural comparison. Figure 9 shows the usual PD for two samples in a GTM map. 187 Bienstock et al.; Frontiers in Molecular Design and Chemical Information Science - Herman Skolnik Award Symposium 2015: ... ACS Symposium Series; American Chemical Society: Washington, DC, 2016.


Figure 9. GTM PD for two samples in a trained map. When images are compared, all pixels have important information and therefore all three indexes are relevant, trying to extract all tiny nuances in both images. GTM PDs are much more crude and limited in their representation, occupying specific regions of the map. Dissimilar samples will occupy different regions in the map, resulting in low similarity. Similar samples will overlap and be somewhat similar. This structural assessment depicted in Equation 17 is enough for our applications. It is important to notice that given some mathematical manipulation, Equation 17 corresponds directly to r2. Assessing similarity between all samples in the system, then, is being portrayed as simple r2 calculation. Depending on the map size chosen for GTM, however, each PD vector can easily go to thousands of variables. Even for small map sizes such as 10 x 10, there are 100 features being compared between samples for similarity assessment. Knowing that samples have local PD, calculating solely correlation between samples is not enough, since it would result in very low similarity for any samples marginally different. In order to cope with that, good practices recommend local assessment, by creating a moving window, which slides point by point throughout both latent grids. Once all local values are calculated, an average similarity is calculated and integrated to A, the AM, as shown in Equation 20 for r1 and r2 presented previously.

where r1j and r2j are local vectors and W is the number of local windows. Once similarity assessment is finished for all responsibility combinations, AM can be constructed. Finally, for any AM very low correlation values are recurrent (

Data Visualization & Clustering: Generative Topographic Mapping

Recommend Documents