Visualization and interpretability in probabilistic dimensionality reduction models
Abstract
Over the last few decades, data analysis has swiftly evolved from being a task addressed mainly within the remit of multivariate statistics, to an endevour in which data heterogeneity, complexity and even sheer size, driven by computational advances, call for alternative strategies, such as those provided by pattern recognition and machine learning. Any data analysis process aims to extract new knowledge from data. Knowledge extraction is not a trivial task and it is not limited to the generation of data models or the recognition of patterns. The use of machine learning techniques for multivariate data analysis should in fact aim to achieve a dual target: interpretability and good performance. At best, both aspects of this target should not conflict with each other. This gap between data modelling and knowledge extraction must be acknowledged, in the sense that we can only extract knowledge from models through a process of interpretation. Exploratory information visualization is becoming a very promising tool for interpretation. When exploring multivariate data through visualization, high data dimensionality can be a big constraint, and the use of dimensionality reduction techniques is often compulsory. The need to find flexible methods for data modelling has led to the development of non-linear dimensionality reduction techniques, and many state-of-the-art approaches of this type fall in the domain of probabilistic modelling. These non-linear techniques can provide a flexible data representation and a more faithful model of the observed data compared to the linear ones, but often at the expense of model interpretability, which has an impact in the model visualization results. In manifold learning non-linear dimensionality reduction methods, when a high-dimensional space is mapped onto a lower-dimensional one, the obtained embedded manifold is subject to local geometrical distortion induced by the non-linear mapping. This kind of distortion can often lead to misinterpretations of the data set structure and of the obtained patterns. It is important to give relevance to the problem of how to quantify and visualize the distortion itself in order to interpret data in a more faithful way. The research reported in this thesis focuses on the development of methods and techniques for explicitly reintroducing the local distortion created by non-linear dimensionality reduction models into the low-dimensional visualization of the data that they produce, as well as in the definition of metrics for probabilistic geometries to address this problem. We do not only provide methods only for static data, but also for multivariate time series. The reintegration of the quantified non-linear distortion into the visualization space of the analysed non-linear dimensionality reduction methods is a goal by itself, but we go beyond it and consider alternative adequate metrics for probabilistic manifold learning. For that, we study the role of \textit{Random geometries}, that is, distributions of manifolds, in machine learning and data analysis in general. Methods for the estimation of distributions of data-supporting Riemannian manifolds as well as algorithms for computing interpolants over distributions of manifolds are defined. Experimental results show that inference made according to the random Riemannian metric leads to a more faithful generation of unobserved data.
Full text
Visualization and Interpretability in Probabilistic Dimensionality Reduction Models Alessandra Tosi Department of Software Universitat Politècnica de Catalunya This dissertation is submitted for the degree of Doctor of Philosophy 2014
Abstract Over the last few decades, data analysis has swiftly evolved from being a task addressed mainly within the remit of multivariate statistics, to an endevour in which data heterogeneity, complexity and even sheer size, driven by computational advances, call for alternative strategies, such as those provided by pattern recognition and machine learning. Any data analysis process aims to extract new knowledge from data. Knowledge extraction is not a trivial task and it is not limited to the generation of data models or the recognition of patterns. The use of machine learning techniques for multivariate data analysis should in fact aim to achieve a dual target: interpretability and good performance. At best, both aspects of this target should not conflict with each other. This gap between data modelling and knowledge extraction must be acknowledged, in the sense that we can only extract knowledge from models through a process of interpretation. Exploratory information visualization is becoming a very promising tool for interpretation. When exploring multivariate data through visualization, high data dimensionality can be a big constraint, and the use of dimensionality reduction techniques is often compulsory. The need to find flexible methods for data modelling has led to the development of non-linear dimensionality reduction techniques, and many state-of-the-art approaches of this type fall in the domain of probabilistic modelling. These non-linear techniques can provide a flexible data representation and a more faithful model of the observed data compared to the linear ones, but often at the expense of model interpretability, which has an impact in the model visualization results. In manifold learning non-linear dimensionality reduction methods, when a high-dimensional space is mapped onto a lower-dimensional one, the obtained embedded manifold is subject to local geometrical distortion induced by the non-linear mapping. This kind of distortion can often lead to misinterpretations of the data set structure and of the obtained patterns. It is important to give relevance to the problem of how to quantify and visualize the distortion itself in order to interpret
iv data in a more faithful way. The research reported in this thesis focuses on the development of methods and techniques for explicitly reintroducing the local distortion created by non-linear dimensionality reduction models into the low-dimensional visualization of the data that they produce, as well as in the definition of metrics for probabilistic geometries to address this problem. We do not only provide methods only for static data, but also for multivariate time series. The reintegration of the quantified non-linear distortion into the visualization space of the analysed non-linear dimensionality reduction methods is a goal by itself, but we go beyond it and consider alternative adequate metrics for probabilistic manifold learning. For that, we study the role of Random geometries, that is, distributions of manifolds, in machine learning and data analysis in general. Methods for the estimation of distributions of data-supporting Riemannian manifolds as well as algorithms for computing interpolants over distributions of manifolds are defined. Experimental results show that inference made according to the random Riemannian metric leads to a more faithful generation of unobserved data.
Abstract (Catalan) Durant les últimes dècades, l’anàlisi de dades ha evolucionat ràpidament de ser una tasca dirigida principalment dins de l’àmbit de l’estadística multivariant, a un endevour en el qual l’heterogeneïtat de les dades, la complexitat i la simple grandària, impulsats pels avanços computacionals, exigeixen estratègies alternatives, tals com les previstes en el Reconeixement de Formes i l’Aprenentatge Automàtic. Qualsevol procés d’anàlisi de dades té com a objectiu extreure nou coneixement a partir de les dades. L’extracció de coneixement no és una tasca trivial i no es limita a la generació de models de dades o el reconeixement de patrons. L’ús de tècniques d’aprenentatge automàtic per a l’anàlisi de dades multivariades, de fet, hauria de tractar d’aconseguir un objectiu doble: la interpretabilitat i un bon rendiment. En el millor dels casos els dos aspectes d’aquest objectiu no han d’entrar en conflicte entre sí. S’ha de reconèixer la bretxa entre el modelatge de dades i l’extracció de coneixement, en el sentit que només podem extreure coneixement a partir dels models a través d’un procés d’interpretació. L’exploració de la visualització d’informació s’està convertint en una eina molt prometedora per a la interpretació dels models. Quan s’exploren les dades multivariades a través de la visualització, la gran dimensionalitat de les dades pot ser un obstacle, i moltes vegades és obligatori l’ús de tècniques de reducció de dimensionalitat. La necessitat de trobar mètodes flexibles per al modelatge de dades ha portat al desenvolupament de tècniques de reducció de dimensionalitat no lineals. L’estat de l’art d’aquests enfocaments cau moltes vegades en el domini de la modelització probabilística. Aquestes tècniques no lineals poden proporcionar una representació de les dades flexible i un model de les dades més fidel comparades amb els models lineals, però moltes vegades a costa de la interpretabilitat del model, que té un impacte en els resultats de visualització. En els mètodes d’aprenentatge de varietats amb reducció de dimensionalitat no lineals, quan un espai d’alta dimensió es projecta sobre un altre de dimensió menor, la varietat immersa obtinguda està subjecta a una distorsió geomètrica local induïda
vi per la funció no lineal. Aquest tipus de distorsió pot conduir a interpretacions errònies de l’estructura del conjunt de dades i dels patrons obtinguts. Per això, és important donar rellevància al problema de com quantificar i visualitzar aquesta distorsió en sí, amb la finalitat d’interpretar les dades d’una manera més fidel. La recerca presentada en aquesta tesi se centra en el desenvolupament de mètodes i tècniques per reintroduir de forma explícita a l’espai de visualització la distorsió local creada per la funció no lineal. Aquesta recerca se centra també en la definició de mètriques per a geometries probabilístiques per fer front al problema de la distorsió de la funció en els models de reducció de dimensionalitat no lineals. No proporcionem mètodes només per a les dades estàtiques, sinó també per a sèries temporals multivariades. La reintegració de la distorsió no lineal a l’espai de visualització dels mètodes de reducció de dimensionalitat no lineals analitzats és un objectiu en sí mateix, però aquesta anàlisi va més enllà i considera també les mètriques probabilístiques adequades a l’aprenentatge de varietats probabilístiques. Per això, estudiem el paper de les Geometries Aleatòries (distribucions de les varietats) en Aprenentatge Automàtic i anàlisi de dades en general. Es defineixen aquí els mètodes per a l’estimació de les distribucions de varietats de Riemann de suport a les dades, així com els algorismes per calcular interpolants en les distribucions de varietats. Els resultats experimentals mostren que la inferència feta segons les mètriques de les varietats Riemannianes Aleatòries dóna origen a una generació de les dades observades més fidel.
Alla mia famiglia.
Acknowledgements First I would like to acknowledge my advisor Alfredo Vellido for the great dedication that he put in the supervision of this thesis: thanks to his patient and guidance I have been able to reach this achievement. When I started this journey I had no experience as a researcher, but since the beginning he has shown interest in my opinion, making me feel as an active part of the research group. Alfredo has always been open to new ideas coming from his students and he encouraged me to think out of the schemes. He has been present as a supervisor, adapting to my rhythm, even when I needed some last-minute revision under deadlines. For these and other reasons I will tell him a big Thank You. I would also thank Neil Lawrence for his supervision and great help during my research stay in Sheffield. His priceless guidance had a great impact in my student path, stimulating me to reach challenging goals. From my first day in his group Neil made me feel like home, and his passion and enthusiasm (not only in research, but in all aspects of life) made my time in Sheffield to be an unforgettable life experience. Thank You Neil. Thanks to my co-authors for their collaboration and help. Thank to Lluís Belanche for giving value to my ideas when I was just a first year (actually, first-days) PhD student. Thanks to Ivan Olier for the patience and attention that he put in our collaboration. A special Thank You to Søren Hauberg, his enlightening ideas and his exceptional enthusiasm have been fundamental to develop our joint work; I am looking forward to sit around a table and discuss with you about extraordinary mathematical theories. I acknowledge the projects that partially founded my research: FP7 HEALTH 2013.2.4.2-1 Shockomics, TIN2012-31377 - Kappaaim, and TIN2009 13895-C02-01 - AIDTumour. Thank you to the colleagues that started with me this journey and made my experience a lot more fun. Thanks to my group-mates Sol and Albert. Thanks to all the crazy and great people who populates the floor S1 of the CS department in Barcelona, thanks for the many lunches and coffees together. Thanks to Sign’n’Forget for the acrobatic problem solving. Thanks to all the people that passed from SITraN during the last months, you made me not want to leave Sheffield (and, in fact, I did not). Grazie ad Alessandro per la pazienza (infinita) e per essere stato un amico prezioso e grazie a Luca per i tanti caffè al Blati e le lunghe chiacchierate. Grazie anche alla Maestra Anna, che mi ha fatto amare la matematica sin dal principio del mio percorso educativo, ed al Professor Bigoni, che ha accresciuto la mia passione per la scienza. The last thanks goes to Andreas, for his patience and support. Thank you for reading the whole thesis, thank you for standing my crazy moments, but mostly thanks for being by my side donating me a smile every day. Dedico questa tesi a coloro che hanno contribuito maggiormente al raggiungimento di questo obiettivo: la mia famiglia. Senza di voi non sarei mai arrivata a questo punto, mi avete sostenuta in ogni modo possibile e mi avete costantemente incoraggiata e spronata, nonché aiutata e consolata nei momenti difficili. Tutto ciò che ho ottenuto finora e che vedo per il mio futuro lo devo a voi, GRAZIE! Grazie a mamma Giulia, papà Carlo e mio fratello Raffaele, che è e sempre sarà un insostituibile punto di riferimento.
xvi List of Figures 4.7 GTM and t-student GTM diagrams for toy data of B-type: prototypes’ grid, MFs and corresponding Cartograms. . . . . . . . . . . . 55 5.1 VBGTM-TT over an Artificial dataset: dataset visualization, statemembership map and magnification factors. . . . . . . . . . . . . . 63 5.2 VBGTM-TT over the Shuttle dataset: dataset visualization, statemembership map and magnification factors. . . . . . . . . . . . . . 64 5.3 A 3−Drepresentation for the CSTP in VBGTM-TT over Artificial data and the Shuttle data. . . . . . . . . . . . . . . . . . . . . . . . . 65 6.1 Magnification Factor colormap on the GP-LVM latent space for the Shuttle data. ............................... 71 6.2 GTM latent space trained over a 3-D artificial dataset (a spiral) . . . 73 6.3 Geodesic via discretisation on GTM latent space: graph-based distance 74 6.4 GP-LVM latent space for a dataset of an artificially rotated digits dataset, together with Euclidean and Geodesic interpolants. . . . . 76 6.5 Inference over rotated digit after sampling along the Geodesic and the Euclidean distances. . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.6 Objects from COIL 100 dataset. . . . . . . . . . . . . . . . . . . . . . 77 6.7 Inference over the GP-LVM latent space for COIL example images after sampling along the geodesic and the Euclidean distances. . . . 78 6.8 Reconstruction error after inference over the GP-LVM latent space for COIL motion capture data . . . . . . . . . . . . . . . . . . . . . 78 6.9 GP-LVM latent space for the CMU motion capture data. . . . . . . . 80 6.10 GP-LVM latent space for the CMU motion capture data. . . . . . . . 81 6.11 Length of the forearm over reconstructions of data after inference over the GP-LVM latent space for CMU motion capture data . . . . 81 6.12 Examples of poses from the CMU motion capture data. . . . . . . . 83 6.13 Inference over the GP-LVM latent space for CMU motion capture data after sampling along the Euclidean distances. . . . . . . . . . . 83 6.14 Inference over the GP-LVM latent space for CMU motion capture data after sampling along the Geodesic distances. . . . . . . . . . . 83 7.1 Samples from the joint distribution over a manifold with 3-F MF visualisation................................. 88
Chapter 1 Introduction Over the last few decades, data analysis has swiftly evolved from being a task addressed mainly within the remit of multivariate statistics, to an endevour in which data heterogeneity, complexity and even sheer size, driven by computational advances, call for alternative strategies, such as those provided by pattern recognition and machine learning. These new data requirements, nowadays under the fashionable concept of Big data, come not only from business enterprises as an extension of traditional data mining, but also from scientific fields such as, for instance, biology [Marx, 2013]. The ensuing big challenge for pattern recognition and machine learning is the translation of raw data into useful knowledge that can be acted upon in practical terms. Any data analysis process, in the end, aims to extract new knowledge from data. Knowledge extraction is not a trivial task and it is not limited to the generation of data models (regardless their sophistication) or to the recognition of (possibly relevant) patterns. Those patterns and models and, in fact, any other results stemming from our analyses require interpretation to become knowledge. Consequently, the use of machine learning techniques for multivariate data (MVD) analysis should aim to achieve a dual target: interpretability and good performance. At best, both aspects of this target should not conflict with each other. A gap between data modeling and knowledge extraction must thus be acknowledged, in the sense that we can only extract knowledge from models through a process of interpretation [Vellido et al., 2012]. Models are often built with the sole goal of achieving high accuracy or precision, even though in many practical applications an optimum performance is likely to be less relevant than achieving interpretability. In this context, exploratory information visualization becomes an useful tool. When exploring MVD through visualization, high data dimensionality can be a big
2Introduction constraint, and the use of dimensionality reduction techniques becomes almost compulsory. Dimensionality reduction techniques are in fact a key tool in high-dimensional MVD analysis, and a large corpus of literature addressing this problem (mostly from the viewpoints of feature selection and feature extraction) is currently available to us, tracing back to over a century ago. The best known and most widely used linear feature extraction dimensionality reduction method is Principal Component Analysis (PCA), introduced by Pearson [1901]. In essence, PCA assumes that a low dimensional latent space, with Gaussian distributed variables, is mapped into the observed data space under a linear transformation. The reduction of dimensionality is operated by finding a few orthogonal linear combinations (principal components) of the original variables with the largest variance. The key to the resilience of this method after more than a century is, probably, its easy interpretability, as the extracted features are just linear combinations of the original ones in the data set. The need to find more flexible (and hopefully better performing) methods for MVD modeling has led to the development of non-linear techniques for dimensionality reduction, which are slowly growing in popularity [Lee and Verleysen, 2007]. A modern approach to non-linear dimensionality reduction (NLDR) of relevance to the current thesis and that involves probabilistic modeling (and, therefore, statistical machine learning) is latent variable modeling (LVM), which works by defining a subset of latent (or hidden) variables to accompany and explain the observed ones. NLDR methods can provide a flexible data representation and a more faithful model of the observed MVD than linear ones. This target is too often reached at the expense of model interpretability, which has an impact in the model visualization results. Both linear and non-linear DR methods aim, in one way or another, to preserve the structure of the observed data as much as possible in the low dimensional data representation that they generate. Unfortunately (from the point of view of interpretation) NLDR methods usually generate different levels of mapping distortion, geometrical and topological, including: manifold compression, stretching, gluing and tearing [Aupetit, 2007]. Many distortion measures (often associated with specific models and specific visualization techniques) have been proposed for different NLDR methods. In manifold learning, when a high-dimensional space is mapped onto a lowerdimensional one, the obtained embedded manifold is subject to some kind of local geometrical distortion induced by the non-linear mapping. This means that there is no guarantee that the inter-point distances in the observed data space will be uniformly reflected in the visualization space. This kind of distortion can often lead
3 to misinterpretations of the data set itself. At best, these NLDR methods can aspire to minimize the distortion of the observed data introduced in their representation, according to some objective function. But, given that it is almost impossible to completely avoid geometrical distortions while reducing dimensionality, it is important to give relevance to another aspect of the problem: how to quantify and visualize this distortion itself in order to interpret data in a more faithful way. The research reported in this thesis focuses on the development of methods and techniques for explicitly reintroducing the local distortion created by NLDR models into the low-dimensional representation of the MVD for visualization that they produce, as well as on the definition of metrics for probabilistic geometries to address this problem. For part of this research, we draw inspiration from a technique originally devised for the analysis of geographic information, namely density-equalizing maps, or Cartograms [Gastner and Newman, 2004]. These maps were originally defined as geographic maps in which the sizes of delimited regions appear distorted in proportion to underlying quantities such as their population. Cartograms were later redefined, using diffusion techniques from physics, to avoid drawbacks such as the undesired overlapping of regions, or a too strong dependence on the choice of coordinate axes. The interpretation leap in the use of Cartograms for NLDR model visualization consists on extrapolating from geographical maps to the latent visualization spaces of NLDR models (particularly manifold learning methods in this thesis), as well as on substituting geography-distorting quantities such as population density by quantities reflecting the mapping distortion introduced by the non-linear methods. We do not aim to provide methods only for static data in the thesis, but also for multivariate time series (MTS). Again, MTS visualization may become difficult to interpret when data are modelled using non-linear techniques. This is the case, for instance, when modelled using Variational Bayesian Generative Topographic Mapping Through Time (VB-GTM-TT) [Olier and Vellido, 2008a], a variational Bayesian variant of the manifold learning family defined for MTS visualization. Its interpretability will be improved through the explicit estimation of probabilities of transition between states described in the visualization space and the quantification of the non-linear mapping distortion. The reintegration of the quantified non-linear distortion into the visualization space of the analysed NLDR methods is a goal by itself, but we want to go beyond that and consider alternative adequate metrics for probabilistic manifold learning.
4Introduction To accomplish this, we study the role of Random Geometries, that is, distributions of manifolds in machine learning and data analysis in general. Methods for the estimation of distributions of data-supporting Riemannian manifolds, as well as algorithms for computing interpolants (geodesics) over distributions of manifolds are defined. In this thesis we propose methods to increase the interpretability of NLDR methods in different instances of MVD analysis using visualization. It is important to stress that this analysis could quite straightforwardly be extended not only to other variants of the methods we investigate (variants of Self-Organizing Maps (SOM) [Kohonen, 2001], Generative Topographic Mapping (GTM) [Bishop et al., 1998a; Svensén, 1998] and Gaussian Process LVM (GP-LVM) [Lawrence, 2005]), but also to other alternative NLDR visualization-oriented methods, provided a local distortion measure, or some approximation for it, could be calculated. 1.1 Summary of the main goals of the thesis The generic goals of the current thesis could be summarily listed as follows: •GG1: Exploration of the concept of local mapping distortion in non-linear dimensionality reduction methods (with specific attention paid to manifold learning techniques) from the viewpoint of the analytical quantification of such distortion. •GG2: Exploration of the Cartogram representation in bounded and partitioned visualization spaces as a tool for increasing the interpretability and usability of such multivariate data visualizations. Definition and implementation of Cartogram-based algorithms for visual representation of multivariate data for batch-SOM and GTM, based on magnification factor (MF) measurements of the mapping distortion they generate. •GG3: Explicit estimation of probabilities of transition between states described in the visualization space and quantification of the non-linear mapping distortion for VB-GTM-TT in the analysis of multivariate time series. •GG4: Definition of adequate metrics for probabilistic manifold learning as an alternative to the standard Euclidean metrics through the study of Random Geometries, including definition of methods for the estimation of distributions of data-supporting Riemannian manifolds as well as algorithms for computing interpolants over distributions of manifolds.
1.2 Structure of the document 5 1.2 Structure of the document The thesis document is structured in the following chapters: Chapter 1 The document starts with a general introduction to the field of interest, mentioning problems that will be tackled later on in this work. We then provide a summary of the notation and symbols used over the document. Chapter 2 This chapter provides the introductory technical background about probabilistic data modelling with a focus on non-linear dimensionality reduction methods for multivariate data visualization. We devote some special attention to manifold learning techniques such as the Generative Topographic Mapping (GTM) and the Gaussian Process Latent Variable Model (GP-LVM). Chapter 3 This is again a technical background chapter, which focuses on the general theme of distortion measures in non-linear dimensionality reduction methods, paying special attention to the concept of Magnification Factors. We also include some necessary basics about Riemannian geometry. Chapter 4 In this chapter, we present research results in the topic of Cartogram-based visualization of multivariate data using manifold learning models. This includes self-contained introductions to Cartogram methods and to Self-Organizing Maps (SOM) models in different variants. We show how to apply Cartograms to the representations of non-linear dimensionality reduction distortion measures in the visualization space with experiments including batch-SOM and t-GTM. Chapter 5 Our analysis departs from static i.i.d. data to address methods to improve visualization-based analysis of multivariate time series using dynamic variants of manifold learning models. This chapter includes a self-contained definition of Variational Bayesian GTM through time (VB-GTM-TT) and an experimental set. Chapter 6 It provides a study of adequate metrics for probabilistic geometries and their
6Introduction impact on model interpretability. It includes a definition of a probabilistic Riemannian metric for GP-LVM and algorithms for the calculation of geodesics distances for this model. A battery of experiments to evaluate the proposed methods is reported. Chapter 7 The final chapter summarises some conclusions of the thesis, highlighting its novelties. It also includes and lists some advanced themes and expected potential future avenues of research that we envisage beyond the advances presented in the thesis. 1.3 Refereed publications directly related to the thesis [1] A. Tosi, S. Hauberg, A. Vellido, N.D. Lawrence. Metrics for probabilistic geometries. In The 30th Conference on Uncertainty in Artificial Intelligence (UAI 2014), pp. 800–808. Quebec City, Canada. [2] A. Tosi, A. Vellido. Probabilistic Geometries as a tool for Interpretability in Dimensionality Reduction Models. The 8th WiML Workshop, Advances in Neural Information Processing Systems (NIPS 2014). Montreal, Canada. [3] A. Tosi and A. Vellido. Local metric and graph based distance for probabilistic dimensionality reduction. The Workshop on Features and Structures (FEAST 2014) International Conference on Pattern Recognition (ICPR 2014), Stockholm, Sweden. [4] A. Tosi, I. Olier, A. Vellido. Probability ridges and distortion flows: Visualizing multivariate time series using a variational Bayesian manifold learning method. In Advances in Intelligent Systems and Computing, Vol.295, pp.55-64, procs. of the 10th Workshop on Self-Organizing Maps (WSOM 2014), Mittweida, Germany. [5] A. Tosi, A. Vellido. Robust cartogram visualization of outliers in manifold learning. In Proceedings of the 21st European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN 2013), Bruges, Belgium, pp.555-560.
1.4 Other refereed publications 7 [6] A. Tosi, A. Vellido. Cartogram representation of the batch-SOM magnification factor. In Proceedings of the 20th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN 2012), Bruges, Belgium, pp.203-208. Description of contributions in relation to the refereed publications In Chapter 3 we introduce the Cartogram-based method which has been presented in [6] and [5] for the batch-SOM algorithm and the t-GTM respectively; in Sec. 4.3 and Sec. 4.4 we illustrate the experimental results as reported in the two papers. The analysis on multivariate time series reported in [4] is presented in Chapter 5, together with the experimental results; the design and implementation of the algorithms has been done using software tools for VB-GTM-TT [Olier and Vellido, 2008b] provided by Dr. Iván Olier. The idea of probabilistic geometries presented in Chapter 6 is the result of a joint project carried out at the Machine Learning group of the University of Sheffield, with Professor Neil Lawrence. The experimental results of [1] are reported in Sec. 6.4. The tools used to design the experiments and train the models are built on the GPLVM [Lawrence, 2005] software provided by the Machine Learning group of the University of Sheffield1; the tools used to compute manifold structures and geodesics via ODE’s solutions, as described in 6.3.2, use the software provided by Dr. Søren Hauberg and previously used in [Hauberg et al., 2012]. The description of the geodesic computation via discrete graphs presented in 6.3.1 refers to [3]. Finally, the work presented in [2] provides an overall summary of the topics of this thesis, focusing on the problem of interpretability in dimensionality reduction introduced in Chapter 2 and Chapter 3 and presenting the advances detailed in Chapter 6. 1.4 Other refereed publications It is worth mentioning other publications that, although not included in the thesis, have inspired its development while exploring new and interesting topics of research: [7] L. A. Belanche, A. Tosi. Averaging of Kernel Functions. Neurocomputing 112 (2013), pp.19-25. 1Software can be downloaded here: https://github.com/SheffieldML/, both in the Matlab and in the Python version
8Introduction [8] L. A. Belanche, A. Tosi. Averaging of kernel functions. In Proceedings of the 20th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN 2012), Bruges, Belgium, pp.363368. 1.5 Symbols and notation •Notation: The matrix Y∈RN×prepresents the observed data space, where each row corresponds to an observed data point and each column to a dimension. We denote with Y:,j the columns of the data matrix, with yithe rows of the data matrix and with yi,j a single scalar element. Y= [Y:,1Y:,2. . . Y:,p] = data features z }| { y1,1y1,2· · · y1,p y2,1y2,2· · · y2,p · · yN,1yN,2· · · yN,p data points Similarly, we denote with xithe rows of the matrix X. In this document we use the symbol xto represent a latent vector of dimension q, and the symbol yto represent an observed vector (data point) of dimension p. We always consider q < p. Given a differentiable function f:Rq−→ Rp x7→ y=f(x) we call Jacobian the p×qmatrix Jcontaining all the partial derivatives J= ∂y(1) ∂x(1) . . . ∂y(1) ∂x(q) . . .. . . ∂y(p) ∂x(1) . . . ∂y(p) ∂x(q) (1.1)
1.5 Symbols and notation 9 •Symbols: Rthe set of real numbers Iidentity matrix Φ(·)a vector of function values, the mth element corresponds to ϕm(·) f(X)a vector of function values, the ith element corresponds to f(xi). Kf,fcovariance matrix whose elements are given by k(xi,xj) x(i)the ith component of the vector x ∂ ∂x(i)the partial derivative with respect to x(i) Jthe Jacobian of a function ▽2the Laplacian operator of a function ∼distributed according to the following probability distribution GP(·,·)Gaussian Process N(·,·)Gaussian Distribution Γ(·,·)Gamma Distribution E[x]expectation of the random variable x
16 Probabilistic Modelling centres µmare constrained to lay on an intrinsically low-dimensional space. In the standard case, the basis distribution is chosen to be a set of Gaussian radial basis functions with the same lengthscale γ: ϕm(x) = exp −γ 2∥x−µm∥2,(2.6) but other distributions can be considered for different types of data; one example will be provided later on in § 4.4.1, in which a mixture if Student-t distributions is used to define a GTM variant that behaves robustly in the presence of atypical data or outliers. The centres of the mixture of distributions can be interpreted as data prototypes or cluster centroids that can be further agglomerated in a full-blown clustering procedure. In this manner, GTM combines the functionalities of Self-Organising Maps (c.f. section §4.3) and mixture models by providing both data visualisation over the latent space and data clustering [Olier and Vellido, 2008c]. Provided a prior distribution over the latent space, this model leads, in a similar way to probabilistic PCA (PPCA) [Tipping and Bishop, 1999], to a Gaussian conditional distribution of the data p(yi|x,W, β) = N y M X m=i w⊤ jϕm(xi), β−1I!(2.7) =β 2πD/2 exp −β 2 p X j=1 yi,j − M X m=i w⊤ jϕm(xi)!2 .(2.8) Model likelihood and expectation maximization The GTM algorithm aims to find the probability p(y|W, β)of a data point given the adaptive weight parameters Wand the noise variance β. To do so, the latent vectors xare integrated out of the model. To make the computation analytically tractable, the prior distribution p(x)is defined by a set of Kequally weighted delta functions p(x) = 1 K K X k=1 δ(x−xk).(2.9) The Kcentres xkare distributed in a predefined regular lattice, which is taken to be squared in the standard case (other choices are allowed, for example hexagonal grids). This discrete choice of the prior distribution simplifies the integration and,
2.2 Generative Topographic Mapping 17 as a result, the data distribution becomes p(y|W, β) = Zp(y|x,W, β)p(x)dx =1 K K X k=1 p(y|xk,W, β),(2.10) and assuming the data points i.i.d., we obtain the following final expression of the model likelihood: L= N Y n=1 p(yn|W, β).(2.11) To estimate the parameters Wand β, we can use a maximum likelihood approach, which is equivalent to consider the maximum of the log-likelihood ℓ=log(L) = N X n=1 log 1 K K X k=1 p(yn|xk,W, β)!.(2.12) The optimization of the adaptive parameters can be achieved by any standard non-linear optimization technique (see e.g. [Press et al., 1988]) but, since we are working with a mixture of Gaussians, the most common choice is to use the expectationmaximization (EM) algorithm [dem; Bishop, 1995]. Given the initial values for W and β, the E-step for the standard GTM formulation is the same as for the general Gaussian Mixture model, where the the conditional probability of each latent point given each observed data point is computed using Bayes’ theorem. The probabilities are usually referred to as the responsibilities rkn rkn ≡p(xk|yn,W, β, ) = p(yn|xk,W, β)p(xk) PK k′=1 p(yn|xk′,W, β)p(xk′).(2.13) Notice that, for the choice of prior distribution made in Eq: (2.9), the effect of the term p(xk)is cancelled. Considering now the choice of the GTM mapping made in Eq. (2.5), we obtain that the M-step of the EM algorithm reduces to the solution of a set of linear equations. For more details about the EM algorithm and the update equations for the standard GTM see [Bishop et al., 1998a, § 2.2].
18 Probabilistic Modelling Fig. 2.3 Visualization of modes projections (black dots, left diagram) and means projections (black dots, right diagram) on the 2-D GTM latent space for the 3-D artificial dataset in Fig. 2.2. Data Visualization Since we are especially interested in data visualization, we can use the conditional probability defined by Eq: (2.13) to obtain both a posterior mode projection of yn kmode n=argmax k rk,n (2.14) (which implies assigning each observed data point to that latent point with the highest responsibility for its generation), or a posterior mean projection xmean n= K X k=1 rknxk(2.15) (locating the observed data point at a location in latent space that results from a responsibility-weighted combination of all latent point locations). We can now visualize (see Fig. 2.3) the observed data points over the low-dimensional latent space using the posterior mean projection, which also provides an assignment of each data point to a representative cluster. For a visualisation purpose, the typical setting is a latent dimension q= 2 or 3. 2.2.2 Time-series analysis with GTM-TT When the observed data space is known to be in the form of a time series, the timedependent nature of the observations makes the assumption of i.i.d. inappropriate. In order to make the GTM model suitable to the analysis of temporal data, we consider here an extension within the framework of hidden Markov models (HMMs)
2.2 Generative Topographic Mapping 19 [Rabiner, 1989]. This model is known as GTM Through Time (GTM-TT) [Bishop et al., 1997b] and was proposed almost in parallel to the standard GTM. The GTM-TT can be interpreted as a standard GTM model in which the latent points are considered as hidden states Z={zt}t=1,...,T for every time step t. Similarly to HMMs, the states are connected by a transition probability aij =p(zj|zi), which represents the probability of making a transition to the state jfrom the current state i. We denote with Athe matrix which describes the transition states A={aij}:aij =p(zt=xj|zt−1=xi), i, j = 1, . . . , K (2.16) where Kis the number of allowed hidden states xk(which is the number of vector prototypes). Given the initial state probabilities on each of the latent points at the first time step t= 1 π={πk}:πk=p(z1=xk),(2.17) then the parameters governing the GTM-TT model are Θ= (π,A,W, β).(2.18) Notice that the parameters Wand β, together with the transition probabilities A, are common to all time steps of the GTM algorithm, so that the number of adaptive parameters in the model is independent of the length of the time series. The adaptive parameters of the model can now be computed (in a similar way to GTM) using a maximum likelihood approach, via EM algorithm. In the context of HMMs, this is generally known as the Baum-Welch algorithm. Given an observed p-variate time series Y={yi}i=1,...,p, the complete data log-likelihood is given by ℓ= K X k=1 vk,1log πk+ N X n=2 K X i=1 K X k=1 vi,n−1vk,n log aik +pN 2log β 2π−β 2 N X n=1 K X k=1 vk,n ∥yn−f(xk,W, β)∥2,(2.19) where the binary vector vnsuch that its component vk,n returns 1 if znis in state k, and zero otherwise (this indicators are suitable to simplify the expression of the
20 Probabilistic Modelling likelihood in order to perform the EM steps). Moreover, vnsatisfies PK k=1 vk,n = 1. A detailed description of the updating equations of the EM algorithm for GTM-TT can be found in [Bishop et al., 1997b, § 4]. Visualization of Time Series In the same fashion as GTM, the GTM-TT allow data visualization simultaneously to data clustering. In this way we have a low-dimensional (usually q= 2) latent space where the multivariate time series is represented by the means of the posterior-mode projection, defined as k(mode) n=argmax k rk,n (2.20) where rk,n are the responsibilities probabilities rk,n ≡p(zn=xk|yn,Θ)(2.21) This model, even if useful for MTS clustering and visualization, does not involve any regularization process. 2.3 Gaussian Processes Latent Variable Models In this section we present the GP-LVM model, which is a Gaussian Process based dimensionality reduction model. To do so, in this section we firstly review the theory of Gaussian Processes (GPs) in section, including GPs for regression and some intuition about covariance functions. After this, we describe the details of the GP-LVM model. 2.3.1 Introduction to Gaussian Processes A GP is used to describe distributions over functions and it is defined as a collection of random variables, any finite number of which have a joint Gaussian distribution [Rasmussen and Williams, 2006]. Let the vector x∈Rqand the function f:Rq→R. A GP is a stochastic process determined by its mean function µ(x)and its covariance function k(x,x′), and it is denoted as f(x)∼ GP(µ(x), k(x,x′)),(2.22)
2.3 Gaussian Processes Latent Variable Models 21 The function ftakes values over continuum on the input space (infinite input values). But, in practical applications, we only consider a finite set of instantiations of the function, since (computationally) we can only have access to a finite number of input vectors. If we consider a collection of inputs X={xn}n=1,...,N , we can generate a random vector of function values f=f(X) = (f(x1), . . . , f(xN)) ∈RN(2.23) which is Gaussian distributed with covariance matrix Kgiven by the gram matrix of the covariance function k, denoted with Kf,f Kf,f= k(x1,x1)k(x1,x2)· · · k(x1,xN) k(x2,x1)k(x2,x2)· · · k(x2,xN) · · · · · · · · · · · · k(xN,x1)k(xN,x2)· · · k(xN,xN) (2.24) In this document we equivalently refer to the covariance function kas the kernel function or simply the kernel. The kernel kdefines the correlation between two inputs and can be interpreted as a similarity or as a distance into the functional space of positive semidefinite kernels k. There is, in fact, a relation between GP prediction and the regularization theory in reproducing kernel Hilbert spaces (RKHSs), as described in detail by Rasmussen and Williams [2006, § 6]. Without loss of generality, the mean function is typically chosen to be equal to zero. An intuitive understanding of the distribution of f(X)can be gained if we think of it as a marginal distribution. In fact, the input space Xcan be divided in two sets of input vectors, the observed ones Xand the unobserved (potentially infinite) ones: due to the properties of joint Gaussian distributions, we can integrate over the unobserved variables. A sample of the function fcan be obtained by sampling from its distribution following the standard procedure for Gaussians. (For more technical details about Gaussian distributions and mathematical identities check Appendix A). Gaussian Processes for regression Let’s consider a regression problem in which a set of observations X={xn}n=1,...,N is mapped into a set of outputs Y={yn}n=1,...,N . We want to find the function f that maps the inputs into the outputs. Using Bayesian statistics, we can perform in-
22 Probabilistic Modelling ference over the function values by combining two pieces of information: our prior belief over the properties of f, encoded in the prior distribution, and the information given by the input data, encoded in the likelihood distribution. We assume that Y constitutes a noise-corrupted version of f(X), according to Eq. 2.1. That is, we assume that ynit is obtained by the correspondig f(xn)with the addition of a Gaussian noise of β−1variance. Then we have p(f|X,Y) = prior z }| { p(f|X) likelihood z }| { p(Y|X,f) p(Y|X) | {z } marginal likelihood (2.25) Assuming i.i.d. inputs Xand a likelihood following a Gaussian distribution, if we use a GP prior over the mapping fwe have p(f|X,Y) = N(f|0,Kf,f) N Y n=1 N(yn|f, βI).(2.26) The computation of the marginal likelihood has been solved analytically using the matrix determinant lemma and the Woodbury identity. Prediction of an unobserved function value f∗=f(x∗)computed in a test point x∗is obtained considering the joint distribution "f f∗#∼ N 0,"Kf,fKf∗,f K⊤ f∗,fKf∗,f∗#!,(2.27) where the elements of the covariance matrix Kf,fare given by 2.24 and Kf∗,fand Kf∗,f∗are Kf∗,f= k(x∗,x1) k(x∗,x2) · · k(x∗,xN) (2.28) Kf∗,f∗=k(x∗,x∗) and, applying the properties of the conditional Gaussians (see Appendix A), it follows p(f∗|f,X) = N(K⊤ f∗,fK−1 f,ff |{z } mean ,Kf∗,f∗−K⊤ f∗,fK−1 f,fKf∗,f | {z } covariance ).(2.29)
2.3 Gaussian Processes Latent Variable Models 23 In applications to real datasets we are typically interested in prediction of noisy observations (as seen in Eq. 2.1). To do so, we can assume an additive i.i.d. Gaussian noise ϵ∼ N(0, β−1)and we use a GP covariance ˜ K=K+β−1I. More details will be provided in section §2.3.2 (in particularly Eq. 2.34) for GP-LVM. Linear White noise −2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2 −2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2 Fig. 2.4 Examples of different GP prior distributions over faccording to different covariance functions: linear and noisy. Hear each function in different colors represents a different sample of f. The mean function is equal to zero, due to our choice of m(x) = 0. The variance of fis here represented by the pink area (corresponding to the 95% confidence interval). Covariance functions In general, a GP prior is fully defined by its covariance. We present here the basic concepts about covariance functions needed for the understanding of the next chapters, and we refer to [Rasmussen and Williams, 2006, § 4] for a more detailed analysis. Let’s consider an input domain X ∈ Rq. A covariance function k(x,x′)(also known as kernel function) is a positive definite function which, intuitively, it is used to describe the similarity between two inputs. A wide range of covariance functions can be used in GP models, and a correct choice might be crucial for specific applications. Each covariance function encodes, in a different way, the properties of the function we wish to learn. We can see in Fig. 2.4 and Fig.2.5 some examples of GP samples from different covariance functions: linear, white noise and periodic. A widely used covariance function is the exponentiated quadratic (also known as
24 Probabilistic Modelling Periodic Covariance −2 −1 0 1 2 10 20 30 40 50 5 10 15 20 25 30 35 40 45 50 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fig. 2.5 Example of a Gaussian distribution of functions according to period kernel function. The corresponding covariance matrix is displayed on the right. We use 50 inputs equally distributed over the horizontal axis, with values taken between -1 and 2. High values of the covariance correspond to high correlation: notice that input points are correlated with a periodic structure (for example, the input x1is highly correlated with the inputs x17, x33). squared exponential or RBF kernel): k(x,x′) = αexp −ω 2∥x−x′∥2 2.(2.30) The popularity of the exponentiated quadratic covariance function is due to the fact that it is capable of smoothly modelling different kind of functions only by varying the value of the lengthscale 1/ω. We can, in fact, find sensible initialization of the model by considering that the lenghtscale is proportional to the number of points where the function is crossing the zero axis. In the left column of Fig. 2.6 we see different examples of GP priors over f. In the central and right columns of Fig. 2.6 we see the posterior distribution of f, after the observation of some data: notice that, far from the observed data, the distribution of the function tends to revert to the prior. The variance of the random variable is here represented by the pink area (corresponding to the 95% confidence interval). 2.3.2 Dimensionality reduction with GPs: the GP-LVM Considering the mapping given by the noise model in Eq. (2.1), GPs have been used in probabilistic non linear dimensionality reduction to define a prior distribution over the mapping f, leading to the formulation of the Gaussian process latent variable model (GP-LVM).
2.3 Gaussian Processes Latent Variable Models 25 −2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2 −2 −1 0 1 2 −2 −1 0 1 2 −2 −1 0 1 2−2 −1 0 1 2 −2 −1 0 1 2 −2 −1 0 1 2 −2 −1 0 1 2−2 −1 0 1 2−2 −1 0 1 2 −2 −1 0 1 2 −2 −1 0 1 2 Fig. 2.6 This diagram shows, in each row, the effect of changing the lenghtscale in the exponentiated quadratic covariance function (from top to bottom: ω= 0.2,1,5). In the first column we can see how the GP prior distribution is affected by this change. The second and third columns show some samples from the posterior after seeing 2 and 9 data points. Each function in different colors represents a different sample of f. The mean function is in black. The variance of fis here represented by the pink area (corresponding to the 95% confidence interval). Lawrence [2005] introduces GP-LVM pointing out its dualistic relation with probabilistic principal component analysis (PPCA). The expression of the marginal likelihood of the PPCA model, as defined in Eq. 2.4, is a tractable integral which leads to the following solution p(Y|W, β) = N Y n=1 N(yn|0,WW⊤+β−1I).(2.31) We now want to find the values of the parameters Wwhich give the maximum value for the marginal likelihood of the data. As suggested in [Roweis, 1997; Tipping and Bishop, 1999], the solution of this problem can be treated as an eigenvalue problem, resulting in a fast computation. In the standard dimensionality reduction approach described before, the likeli-
32 Some Tools for Improving Interpretability dealing with noise and observations with complex structure and non linearities, the use of the tools of Riemannian geometry are becoming more popular among the machine learning community. In fact, we show that we can interpret an embedded Riemannian manifold as the underlying support of the data distribution. For every point on the manifold, all the geometrical properties are specified by a local positive definite matrix, called the metric tensor. Once this tensor is known, it is then possible to compute interesting objects as length minimizing curves (known as geodesics) and magnification factors (introduced in chapter § 3.3). 3.1 Dimensionality reduction and distortion measures The use of dimensionality reduction techniques is an essential tool when dealing with visualization oriented applications involving high-dimensional data. Examples of popular models are Principal Component Analysis (PCA) and multidimensional scaling (MDS). PCA, independently introduced by Pearson [1901] and Hotelling [1933], provides a low dimensional linear representation of the data set which captures the most variation in the observed high-dimensional variables, while MDS [Mardia et al., 1979] aims to preserve distances and proximities between pairs of observations by the definition of a similarity matrix. These methods are easy to interpret for practical purposes, but they also suffer from some limitations. For example, PCA is very useful for reducing redundancy of features in the original data set, but is restricted to the assumption of linearity, and MDS relies on the preservation of local distances, which is not always the best approach. The need to find less constrained (more flexible) methods for multivariate data modelling has led many researchers to explore and define non-linear techniques of dimensionality reduction (NLDR), which are becoming increasingly popular [Lee and Verleysen, 2007]. The most interesting contributions to this area range from spectral-based methods to manifold learning techniques. Some examples of spectral approaches include kernel PCA [Schölkopf et al., 1997], local linear embedding (LLE) [Roweis and Saul, 2000], isometric feature mapping (ISOMAP) [Tenenbaum et al., 2000], Laplacian eigenmaps (LE) [Belkin and Niyogi, 2003], or maximum variance unfolding (MVU) [Weinberger and Saul, 2006]. Advantages of spectral approaches include the reach of a global optimum and a smooth mapping from the data space to the low-dimensional space. On the other
3.1 Dimensionality reduction and distortion measures 33 side, when considering a probabilistic framework, the advantages typically include the explicit access to a (smooth) mapping from the low dimensional space to the high dimensional observed space. Other relevant advantages include marginalization of missing data, model selection through Bayesian inference and integration with other models (such as mixture models or temporal models). Specific attention will be paid, in the following sections, to the analysis of probabilistic generative models of the latent variable models family. In particular, we will focus part of our research on models which can be expressed by the mapping given in Eq. (2.1). These models include the aforementioned models in § 2.1. In general, NLDR techniques attempt to minimize the unavoidable distortion that they introduce in the mapping of the high-dimensional data from the observed space onto lower-dimensional spaces. For a more faithful interpretation of models, a large number of distortion measures have been proposed and adapted to visualization techniques for different NLDR methods. While reducing dimensionality, NLDR generate different levels of local mapping distortion, that lead to a loss of information that we aim to recover, to some extent, into the visualization space. Stretching or compressing a space affects the preservation of pairwise distances between points. An example of such non-linear mapping is displayed in Fig. 3.1. An interesting contribution comes from research by Aupetit [2007], who identifies different types of distortion, classified as geometrical and topological (including: manifold compression, stretching, gluing and tearing) and proposes the use of Voronoi diagrams and colour scales to visualize manifold-based measures such as point-based, segment-based and triangle-based measures. The non-linearity of dimensionality reduction methods such as the Generative Topographic Mapping (GTM) and the Gaussian Process Latent Variable Model (GPLVM) entails the existence of local distortion in the mapping of the data from the observed space onto the visualization space. This fact limits the direct interpretation of the visual data representation and there have been efforts to provide visual solutions to this limitation by defining and visualizing DR quality measures that, embedded in the method, can be associated to each data point, using colouring procedures for the data-corresponding cells in the Voronoi tessellation of the projection space. A widespread used method for Self-Organising Maps (SOM) (which is not a true latent or generative model), called Unified distance Matrix (U-matrix), proposed in [Ultsch, 2003], allows the visualization of the pairwise distances between corresponding points in the original data space on the pseudo-latent space of the model.
34 Some Tools for Improving Interpretability The values of these distances can be represented with a color map accompanying the SOM topographic grid. Another interesting approach, proposed in [Bishop et al., 1997a] for GTM (and extended to SOM), is the calculation of a magnification measure in a continuous way over the representation map. We delve into this technique in §3.3. We show in the following section that, while dealing with generative LVMs where the mapping is continuous and differentiable, the expression of magnification factors (MFs) can be analytically derived. To do so, we will will use the tools of Riemannian geometry in order to formalize this approach. Fig. 3.1 Mapping between a 2-D plane and a surface embedded in a 3-D space. A straight line is subject to distortion under the non-linear projection. 3.2 Concepts of Riemannian Geometry We study latent variable models as embeddings of uncertain surfaces (or manifolds) into the observation space. From a machine learning point of view, we can interpret this embedded manifold as the underlying support of the data distribution. To this end, we review the basic ideas of differential geometry, which study surfaces through local linear models. Gauss’ study [1827] of curved surfaces are among the first examples of (deterministic) latent variable models. He noted that a 2-dimensional surface Membedded in a3-dimensional Euclidean space is well-described through a collection of mappings xαfor each point ξ∈ S xα:U∈R2→V∩ S ∈ R3(3.1) α7→ ξ(3.2) Intuitively, the regular surface is defined by a collection of open sets of R2(small
3.2 Concepts of Riemannian Geometry 35 pieces of planes) stick together in a way that the transition from a set to an other is made in a smooth way (with no sharp points, no self-intersections and no edges). Historically, Gauss considered the case of two-dimensional surfaces embedded in R3, while the extension to higher dimensional manifolds is due to his student Bernhard Riemann [1854]. The pair (xα, Uα)is called parametrization (or system of coordinates) on the manifold. A differential manifold is characterized by mappings xαthat are smoothly varying (i.e., varying in a differentiable way) between the open sets Uα. Given a smooth manifold, we can take advantage of its differential structure in order to make computations wit its elements. We can, in fact, define a smoothlyvarying inner product on the tangent space TM: this is a Riemannian metric. Definition (Riemannian Metric): A Riemannian metric1Gxon the differential manifold Mis symmetric and positive definite matrix associated with a smoothly varying inner product ⟨·,·⟩x ⟨a,b⟩x=a⊤Gxb(3.3) on the tangent space TxM, for each point x∈ M and a,b∈TxM. The matrix G is called the metric tensor. The pair (M,⟨·,·⟩x)is called Riemannian manifold. Example 3.2.1: The most common and fundamental example of Riemannian manifold is the Euclidean space equipped with the canonical inner product (R, can). Example 3.2.2: Suppose that we have an embedding f:M −→ N, and that (N,⟨·,·⟩y) is a Riemannian manifold. We can construct a Riemannian metric on Mby pulling back ⟨·,·⟩yto ⟨·,·⟩x=f∗⟨·,·⟩yon M. In other words we have: a⊤Gxb=⟨a,b⟩x= ⟨∂ ∂x f(a),∂ ∂x f(b)⟩ywhere G=J⊤J(3.4) Remark: The Riemannian metric needs not be restricted to G=J⊤Jand can be any smoothly changing symmetric positive definite matrix [do Carmo, 1992]. In this research we restrict ourselves to the more simple definition of Eq. 3.4 as it suffices for our purposes, but note that the more general approach has been used in machine learning, e.g. in metric learning [Hauberg et al., 2012] and information geometry [Amari and Nagaoka, 2000]. 1For simplicity we use xinstead of xα(α)to denote a point on the manifold. The subscript xin Gxis omitted when there is no ambiguity.
36 Some Tools for Improving Interpretability The Riemannian metric encodes the geometrical properties of the manifolds and can be employed to compute useful quantities, such as curvatures, arc lengths and angles. In particular, we are interested in computing distances. Definition (Geodesic Curve): Given two points x1,x2∈ M, a geodesic is a length-minimising curve connecting the points γg= arg min γLength(γ), γ(0) = x1, γ(1) = x2,(3.5) where the length of a general curve γ: [0,1] →Rqis computed as Length (γ) = Z1 0q⟨γ′(t), γ′(t)⟩γ(t).dt (3.6) It can be shown [do Carmo, 1992] that geodesics satisfy the following second order ordinary differential equation γ′′ =−1 2G−1∂vec G ∂γ ⊤ (γ′⊗γ′),(3.7) where vec Gstacks the columns of Gand ⊗denotes the Kronecker product. The Picard-Lindelöf theorem [Tenenbaum and Pollard, 1963] then implies that geodesics exist and are locally unique given a starting point and an initial velocity. 3.3 Magnification Factors The metric tensor defines the local geometrical properties of the considered generative LVM model and it can be used as a tool to data exploration. One way to visualise the tensor metric is through the differential volume of the high dimensional parallelepiped spanned by the mapping; this, for a latent dimension q= 2 is known as MF and it was introduced by Bishop et al. [1997a] for GTM (and standard SOM). Its explicit formulation, using the notation given in Eq. 3.4, is given by MF = pdet (J⊤J).(3.8) The MF points out the non-linear distortion generated by the projection of the observed data onto the representation map. Being computed over continuum on the input space, the MF values can be visualised as a colormap over the 2-D display of the latent space of the chosen model. We
3.3 Magnification Factors 37 display in Fig. 3.2 a colormap representing MF values compute over a fine regular grid over the 2-D latent space of a GP-LVM model. The explicit formulation of the MF for the GP-LVM model is detailed later on in this document, in section §6.2. MF 50 100 150 200 250 300 Fig. 3.2 Magnification Factor colormap on the GP-LVM latent space for a jogging motion from the CMU motion capture database. Other examples of this type of this visualisation are provided throughout this document in the sections describing the experimental results. See Fig. 4.5 (top left), Fig. 4.6 (center), Fig. 4.7 (center), Fig. 5.1 (bottom right), Fig. 5.2 (bottom right), Fig. 6.1, Fig. 6.3 and Fig. 6.9. We also introduce a 3-D plot of the MF colormap, in which the vertical axes represents the values of the MF. In Chapter 6 we define a distribution over the local metric tensor for generative latent variable models and in this context the 3-D display is used to visualise the variation of MF values between different samples from the metric, as in Fig. 7.1. Another interesting example of visualization is provided in [Svensén, 1998, § 4.5], where the author visualizes the different levels of magnitude and directions of stretch of the MF over the GTM latent space. The concept of MF has its origin in the field of computational neuroscience, where it evaluates the mapping distortion between the spatial density of biological
38 Some Tools for Improving Interpretability sensors and the two-dimensional spatial density of the corresponding topographic maps in the visual and somatosensory areas of the cortex. More specifically, the cortical MF would indicate the linear distance along the primary visual cortex concerned with each degree of visual field [Pointer, 1986], although controversy remains on whether the cortical magnification of the central visual field reflects its selective amplification, or merely reflects the ganglion cell density of the retina [Wässle et al., 1990]. As expressed in the context of vector quantization models [Hammer et al., 2007], local magnification is the result of a specific connection of the density of model prototypes and stimuli. One of the most interesting facts is that the distortion caused by the non-linear mapping can be explicitly quantified as a MF over the latent space used for data visualization. For this property, the concept of magnification has recently been applied to manifold learning methods for NLDR, in order to visualize the distortion due to the embedding of a manifold in a high-dimensional space [Tosi and Vellido, 2012; Vellido et al., 2013; Tosi and Vellido, 2013; Gianniotis, 2013; Tosi et al., 2014b,a]. More details about novel visualization techniques involving MF are presented in the following chapter.
Chapter 4 Advances in Mapping Distortion Visualization for Non-linear Dimensionality Reduction Using Cartograms One of the main advantages of non-linear dimensionality reduction methods, as mentioned in previous chapters, is their flexibility in the process of multivariate data modelling. Unfortunately, this flexibility is also accompanied by limitations, such as the difficult interpretation of the data (visual) representations they generate. Even latent variable manifold learning models, which represent multivariate data in low-dimensional representation spaces, can be difficult to make sense of, due to the fact that their coordinates in latent space are complex non-linear transformations of the observed ones, so that heterogeneous levels of local distortion are generated. These locally varying distortions and the loss of straightforward meaning in the low dimensional variables make tools to assist their interpretation an almost compulsory requirement. Linear dimensionality reduction methods, on the other hand, are less flexible and their data representation can be less faithful as a result, but they compensate for this with the straightforward interpretation of their representation coordinates. This comes a long way to explain the popularity of linear dimensionality reduction techniques and the difficulties for non-linear ones to become mainstream. Given that some of the modelling techniques considered in this thesis are of non-linear nature, we are faced with an interpretability challenge. In the current chapter, we respond to this challenge using a technique for data visualization in non-
40 Advances in mapping distortion visualization for NLDR using Cartograms linear latent models that is inspired in geographical representation methods. This technique is suited to both static data and multivariate time series representation. The technique in question, known as Cartogram, draws inspiration from both geographic representations and physics. We introduce novel variants of Cartogrambased visualization for NLDR techniques. We illustrate the proposed method providing a self-contained description of the algorithms for the visualization of the BatchSOM algorithm in § 4.3 and the robust tGTM algorithm in § 4.4. 4.1 Cartograms In the area of thematic representations for geography, there is a particular type of mapping known as Cartogram. In this mapping, specific areas, often delimited by political borders, are locally distorted (both stretched or compressed) to convey the information on locally-varying underlying quantities of interest such as population density or socio-economic data. Examples of applications of Cartogram techniques include the visualization of census data, disease incidence, birth rate, and annual income. Fig. 4.1 Cartogram representation of the world map according to the population density: countries with a higher population density are represented with an area that is bigger than the geographical one. In this way it is easy to visually notice the contrast between countries with high population densities (such as India or Japan) and countries with lower population density (such as Australia or Canada). Source of image: http://www.worldmapper.org/ display.php?selected=2. In the last decades, the development of feasible computer-based Cartograms has been a challenging undertaking and several methods have been proposed to accomplish it. An extended review of the recent history of this subject is given by
4.1 Cartograms 41 Tobler [2004], where the author completes and unifies pioneering work started in the 1960’s. Nowadays, the use of cartograms for the visual representation of socioeconomic data in geographical maps has become popular thanks to public resources such as Worldmapper1(An example is provided in Fig. 4.1). From a mathematical point of view, the geometrical distortion of Cartograms takes the form of a continuous transformation from R2to itself. Fig. 4.2 represents the distortion of a square patch on the original plane: here a vector x= (x(1), x(2)) is mapped onto a vector x′= (x′(1), x′(2))according to a transformation Tin such a way that the Jacobian of the transformation is proportional to an underlying distortion variable d, which describes the deformation of the patch quantitatively : T:R2→R2 x→ T (x) = x′ ∂Ti ∂x(j)i,j ∝d. (4.1) Fig. 4.2 A continuous transformation is applied to a 2-D square, which is mapped into a distorted patch (dashed line). The Cartogram projection is not determined uniquely by the condition 4.1, since the problem has two degrees of freedom. To fix the projection we need one more constraint and, to do so, we can operate in many different ways. Since there is not a single choice for the constraints of the projection T, some choices may result in loss of connectivity between the fragment borders or overlapping of neighbours areas. A method for building Cartograms, based on the physics’ principle of linear diffusion processes, has been proposed by Gastner and Newman [2004]. The principle of diffusion applies the theory of parabolic partial differential equations to the problem of diffusion of a large number of particles, knowing their density as a function of the position and the time. In a natural system, the particles will flow, over time, from areas of higher density to areas of lower density, resulting in a final state where the overall density is homogeneous. In order to apply this model to the problem of Cartograms, the distorting variable dis interpreted as a diffusion function ρ(x, t) depending on the position x∈R2and time t∈[0,∞). The diffusion function is 1www.worldmapper.org
48 Advances in mapping distortion visualization for NLDR using Cartograms three clusters which, in fact, are far from each other, as evidenced by the direct data visualization. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Fig. 4.5 Top row: left) Map of MF values together with a colorbar on the right-hand side of the map; right) corresponding Cartogram. Bottom row: left) U-matrix map; right) corresponding Cartogram. We compute the MF value for each node xkusing Eq. 4.8. In the experiments, the batch-SOM maps are transformed into a Cartogram by using the rectangular grid, defined by the nodes xk, as map internal boundaries (effectively, defining a centroidal Voronoi tesselation [Du et al., 1999]). As described in section § 4.2, we assume that the level of distortion dexternal in the space beyond this rectangle is uniform and equal to the mean distortion over the complete map. A similar procedure is applied to compute the Cartogram using the unified distance matrix as the underling distortion measure. The U-matrix [Ultsch, 2003] allows the visualization of the pairwise distances between corresponding points in the original data space on the low-dimensional visualization space of the SOM topographic
4.4 Cartogram representations for GTM 49 grid. Results are displayed in Fig. 4.5. The overlaid grid of reference vectors seen in Fig. 4.3 explains the fact that many reference vectors are squashed in data-dense regions whereas only a few are stretched to cover the empty space in-between. This varying distortion is nicely reflected by both the MF and the U-matrix, on the left column of Fig. 4.5. The batch-SOM map and the distortion measures finally come together in the Cartograms (MF: top-right, and U-matrix: bottom-right of Fig. 4.5. The empty spaces between clusters are now fairly stretched, providing a clear view of the separation. Interestingly, part of the data reside in stretched areas: These are the ones further from the cluster centres. This effect should warn us against a too straightforward interpretation of high-distortion areas as completely empty spaces. The visualization of the MF on the batch-SOM map may inform us of the existence of data clusters and the sparsely populated spaces that separate them, as they undergo different levels of distortion: low in dense areas, while high in empty ones. In this task, it is a principled alternative to the widely used U-Matrix [Ultsch, 2003]. This direct visualization is not always intuitive. Instead, the Cartogram-based representation of the batch-SOM map retains its simplicity while visually factoring out the non-linear distortion as measured by the MF. 4.4 Cartogram representations for GTM The same idea proposed above for batch-SOM algorithm can be applied to the Generative Topographic Mapping (GTM), cfr section 2.2. In the following section, we propose here to explore some artificial data sets with simple statistical properties, in order to assess the properties of the method in a controlled setting. We will work with 3-D data divided into three neatly defined clusters, to which atypical cases or outliers will be added. The dimensionality of data is chosen, again, to allow the direct visualization of prototypes embedded in the observed data space. Following the line of the experiment described in § 4.3, we explicitly calculated the MF for the t-GTM variant of the standard GTM and we introduce a Cartogram visualization of the latent space. The standard GTM algorithm was implemented in Matlab®, using the drtoolbox3. 3http://homepage.tudelft.nl/19j49/Matlab_Toolbox_for_Dimensionality_ Reduction.html
50 Advances in mapping distortion visualization for NLDR using Cartograms 4.4.1 Robust topographic mapping and its magnification factors One constraint of the basic GTM model is due to the fact that the centres of the mixture components do not move independently from each other, as they are limited by definition to lie in a low-dimensional embedded manifold. The basic GTM model presented in § 2.2.1, has some obvious limitations when dealing with atypical data or outliers, due to the narrowness of the tails of the Gaussian distributions. The point is that the presence of outliers is likely to bias the estimation of parameters Wand β, so other more robust formulations of GTM has been proposed using a mixture of Student’s t-distributions (the t-GTM model). This model, when applied to multivariate data clustering and visualization, provides a more accurate imputation of missing values and it is robust when dealing with outliers [Vellido et al., 2006; Vellido, 2006b]. To introduce the multivariate t-GTM model, we assume that the basis functions ϕmin the GTM mapping given by Eq. 2.5 are replaced by Student t-distributions. The conditional distribution of the data given the latent variables and the adaptive parameters is: p(y|x,W, βν) = Γν 2+p 2βp/2 Γν 2(νπ)p/2 1 + β ν p X j=1 yj− M X m=i w⊤ jϕm(x)!2 −(ν+p)/2 (4.10) where Γis the gamma function Γ(t) = Z∞ 0 xt−1e−xdx(4.11) and νis an adaptive parameter such that the multivariate t-distributions converges to a multivariate normal one when ν→ ∞. We can now integrate the latent variables out. After this, we obtain new expressions of the log-likelihood and again the model can be fitted to data using the EM algorithm to obtain, along with other results, the responsibilities rkn. For more details, see Vellido [2006a,b]. To update the parameters Wand β, we apply the maximization step of the algorithm by maximizing the expected log-likelihood, but a similar update cannot be applied to the parameter ν. To update ν, we have to consider alternative approaches such as, for instance, running experiments for a range of its possible values, selecting
4.4 Cartogram representations for GTM 51 the best choice out them. Svensèn [1998] gives an interpretation of the update expression of βas the offmanifold variance of the model being updated to the average weighted distance between original data and prototypes. This update formula, together with the one of the responsibility, implies the existence of a further weighting term for the t-GTM, which will be small for data outliers, as stated in [Peel and Mclachlan, 2000]. As a result, the influence of outliers on the model parameters will be effectively minimized [Vellido et al., 2006]. Magnification factors The t-GTM generates, as its standard counterpart, a varying local distortion that can make exploratory data visualization difficult. This distortion can again be quantified over the latent space continuum with MF, as explained in § 3.3. The analytical quantification of the MF can be expressed in terms of the derivatives of the basis functions ϕmas MF =pdet ((WΨ)⊤WΨ),(4.12) where J=WΨ is the Jacobian of the mapping transformation and Ψ∈RM×qhas elements ψmi defined as ∂ψm(x) ∂x(i)=Γ(ν+p 2)(−ν−p)βp+2 2 Γ(ν 2)πp/2νp+2 2x(i)−µ(i) m1 + β ν∥x−µm∥2−ν+p−2 2 (4.13) where µm,m= 1, . . . , M are the centres of the Student t-distributions. Remark Notice that the computation of the MF for the basic GTM can be done in the same way according to Eq. 4.12. The only change that need to be done is to replace the derivatives of the the Student-t functions with the derivatives of the Gaussian radial basis functions. 4.4.2 Cartogram visualization for t-GTM The following experiments compare the effect of outliers on the Cartogram representations of the MF for the standard GTM and for t-GTM. An artificial data set of
52 Advances in mapping distortion visualization for NLDR using Cartograms 3-D points is used to make possible the direct visualization of the data vectors yk in the observed data space. A total of 1,500 3-D points are randomly drawn from 3 spherical Gaussian distributions (500 points each), all with unit variance and with centres set at the vertices of an equilateral triangle. Two different subsets of outliers are added to this data set: •A-type): three outliers located on the normal to the imaginary plane defined by the cluster triangle that passes through its barycenter; •B-type): three outliers located on the normal to one vertex of the imaginary triangle. We choose a 15 ×15 regular grid of 2-D latent points. Both GTM and t-GTM are trained with the same initialization. The MF is calculated for both methods and Cartograms are generated using these values. In all cartograms, we assume assumed a homogeneous level of distortion in the space beyond the grid: this value is chosen to be equal to the mean distortion over the complete map 1/K PK k=1 MFxk, as described in section § 4.2. Likewise, we assume that the level of distortion within each of the squares associated to xkis itself uniform. The first experiment, displayed in Fig. 4.6, corresponds to the inclusion of Atype outliers, while the second, displayed in Fig. 4.7, corresponds to the inclusion of B-type outliers. Despite the fact that most GTM prototypes concentrate in the three clusters, it is clear from the image in Fig. 4.6 (top row, left) that, in the case of standard GTM, the A-type outliers force the manifold towards them in an undue manner. This causes a distortion that is more controlled by the outliers than by the empty space between clusters. Even though, the Cartogram visualizations generated by GTM and t-GTM are rather similar. The reason for this is the artificial symmetry of the outliers location. The maps in Fig. 4.7, corresponding to the second experiment with added B− type outliers, tell a very different story. Now, the symmetry is lost and the MF of the standard GTM reflects the fact that the model stretches one of the sides of the manifold in its attempt to cover the outliers (top row, left). As a result, an artefactual high distortion appears in the top-righ corner of the MF representation map (where outliers are seen to be mapped) and biases the Cartogram representation. The tGTM, instead, ignores the outliers and respects the symmetry of the representation while restricting the manifold to the imaginary triangle defined by the three clusters. This is clearly reflected in the corresponding Cartogram.
4.5 Discussion 53 Notice that, in both experiments, the extra MF distortion introduced by the outliers makes the data representation of the data of all clusters far more compact for GTM than for t-GTM. In any case, this simple preliminary experiments illustrate how modelling methods that behave robustly in the presence of outliers are more likely to produce more faithful representations of the non-linear mapping distortion and, as a result, more faithful data visualizations. 4.5 Discussion In this chapter we have presented a novel technique that allows to introduce into the visualization space a quantitative representation of the non-linear mapping distortion due to the chosen dimensionality reduction model. One advantage of this technique is its portability, since a Cartogram visualization is possible for any non-linear dimensionality reduction method in which a distortion measure can be defined. The experiments presented here are meant to be exploratory, aiming to asses the validity of the proposed model. Further investigation can be done by using more complex datasets, with different characteristics and different statistical properties. One example of this is presented in the recent work of Vellido et al. [2013], where the results over diverse artificial datasets provide some guidelines for the use of Cartograms in NLDR. It is shown that some underlying quantities, such as the dimensionality of the data features, affect the visualization space more than others, such as the density of data clusters or the size of the latent grid. It is important to stress that the use of Cartogram visualization does not affect or modify the performance of the model during its training. It is just an a posteriori method for the visual display of the results. The Cartograms provide added value to the visualization of the low-dimensional space by enabling a deeper analysis of the model capabilities. One example of this has been given in section 4.4, where we have shown that the Cartogram representation reflects the model capability of behaving robustly in the presence of atypical data.
54 Advances in mapping distortion visualization for NLDR using Cartograms −30 −20 −10 010 20 −20 0 20 40 −20 −15 −10 −5 0 5 10 15 20 GTM −30 −20 −10 010 20 −20 0 20 40 −20 −15 −10 −5 0 5 10 15 20 t−GTM 500 1000 1500 2000 500 1000 1500 2000 Fig. 4.6 Top row) Representation of the original data clusters (1,500 points, 500 in each cluster, plus A−type outliers, all represented with crosses) with standard GTM (left) and t-GTM (right) The generated manifold is superimposed; it is represented as a grid whose knots are the model prototypes yk; central row) Maps of MF values together with a colorbar for interpretation on the right-hand side of the maps; bottom row) Corresponding Cartograms, based on the MF, to which the mean projections of the data are superimposed. The mapping locations of outliers are highlighted with circles.
4.5 Discussion 55 −40 −20 020 40 −20 −10 0 10 20 30 −20 −15 −10 −5 0 5 10 15 20 GTM −40 −20 020 40 −20 −10 0 10 20 30 −20 −15 −10 −5 0 5 10 15 20 t−GTM 200 400 600 800 1000 1200 1400 1600 200 400 600 800 1000 1200 1400 1600 Fig. 4.7 Representation of data, manifold grid, MF maps and Cartograms for the second experiment with B−type outliers as in Fig. 4.6. Notice the difference in the mapping locations of outliers (again highlighted with circles) as compared to Fig. 4.6. In this case, GTM maps the outliers in a high-distortion area that is generated by the own outliers and not by the cluster data points (note that this high distortion appears because of just three outlier points), whereas the t-GTM maps them correctly to the closest cluster, without any artefactual distortion.
Chapter 5 Increasing Interpretability of MTS Modelling through Visualization Using Manifold Learning Time-dependent natural phenomena and artificial processes can often be quantitatively expressed as multivariate time series (MTS). As in any other process of knowledge extraction from data, the analyst can benefit from the exploration of the characteristics of MTS through data visualization. This visualization often becomes difficult to interpret when MTS are modelled using non-linear techniques. Despite their flexibility, non-linear models can be rendered useless if such interpretability is lacking. The methods described in previous chapters have mostly focused on static i.i.d. data. In this chapter, we model MTS using VB-GTM-TT, a variational Bayesian variant of a constrained hidden Markov model (HMM) of the manifold learning family defined for MTS visualization. We aim to increase its interpretability by taking advantage of two results of the probabilistic definition of the model: the explicit estimation of probabilities of transition between states described in the visualization space and the quantification of the non-linear mapping distortion. 5.1 Exploring MTS Most applied analysis of multivariate time series involves, in one way or another, problems with specific targets such as prediction, forecasting, or anomaly detection. A less explored avenue of research is the exploratory analysis of multivariate time series using machine learning and computational intelligence methods [Fu, 2011].
64 Increasing MTS Model Interpretability through Visualization Using Manifold Learning 0 200 400 600 800 1000 A B CD E A C D E 0 200 400 600 800 1000 0 5 10 15 20 25 Magnification Factors MF Magnification Factors 0 5 10 15 20 25 30 35 40 45 Fig. 5.2 Top left: Plot of Shuttle-data; the five intervals or regimes separated by sudden transitions are identified as A, B, C, D and E. Top right: State-membership map generated by VB-GTM-TT, with a 13 ×13 grid of hidden states represented as squares; the relative size of these squares is proportional to the time data points assigned to them; the starting point of the MTS is represented as a star, the ending point as a circle. Bottom Left: The Magnification Factors as a function of time, including the mean MF over all states (represented as a dashed line); narrow peaks of distortion are detected precisely in the areas of sudden transitions. Bottom Right: MF gray-shade color map, represented in the VB-GTM-TT latent space visualization grid; white areas correspond to high distortion. (bottom row, left): MF narrow spikes of varying magnitude (particularly strong in the transition from Bto C) appear in the transitions between time intervals. These spikes take values well over the mean MF of the map. This result suggests that the evolution of the MF over time could directly be used to detect sudden regime transitions in MTS. The CSTP maps in Fig. 5.3 are very consistent with their MF counterparts, and complement them. Alternatively displayed as 3−Dmaps over the grid of hidden states, they provide an intuitive illustration of the previously described behaviour.
5.4 Discussion 65 Following a geographical representation visual metaphor, the MTS can be seen to flow across cumulative state transition probability ridges, where rapid transitions between regimes see the MTS moving through relatively lower-valued depressions in those ridges. An opposite graphical metaphor could be used for the MF distortion, with the MTS flowing through its valleys, that is, across areas of the map characterized by low MF values. 0 1 2 0 1 2 Fig. 5.3 A 3−Drepresentation for the CSTP plots. The values in the vertical axis correspond to the CST P values over the latent space. Left: artificial data; right: Shuttle-data. 5.4 Discussion Data visualization can be of great assistance in knowledge extraction processes. High dimensionality is always a barrier for visualization. In the case of MTS, this is compounded by their i.i.d. nature, because the search for patterns over time is often relevant in their study. Dimensionality reduction can make visualization operative for high-dimensional MTS. The use of non-linear dimensionality reduction methods to this purpose poses a challenge of model interpretability due to the existence of locally-varying distortion. In this chapter, we have proposed to use MF and CSTP to improve interpretability for VB-GTM-TT, a manifold learning NLDR method. The model mapping distortion has been explicitly quantified in the latent space continuum and the probabilistic nature of the method has allowed us to define a cumulative probability of state transition. The reported experiments have shown that both metrics can provide interesting insights that enhance the low-dimensional visualization of the MTS provided by the model. This exploration approach is quite flexible and could be extended to other dimensionality reduction models for MTS analysis, provided their local distortion can be
66 Increasing MTS Model Interpretability through Visualization Using Manifold Learning quantified. Examples of this may include GP-LVM [Lawrence, 2005]), GP dynamical models (GPDM, [Wang et al., 2008; Damianou et al., 2011]) or temporal Laplacian eigenmaps [Lewandowski et al., 2010]. It could also be extended to alternative visual display methods, such as the Cartograms presented in chapter § 4, [García et al., 2013; Tosi and Vellido, 2013, 2012] and warped topographic maps [Gianniotis, 2013].
Chapter 6 Metrics for Probabilistic Geometries and Their Impact on Interpretability In many practical applications of multivariate data analysis, probabilistic models are combined with heuristic algorithms that use computed distances and interpolants between observations; this is, for instance, the case of distance-based methods for classification and clustering, methods of data interpolation for the reconstruction of missing frames (in computer vision problems), or techniques for the definition of the optimal path from one data point to another (in robotics, or cartography), to name a few. These heuristic algorithms are often based on the assumption of a deterministic system, without taking into account the uncertainty deriving from the probabilistic framework. Therefore, if the data observations are noisy, should the underlying geometry (and hence distances and interpolants) not be considered noisy as well? This chapter studies the role of Random Geometries, i.e. distributions of manifolds, in machine learning and, more generally, in data analysis. We develop methods for the estimation of distributions of Riemannian manifolds which support observations as well as algorithms for computing interpolants (geodesics) over distributions of manifolds. The new geometrical insight given to the interpretation of probabilistic modelling may not only theoretically advance statistical machine learning, but also improve the usability and flexibility of existing models.
68 Metrics for Probabilistic Geometries and Their Impact on Interpretability 6.1 Metrics for Probabilistic LVMs Manifold learning approaches attempt to learn the underlying support of the data (the manifold). Using the concepts of Riemannian geometry presented in section § 3.2, it is possible to derive the intrinsic geometrical properties of the model by explicitly computing its local metric tensor continuously over the input space. Once the metric has been derived, is then possible to compute geometrical quantities such as distances, angles, or the curvature of the space. Given that manifold learning models can be of different nature, we have decided to restrict the developments in this chapter to smooth generative models. In such models, the local metric varies smoothly across the input space, in contrasts with prior approaches such as those reported in [Bregler and Omohundro, 1994; Tenenbaum, 1997; Tenenbaum et al., 2000], which use metrics that vary discretely across the space. Other relevant approaches related to Gaussian models can be found in [Lawrence, 2012]. In the following sections, the local metric tensor for generative latent variable models is first defined. We then illustrate the specific cases of GP-LVM and GTM, providing two algorithm to compute shortest paths (geodesics). The novelty of this approach is the probabilistic expression of the local metric, which opens to new streams of investigation in the field of probabilistic geometries for latent variable models. 6.2 The distribution of the natural metric Given a noise model as in in Eq. (2.1), we assume the mapping fbetween the latent space and the observed space to be a differentiable function. This appears not to be a strict requirement, since some of the most used models fulfill it. Two examples are the standard GTM and the GP-LVM with the exponentiated quadratic kernel (EQ), since in both cases the differential of freduces to the differential of an exponential. Under this assumption of smoothness, the output of the mapping can be interpreted as a differential manifold (c.f. section § 3.2). It is then possible to explicitly compute the natural Riemannian metric of the given model as follows: let Jbe the Jacobian of the mapping fgiven as in Eq. (1.1), then the tensor G=J⊤J defines a local inner product structure over the latent space according to Eq. 3.3.
6.2 The distribution of the natural metric 69 Being the Gaussian distribution widely used in latent variable models, it is practical to analyse the particular case in which the conditional probability over the Jacobian follows a Gaussian distribution as well. The distribution over the Jacobian is inducing a distribution over the local metric tensor Gin a natural way. In fact, assuming independent rows of J, the distribution of the Jacobian is the product of p multivariate Gaussians p(J|X,f, β) = p Y j=1 N(Jj,:|µJj,:,ΣJ).(6.1) The resulting random variable follows a non-central Wishart distribution [Anderson, 1946] G=Wq(p, ΣJ,E[J⊤]E[J]),(6.2) where prepresents the number of degrees of freedom; the quantity Σ−1 JE[J⊤]E[J] is known as the non-centrality matrix and it is equal to zero in the central Wishart distribution. Intuitively, we can interpret the Wishart distribution as a multivariate generalisation of the Gamma distribution. It is interesting to observe that the expectation of the metric tensor is given by the sum of two terms, a mean term and a covariance term E[J⊤J] = E[J⊤]E[J] |{z } mean term +p·ΣJ |{z} covariance term ,(6.3) whose role will be made more explicit at the end of the section. Given the general formulation of the distribution of the Riemannian metric tensor in generative latent variable models, we can now provide the explicit expression for the models of our interest. GP-LVM local metric As described in section §2.3, a Gaussian process (GP) can be used in dimensionality reduction to describe distributions over a mapping f f(x)∼ GP(m(x), k(x,x′)) yi,j =Kfi,fKY:,j +ϵi, leading to the formulation of the Gaussian process latent variable model (GP-LVM) Lawrence [2005].
70 Metrics for Probabilistic Geometries and Their Impact on Interpretability It follows from Eq. 2.36 and the properties of the GPs that the distribution of the Jacobian of the GP-LVM mapping is the product of pindependent Gaussian distributions (one for each dimension of the dataset) with mean µJ(j,:) and covariance ΣJ. For every latent point x∗, the Jacobian takes the following form: p(J|Y,X,x∗) = p Y j=1 N(Jj,:|µJj,:,ΣJ)(6.4) = p Y j=1 N(∂K⊤ f∗,f˜ K−1 f,fY:,j, ∂2Kf∗,f∗−∂K⊤ f∗,f˜ K−1 f,f∂Kf∗,f), which (c.f. Eq. 6.2) gives a distribution over the metric tensor G G=Wq(p, ∂2Kf∗,f∗−∂K⊤ f∗,f˜ K−1 f,f∂Kf∗,f,E[J⊤]E[J]).(6.5) From this distribution, the expected metric tensor can be computed as E[J⊤J] = E[J⊤]E[J] + p ∂2Kf∗,f∗−∂K⊤ f∗,f˜ K−1 f,f∂Kf∗,f | {z } covariance term .(6.6) Magnification Factors The metric tensor defines the local geometric properties of the GP-LVM model and it can be used as a tool for data exploration that helps increasing the model interpretability. One way to visualise the tensor metric is through the differential volume of the high dimensional parallelepiped spanned by GP-LVM; this, for a latent dimension q= 2 is known as the Magnification Factor (MF), see section § 3.3. The explicit formulation of the MF for GP-LVM is given by MF = pdet (E[J⊤J]).(6.7) An illustrative example of MF visualization, using the Shuttle dataset presented in section § 5.3.1, is displayed in Fig. 6.1. The colormap is computed over a fine regular grid defined over the latent space. An additional example has been displayed in the background section § 3.3, using the jogging motion of the CMU motion capture database described in section § 2.3.3. Further examples are provided in the experimental results of this chapter.
6.3 Computing geodesics 71 MF 5 10 15 20 25 30 35 40 45 50 55 Fig. 6.1 Magnification Factor colormap computed over the GP-LVM latent space using the Shuttle dataset described in section § 5.3.1. 6.3 Computing geodesics Given a latent space endowed with an expected Riemannian metric, we now consider how to compute geodesics (shortest paths) between given points. Once a geodesic is computed, its length can be evaluated through the numerical integration of Eq. 3.6. An obvious solution to the shortest path problem is to discretise the latent space and compute shortest paths on the resulting graph using, e.g., Dijkstra’s algorithm [Cormen et al., 1990]. The computational complexity of this approach, however, grows exponentially with the dimensionality of the latent space and the approach quickly becomes unfeasible. Moreover, this approach (presented in section § 6.3.1) will also introduce discretisation errors due to the finite size of the graph. Instead, we propose solving the geodesic differential equation (3.7) numerically. This scales more gracefully as it only involves a discretisation of the geodesic curve which is always one-dimensional independently of the dimension of the latent space. This approach is presented in section § 6.3.2.
72 Metrics for Probabilistic Geometries and Their Impact on Interpretability Remark Note that the expectation of the GP-LVM metric tensor includes a covariance term depending on the covariance function of the GP prior. This implies that the metric tensor expands as the uncertainty over the mapping increases. Hence, curve length also increases when traversing uncertain regions and, as a consequence, geodesics will tend to avoid these regions. This nice effect comes from on the fact that in GPLVM the covariance term depends on the values of the latent mapping. Not all the latent variable models have this property. For instance, in GTM (c.f. Eq. 2.10) the covariance of the data posterior distribution is constant and equal to β−1. It follows that the derivatives of the covariance terms with respect to the latent space are equal to zero, which means that also the covariance term which appears in the formulation of the expected metric in Eq. 6.3 is equal to zero. This observation shows that the formulation of the expected metric tensor for the GTM reduces to the original formulation of the metric given by Bishop et al. [1997a], which is computed without taking into account the covariance term. In the case of GTM, this omission makes no difference because the covariance of J is constant, due to the fact that the covariance of the conditional distribution given by Eq. 2.10 does not depend on the latent space. 6.3.1 Geodesics via discretisation As a first approach, we propose [Tosi and Vellido, 2014] to discretise the space into a square grid of nodes ν(i,j)∈Rqand build a weighted graph where every node is connected to its eight neighbours. A graphical visualisation of the connections of ν(i,j)is the following ν(i−1,j−1) ν(i−1,j)ν(i−1,j+1) .... . . ... ν(i,j−1) · · · ν(i,j)· · · ν(i,j+1) ... . . .... ν(i+1,j−1) ν(i+1,j)ν(i+1,j+1) The weights ω(i,j)of the graph are computed according the local value of the MF
6.3 Computing geodesics 73 evaluated in the node ν(i,j) ω(i,j)=MFx∗=pdet ((WΨ(x∗))⊤WΨ(x∗)),x∗=ν(i,j).(6.8) Notice that the MF expression presented in Eq. 4.12 for the t-GTM can be used for the computation of the MF in GTM by replacing the derivatives of the Student-t functions with the derivatives of the Gaussians. The distance between any two nodes is then defined as the shortest path over the graph using Dijkstra’s algorithm1[Cormen et al., 1990]. We display an illustrative example using 3-D data points sampled from a spiral. A standard GTM is trained wit it in order to learn a 2-D representation of the given data (see Fig. 6.2). We show in Fig. 6.3 that geodesic distances computed as shortest paths over the graph provide a more faithful interpolation between data points: in fact, the resulting interpolating path follows the natural structure of the data, giving a more faithful distance than to the Euclidean straight line. Fig. 6.2 A 3-D artificial dataset (a spiral) is used for training. The Generative Topographic Mapping (GTM) is used here as an illustrative technique; the model is trained over a 15×15 grid of nodes, laying on a 2-D latent space (left); training data are represented as blue dots, connected by a continuous line. The GTM latent grid is projected into the 3-D observed space to visualise the embed of the observations. The computational complexity of the algorithm is that of the Dijkstra algorithm. To achieve a better accuracy in the geodesic computation, we can consider a finer grid, but this results in a growing computational cost. 1The algorithm has been implemented using the function graphshortest in Matlab®.
80 Metrics for Probabilistic Geometries and Their Impact on Interpretability MF 500 1000 1500 2000 2500 Fig. 6.9 GP-LVM latent space representation for the motion capture data. White dots denote latent points xn, whereas the background colour is proportional to the values of the MF (6.7), according to the code expressed by the colorbar. It is quite obvious that latent points are spread across paths of low MF, whereas latent space areas of high MF are avoided. 6.5 Discussion When the mapping between a latent space and the observation space is not isometric (the common case for non-linear mappings), a Euclidean distance measure in the latent space does not match that of the original observation space. In fact, the distance measures in the latent and observation spaces can be arbitrarily different. This makes it difficult to perform any meaningful statistical operation directly in the latent space as the used metric is difficult to interpret. We solve this issue by carrying the metric from the observation space into the latent space in the form of a random Riemannian metric. This gives a distribution over a smoothly changing local metric at each point in the latent space. We then provide an expression for the expected local metric and show how shortest paths (geodesics) can be computed numerically under the resulting metric. These geodesics provide natural generalisations of straight-lines and they are, as a result, suitable for interpolation under the new metric.
6.5 Discussion 81 DataPoints Geodesic Euclidean Distance Fig. 6.10 GP-LVM latent space for the motion capture data. Black dots denote latent points xn. The green curve denotes the geodesic interpolant, while the dashed brown curve is the straight-line interpolant. This figure is best viewed in colour. 20 22 24 26 28 30 Geodesic Euclid Latent Real Fig. 6.11 Length, in centimetres, of the subjects forearm during latent space interpolation. The blue curve is according to the geodesic interpolant, and the red dashed curve is according to the straight-line interpolant. For reference, the black dots show the true length. For the GP-LVM model, the expected metric depends on its uncertainty, in such
82 Metrics for Probabilistic Geometries and Their Impact on Interpretability a way that distances become longer in regions of high uncertainty. This effectively forces geodesic curves to avoid uncertain regions in the latent space, which is the desired behaviour for most applications. It is worth noting that a similar analysis for the GTM does not provide a metric with this capacity as the uncertainty is constant in this model. The idea of considering the expected metric is practical as it turns the latent space into a Riemannian manifold and this opens up to many applications. E.g. tracking can be performed in the latent space through a Riemannian Kalman filter [Hauberg et al., 2013], classification can be done using the geodesic distance, etc. It is, however, potentially misleading to only consider the expectation of the metric rather than the entire distributions of metrics. Although, if the latent dimension is much lower than the observed data dimension, it can be shown that the distribution of the metric concentrates around its mean. But, in general, random Riemannian manifolds are mathematically less well-understood, e.g. it is known that geodesics are almost surely not length minimising curves under a random metric [LaGatta and Wehr, 2014]. We are suggesting that manifolds derived from data are necessarily uncertain, and there is much to gain from further consideration of these spaces, which then naturally lead to distributions over geodesics, distances, angles, curvature and so forth. In this chapter, we have only considered how geometry can be used to understand an already estimated LVM, but it is also worth considering if this geometry can be used as part of the LVM estimation. That is, it would be worth investigating if a prior on the curvature of the latent manifold is an effective way to influence learning.
6.5 Discussion 83 Fig. 6.12 Example poses from the motion capture data. These poses are temporarily spaced between the end-points of the interpolating curves, i.e. they are comparable to the interpolated reconstructions. Fig. 6.13 Interpolated poses according to the straight-line interpolant. In particular, note the bending of the knees, which does not occur in the training data. Fig. 6.14 Interpolated poses according to the geodesic. These are visually similar to the poses in Fig. 6.12.
Chapter 7 Conclusions The final chapter of the thesis aims to summarise the main contributions presented in previous chapters and sketch some avenues for future research. 7.1 Summary of the thesis and its main contributions In this thesis, we have addressed the problem of visualization of high-dimensional data sets, with the main objective of improving the interpretability (and as result, the usability) of the probabilistic non-linear dimensionality reduction models used to generate such visualization. To do so, we have mainly, but not only, exploited the intrinsic geometrical structure of the aforementioned models. The main contribution in the domain of visualization techniques is given by the Carotgram-based representation, presented in chapter § 4. This novel technique, inspired by cartographic maps in the geography domain, has been used to reintroduce in the visualization space a loss of information in which the non-linear mapping incurs. The analytical quantification of such a distortion has been expressed in the form of Magnification Factors, and then computed and visualised together in the form of the Cartogram maps. Results for the case of Self Organizing Maps and Student-tGenerative Topographic Mapping have been presented. The research carried out in this thesis can be applied to multivariate data of different nature and from different domains. In chapter § 5 we presented experimental results for multivariate time series. We improved interpretability for the VB-GTM-TT model and we introduced the explicit estimation and visualisation of the cumulative probabilities of transition between states CSTP. The main theoretical contribution in the domain of Random Geometries is
86 Conclusions provided in the last chapter, §6. There we reinterpreted latent variable models as Riemannian manifolds, by pulling back the Riemannian metric from the observation space. Differently from previous approaches to metric learning, our local metric is defined in a probabilistic way, providing an explicit expression of its probability distribution for the considered models. Experimental results have shown that inference made following the Riemannian metric leads to a more faithful generation of new data. In addition, we have stated that the algorithms described in this thesis can be extended to other generative latent variable models characterised by a smooth mapping between the latent space and the observed space. This property points out to the portability of the proposed approach. In conclusion, we have carried out an extensive analysis of the problem of interpretability in probabilistic dimensionality reduction, from the differential geometry point of view. We hope that this work will be of help to other researches with related focus and open the way to novel investigation, such as the one suggested in the next session. 7.2 Open questions and future directions Random geometries We have shown how to define a distribution over the metric in latent variable models. In particular, this was achieved by pulling back the metric from the observation space into the latent space in the form of a random Riemannian metric. This research opens to new streams of investigation in the field of Random Geometries in relation to machine learning. Random Riemannian manifolds are mathematically less well-understood than Riemannian manifolds. In fact, it is known that geodesics are almost surely not length minimising curves under a random metric [LaGatta and Wehr, 2014]. We are suggesting that manifolds derived from data are necessarily uncertain, and there is much to gain from further consideration of these spaces, which then naturally lead to distributions over geodesics, distances, angles, curvature and so forth. In this thesis, we have only considered how probabilistic geometry can be used to understand an already estimated generative latent variable model. This work opens the way to promising direction of investigation, namely the applicability of probabilistic geometry as part of the model estimation itself and, if so, it is worth understanding its influence in the learning process.
7.2 Open questions and future directions 87 Distributions of geodesics We have developed a probabilistic framework where the support of the data can be interpreted as a random Riemannian manifold and geodesic distances can be computed by taking the uncertainty of the metric into account. The preliminary results presented in chapter § 6 give rise to some questions, such as: (1) How does uncertainty defined over the the metric tensor impact the geodesic distances computed in the observed space? (2) What are the conditions of existence and uniqueness of geodesics in a more general setting? (3) How do we analytically define, if exists, a distribution over geodesics? (4) What is the behaviour of the lengths of geodesics when the number of features grows? Do distances concentrate as the dimensionality of the feature space goes to infinity? In order to answer these questions, we are currently investigating how to develop an algorithm to explicitly compute distributions over geodesics in probabilistic dimensionality reduction. One approach to be considered is that which entails combining the local distributions of the metric tensor for a given set of points by considering a joint sample: this way we would obtain samples of the metric along the whole manifold (i.e. samples from the distribution of the random manifold). We display examples of these samples in Fig. 7.1, where the diagrams refer to the 3-dimensional visualization of the MF computed over a 2-dimensional GP-LVM latent space (the 3-D visualisation has been introduced in section §3.3). Here the first diagram represents the plot of the MF of the expected metric tensor, defined according to Eq. 6.6; the rest of diagrams show some random joint samples from the Wishart distribution of the metric tensor (rather than showing just the mean). Once the samples of the manifold are computed, we can compute (for each sample) the samples of geodesics using the algorithms proposed in section §6.3.1 and section §6.3.2. This is straightforward for the GP-LVM model and a Wishart distributed metric introduced above, and the results can be extended to similar generative models. This will provide new theoretical insights into the problems of probabilistic geometries, as the conditions of existence and uniqueness of geodesics in a more general setting (i.e. wider class of models) are not well known. One challenge is the fact that, in general, the distribution of the geodesic can be very complex because even small changes in the metric tensor can result in a big change in the geodesic and its length. Moreover, if the dimensionality of the feature space is increased and sent to infinity, this is thought to have an effect on the lengths of interpolating paths, resulting in an effect of concentrated distances. This aspect represents one of the
88 Conclusions ↙↓↘ Fig. 7.1 On the top plot: mean of the distribution of a random manifold (generated after training a GP-LVM over a jumping jacks motion form the CMU database, c.f. Experimental results in section §6.4.3). Three different joint random samples are generated from this Wishart, and values of the MF factors are represented by a colormap. A 3-D visualization of the MF is used: the vertical axes of each plot represents the values of the MF. future directions of research. Probabilistic numerics New sparks of research aim to identify numerical methods as learning problems using probabilistic models. This approach is know as Probabilistic numerics1and addresses classical optimisation algorithms and numerical methods for the solution of differential equations and integrals. Such probabilistic numerical methods, applied to find solutions of ordinary differential equations (ODEs), have an impact on the analysis of statistical Riemannian manifolds [Hennig and Hauberg, 2014]. In particular, the very recent work of Schober et al. [2014] provides a probabilistic model for the solution of ODEs which matches the classical Runge-Kutta method. A future direction of research is to investigate the combination of these recent advances in probabilistic numerics and the propagation of the uncertainty defined over the metric tensor through the geodesic ODE solver presented in section §6.3.2. Big data In this thesis, we have addressed the issue of high-dimensional data as a problem of datasets with a high number of features. But we also need to consider that, given the ever increasing amount of available data generated on daily basis in different 1About: http://probabilistic-numerics.org/
7.2 Open questions and future directions 89 domains, we have at our disposal a growing amount of datasets which are big mostly in terms of number of observations. This has come to be popularised under the name of Big data ant constitutes a very up-to date problem. Performing inference with Gaussian processes-based techniques suffers from a high computational cost, and scaling up GPs is a topic of ongoing research. The complexity can be reduced using appropriated approximation techniques as well as appropriate distributed algorithms. Very recent results [Hensman et al., 2013; Gal et al., 2014] apply GPs to data of the order of N∼106. Following these promising results, we aim to extend the approach presented in this thesis to the analysis of larger datasets. Intrinsic dimensionality of the dataset While performing dimensionality reduction, we have mostly set the dimensionality of the latent space to be q= 2, in order to display the experimental results on the visualisation space. This has been done because one of the scopes of this thesis is to make progress in the visualization techniques in order to improve the interpretability of the considered models. The theoretical results presented, however, are valid for any choice of latent dimension q≤p. By removing the restriction of a 2-D latent space, we face the problem of how to set the value of q. This question has been answered for a certain class of models, and in the context of GP based dimensionality reduction has been solved by the Varitional GP-LVM [Titsias and Lawrence, 2010; Damianou et al., 2014]. In this model the algorithm is capable of computing the intrinsic dimensionality of the data by optimising the lengthscales of the kernel in each latent dimension (and, eventually, switching off the non relevant ones). Extension to other models The focus of this work is mainly on two models of interest, the Generative Topographic Mapping and Gaussian Process Latent Variable Model, both part of the family of generative latent variable models. We have explicitly defined a distribution over the local metric and we have visualised such metric using the values of Magnification Factors. The conditions to extend our approach to a wider class of probabilistic models is to have a smooth mapping between the latent space and the observed space. This is the case of GP-based models that use differentiable kernel functions, which opens to a wide class of models. Examples of a straightforward extension have been done in this thesis using the GP dynamical system (GPDM,
96 References I. Olier and A. Vellido. Advances in clustering and visualization of time series using GTM through time. Neural Networks, 21(7):904–913, 2008c. K. Pearson. On lines and planes of closest fit to points in space. The London, Edinburgh and Dublin Philosophical Magazine and Journal of Science, 2:559–572, 1901. D. Peel and G. J. Mclachlan. Robust mixture modelling using the t distribution. Statistics and Computing, 2000. J. Pointer. The cortical magnification factor and photopic vision. Biological Reviews, 61(2):97–119, 1986. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipies In C: The Art of Scientific Computing. Cambridge University Press, Cambridge, England, 1988. L. R. Rabiner. A tutorial on hidden Markov models and selected applications in speech recognition. In Proceedings of IEEE, volume 77, pages 257–286. IEEE, 1989. C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. Cambridge, MA, 2006. B. Riemann. On the Hypotheses Which Lie at the Foundations of Geometry. 1854. S. T. Roweis. EM algorithms for PCA and SPCA. In M. I. Jordan, M. J. Kearns, and S. A. Solla, editors, Advances in Neural Information Processing Systems (NIPS), pages 626–632. The MIT Press, 1997. S. T. Roweis and L. K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science, 290(5500):2323–2326, 2000. M. Schober, D. Duvenaud, and P. Hennig. Probabilistic ODE solvers with RungeKutta means. In Advances in Neural Information Processing Systems (NIPS). In press, 2014. B. Schölkopf, A. J. Smola, and K.-R. Müller. Kernel principal component analysis. In Proceedings 1997 International Conference on Artificial Neural Networks, ICANN’97, page 583, Lausanne, Switzerland, 1997. E. Solak, R. Murray-Smith, W. E. Leithead, D. J. Leith, and C. E. Rasmussen. Derivative observations in Gaussian process models of dynamic systems. In S. Becker, S. Thrun, and K. Obermayer, editors, Advances in Neural Information Processing Systems (NIPS), pages 1033–1040. MIT Press, 2002. M. Svensén. GTM: The Generative Topographic Mapping . PhD thesis, Aston University, 1998. J. B. Tenenbaum. Mapping a manifold of perceptual observations. In M. I. Jordan, M. J. Kearns, and S. A. Solla, editors, Advances in Neural Information Processing Systems (NIPS). The MIT Press, 1997. J. B. Tenenbaum, V. Silva, and J. Langford. A global geometric framework for nonlinear dimensionality reduction. Science, 290(5500):2319–2323, 2000.
References 97 M. Tenenbaum and H. Pollard. Ordinary Differential Equations. Dover Publications, 1963. M. E. Tipping and C. M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society, 6(3):611–622, 1999. M. K. Titsias and N. D. Lawrence. Bayesian Gaussian process latent variable model. In Proceedings of the 13th international Conference on Artificial Intelligence and Statistics (AISTATS), volume 9, pages 844–851, 2010. W. R. Tobler. Thirty-five years of computer cartograms. Annals of the Association of American Geographers, 94:58–73, 2004. A. Tosi and A. Vellido. Cartogram representation of the batch-SOM magnification factor. In The 20th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN), pages 203–208, Bruges, Belgium, 2012. A. Tosi and A. Vellido. Robust cartogram visualization of outliers in manifold learning. In The 21th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN), pages 555–560, Bruges, Belgium, 2013. A. Tosi and A. Vellido. Local metric and graph based distance for probabilistic dimensionality reduction. In Proceedings of the Workshop on Features and Structures (FEAST 2014) International Conference on Pattern Recognition (ICPR 2014), Stockholm, Sweden, 2014. A. Tosi, S. Hauberg, A. Vellido, and N. D. Lawrence. Metrics for probabilistic geometries. In J. T. Nevin L. Zhang, editor, Proceedings of the 30th Conference on Uncertainty in Artificial Intelligence (UAI), pages 800–808, Quebec City, Canada, 2014a. AUAI Press Corvallis, Oregon. A. Tosi, I. Olier, and A. Vellido. Probability ridges and distortion flows: Visualizing multivariate time series using a variational Bayesian manifold learning method. In Advances in Self-Organizing Maps, the 10th International (WSOM), Advances in Intelligent Systems and Computing, pages 55–64. Springer, 2014b. A. Ultsch. U*-matrix: a tool to visualize clusters in high dimensional data. Technical Report 36, Philipps-University Marburg, Germany, 2003. A. Vellido. Assessment of an unsupervised feature selection method for generative topographic mapping. In S. D. Kollias, A. Stafylopatis, W. Duch, and E. Oja, editors, The International Conference on Artificial Neural Networks, volume 4132 of Lecture Notes in Computer Science, pages 361–370. Springer, 2006a. A. Vellido. Missing data imputation through GTM as a mixture of t-distributions. Neural Networks, 19(10):1624–1635, 2006b. A. Vellido, W. El-Deredy, and P. J. G. Lisboa. Selective smoothing of the generative topographic mapping. IEEE Transactions on Neural Networks, 14(4):847–852, 2003. A. Vellido, P. J. G. Lisboa, and D. Vicente. Robust analysis of MRS brain tumour data using t-GTM. Neurocomputing, 69(7-9):754–768, 2006.
98 References A. Vellido, J. D. Martín, F. Rossi, and P. J. G. Lisboa. Seeing is believing: The importance of visualization in real-world machine learning applications. In The 19th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN), 2011. A. Vellido, J. D. Martín, and P. J. G. Lisboa. Making machine learning models interpretable. In The 20th European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning (ESANN), pages 163–172, 2012. A. Vellido, D. L. García, and Nebot. Cartogram visualization for nonlinear manifold learning models. Data Mininig and Knowledge Discovery, 27(1):22–54, 2013. J. M. Wang, D. J. Fleet, and A. Hertzmann. Gaussian process dynamical models for human motion. IEEE Transactions on Pattern Recognition and Machine Intelligence (PAMI), 30(2):283–298, Feb. 2008. H. Wässle, U. Grünert, J. Röhrenbeck, and B. Boycott. Retinal ganglion cell density and cortical magnification factor in the primate. Vision Research, 30(11):1897– 1911, 1990. K. Q. Weinberger and L. K. Saul. An introduction to nonlinear dimensionality reduction by maximum variance unfolding. In The Conferece for the Association for the Advancement of Artificial Intelligence (AAAI), pages 1683–1686. AAAI Press, 2006.