# Michael Christoph Thrun

# Projection-Based Clustering through Self-Organization and Swarm Intelligence

Combining Cluster Analysis with the Visualization of High-Dimensional Data

Projection-Based Clustering through Self-Organization and Swarm Intelligence Michael Christoph Thrun

# Projection-Based Clustering through Self-Organization and Swarm Intelligence

Combining Cluster Analysis with the Visualization of High-Dimensional Data

Michael Christoph Thrun Marburg, Germany

Philipps-Universität Marburg 2017, Hochschulkennziffer 1180

ISBN 978-3-658-20539-3 ISBN 978-3-658-20540-9 (eBook) https://doi.org/10.1007/978-3-658-20540-9

Library of Congress Control Number:<PHONE_NUMBER>

#### Springer Vieweg

© The Editor(s) (if applicable) and The Author(s) 2018. This book is an open access publication.

**Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, express or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Printed on acid-free paper

This Springer Vieweg imprint is published by Springer Nature The registered company is Springer Fachmedien Wiesbaden GmbH The registered company address is: Abraham-Lincoln-Str. 46, 65189 Wiesbaden, Germany

# **Acknowledgments**

My gratitude goes to my sister Monika Sikora. You have converted with great devotion my often complex phrases and contexts into intelligible words. I thank you for your intuitive understanding of language.

I would like to thank Prof. Dr. Alfred Ultsch for his demanding scientific guidance and continuing education. His leadership and preparatory work have provided me with the tools to meet the demands of the job. Your suggestions have given me the creativity that was often necessary to redirect my research when a solution was not obvious.

Without my student colleagues, Felix Pape and Florian Lerch, some of the ideas sprouted in this work would not have been feasible. I owe you, Felix, because of your selfless voluntary commitment to the realization of manufactured 3D printing models of the U-matrices. The fundamental code for visualizing the U-matrices in the form of a topographic map would not have been possible without your cooperation, for which I thank you, Florian.

The positive working environment, which allowed me to flourish, I largely owe to my colleague Catharina Lippmann. Your collegiality and constructive cooperation created the atmosphere that was continuously aiding my systematic research.

# **List of figures**

# **List of tables**

# **Zusammenfassung**

Die vorliegende Arbeit befasst sich mit einem neuen Ansatz zur Clusteranalyse hochdimensionaler Daten. Die projektionsbasierte Clusteranalyse verbindet in zwei Dimensionen erhaltenen Strukturen mit zugrunde liegenden hochdimensionalen Strukturen.

Hierbei werden Cluster als natürlich definiert, wenn sie auf hochdimensionalen Daten beruhen, welche Diskontinuitäten aufweisen. Solche distanz- oder dichtebasierte Diskontinuitäten bezeichnen entweder kompakte oder verbundene Strukturen. Natürliche Cluster mit kompakten Strukturen werden hauptsächlich durch Inter- und Intra-Cluster-Distanzen definiert, während verbundene Strukturen auf dem Prinzip von Nachbarschaften zwischen Datenpunkten beruhen.

Diese Arbeit geht der Fragestellung nach, wie man einen korrekten Typ von Strukturen herausfinden kann, der Cluster in einem hochdimensionalen Datensatz ohne Vorannahmen definiert. Es wird dargelegt, dass Verfahren der Dimensionsreduktion helfen können, dieses Problem zu lösen.

Projektionsverfahren stellen einen gängigen Ansatz zur Dimensionalitätsreduktion hochdimensionaler Daten dar. Sie werden verwendet, um die Größe des Eingaberaumes zu reduzieren um dadurch eine Visualisierung der hochdimensionalen Daten zu ermöglichen.  Die Projektion kann zu einer irreführenden Interpretation der Strukturen führen. Die Qualitätsmaße (QM) zur Bewertung der Projektion haben Schwierigkeiten Diskontinuitäten in hochdimensionalen Daten korrekt zu erfassen, weil sie unter Umständen auf falschen Annahmen über die zugrunde liegenden hochdimensionalen Strukturen basieren. Andernfalls könnte mittels einer QM eine globale Zielfunktion definiert werden. Es wäre damit immer möglich, eine strukturerhaltende Projektion durch Optimierung dieser Zielfunktion zu erhalten.

Das aus diesen drei Modulen bestehende Verfahren Databionicswarm (DBS) wird in dieser Arbeit vorgestellt. Das erste Modul des hier vorgeschlagenen Ansatzes besteht darin, hochdimensionale Distanzen in der zweidimensionalen Projektion durch eine dreidimensionale topographische Karte mit hypsometrischen Farben zu visualisieren.

Im zweiten Modul wird das neue Projektionsverfahren Pswarm vorgeschlagen. Pswarm nutzt die Konzepte der Schwarmintelligenz, Selbstorganisation, Symmetrieüberlegungen der Physik und das Nash-Gleichgewichtskonzept aus der Spieltheorie.  Dieses Projektionsverfahren erfordert, abgesehen von der Distanz, keine Eingabeparameter für die Projektion. Durch Selbstorganisation können Strukturen von hochdimensionalen Daten durch einen Prozess abgebildet werden, der als Emergenz bekannt ist.  Pswarm wurde mit den üblichen Projektionsmethoden PCA, CCA, t-SNE, ESOM, NeRV und dem MDS-Technik-Sammon-Mapping verglichen. Hierbei wurde ein neues Qualitätsmaß (Delaunay Classification Error, DCE) eingesetzt. Der DCE ermöglicht durch die Verwendung vorgegebener Klassifikationen eine unvoreingenommene Beurteilung der Projektionsqualität für beide Arten von Strukturen. Die Ergebnisse zeigen, dass es mit Pswarm-Projektionen möglich ist Projektionen resultierend aus der Optimierung einer globalen Zielfunktion zu übertreffen.

DBS übertrifft die gängigen Methoden der Clusteranalyse (k-means, PAM, Single-Linkage, Spektralclustering, modellbasierte Clustering und Ward) hinsichtlich Stabilität und Plastizität auf einem künstlichen Benchmark-System von Datensätzen (FCPS). Im Gegensatz zu anderen üblichen Methoden der Clusteranalyse findet DBS keine Cluster, wenn keine natürlichen Cluster vorhanden sind. Die Anzahl der Cluster kann hierbei mit Hilfe einer Visualisierung abgeschätzt werden.

Die Anwendung von DBS auf drei hochdimensionale und multivariate Datensätze für den praktischen Gebrauch (Leukämie, Welt-Bruttoinlandsprodukt, Tetragonula-Bienen) reproduzierten bereits bekannte Erkenntnisse. In zwei aktuellen Anwendungen, Hydrologie und Schmerz-Gene findet DBS plausible und erklärbare Cluster.

   Darüber hinaus können 3D-Drucke der visualisierten Strukturen von hochdimensionalen Datensätzen mit üblichen 3D-Drucktechniken hergestellt werden.

# **Abstract**

This work introduces a new approach for cluster analysis defined as projection-based clustering. The projection based clustering combines structures preserved in two dimensions with underlying high-dimensional structures, if natural clusters exist in high-dimensional data. Clusters are defined as natural, if they are based on patterns in high-dimensional data characterized by discontinuity. Discontinuous patterns, which can either be based on distance or density, are described in this work as compact or connected structures. Natural clusters with compact structures are defined mainly by inter- versus intracluster distance, whereas the connected structures are based on the idea of neighborhoods present between data points.

With the use of basic principles founded on graph theory, this work demonstrated that the objective functions of clustering and visualization are based on the fundamental distinction between connected and compact structures. The derived conclusion is that in a case when the goal is to achieve a structure-preserving visualization or clustering, the optimization of a mathematical objective function could yield misleading results if the underlying structures of the highdimensional data do not coincide with the objective function. The question that arises is how to recognize structures that defines clusters in a high-dimensional data set without prior knowledge. The argument here is that dimensionality reduction methods may help solve this problem.

Projections are common dimensionality reduction methods to visualize high-dimensional data in a two-dimensional space. However, when restricting the Output space into two dimensions resulting in a two dimensional scatter plot (projection) of the data, low dimensional similarities do not represent high dimensional distances coercively. This could lead to a misleading interpretation of the underlying structures. Further, it is argued here that the quality measures (QMs), which evaluate this projection, have difficulties to correctly grasp discontinuities in high-dimensional data; this is because they imply assumptions about the underlying high-dimensional structures. Otherwise, a global objective function could be defined using the best QM, and it would always be possible to obtain a structure-preserving projection or clustering by optimizing this objective function.

Therefore, the first module for a solution proposed here is to visualize high-dimensional distances in the projection through a three dimensional topographic map with hypsometric colors, which is a further development of the generalized U-matrix.

After an extensive review of application of artificial intelligence in data science, two interesting concepts are addressed here, called self-organization and swarm intelligence. The irreducible structures of high-dimensional data can emerge through self-organization in a phenomenon called emergence. If properly applied through the use of a swarm of intelligent agents, the datadriven approach presented in this work can outperform the optimization of a global objective function in the tasks of clustering and dimensionality reduction.

Here, the second module called Pswarm, is presented for projecting high-dimensional data. Pswarm exploits the concepts of swarm intelligence, self-organization, symmetry considerations in physics, and the Nash equilibrium concept from game theory. It eliminates the need for a global objective function and does not require any input parameters for projection besides a distance. The data-driven Pswarm was compared to the common projection methods PCA, CCA, t-SNE, ESOM, NeRV and the MDS technique Sammon mapping. Using the new quality measure (Delaunay classification error) this work showed that the resulting two-dimensional projections of Pswarm are comparable to the state of the art projection methods like NeRV and ESOM. By using prior classifications, the Delaunay classification error allows for an unbiased evaluation of projection quality for both types of structures.

For the third module, the author expands the idea of previous works by using shortest paths between geodesic distances of the abstract U-matrix of projected points in the case of cluster analysis. The whole method is called Databionic swarm (DBS) and it outperforms the common clustering methods (k-means, PAM, single-linkage, spectral clustering, model based clustering and Ward) in terms of stability and plasticity on an artificial benchmark system of data sets (FCPS). Contrary to other common clustering methods, the DBS finds no clusters if no natural clusters exist. The number of clusters can be estimated with the help of the topographic map.

On three different high dimensional and multivariate data sets (types of leukemia, world gross domestic product, Tetragonula bees), the already known insights can be reproduced. In two real world applications of hydrology and pain genes, the DBS retrieves meaningful clusters, which was confirmed by domain experts.

Through the modularization, DBS can be generalized to projection to projection-based clustering. The visualization by the generalized U-matrix and the DBS clustering can be applied to every projection method for both types of structures. Through the use of the topographic map, results of common clustering methods can be verified by the structures found by Pswarm or any other projection method. Additionally, 3D prints of the visualized structures of high dimensional data sets can be manufactured with common 3D printing techniques

# **1 Introduction**

We live in a time when information is cheaply available and saved as data nearly everywhere. The amount of generated data is growing exponentially. By the end of the year 2016 alone, 9000 exabytes of data will have been generated, equal to 9 trillion gigabytes or the capacity of 360 billion Blu-ray Discs [Schiele, 2016]. The goal of the interdisciplinary field of data science is to extract knowledge from these data with the help of statistics, machine learning or data mining. Unlike in physics, a data scientist hardly ever starts with a hypothesis; he also is not interested in the source of the data or how they were collected. The data must be mined to gain knowledge through the identification of consistent patterns, and this is usually a very trying task.

Among the various available methods of analyzing data, the focal point of this work is cluster analysis. In contrast to common approaches, the goal here is not merely to group similar information but also to explain why the grouping of information in a certain context is valid, nontrivial and useful. Only then will the clustering of data be helpful to a domain expert. Cluster analysis "is a discipline on the intersection of different fields and can be viewed from different angles, which may be sometimes confusing because different perspectives may contradict each other" [Mirkin, 2005, p. 33]. From the statistical perspective, some assumption regarding the underlying model is required, and data clusters are viewed as probability distributions whose properties can be estimated from the data themselves [Mirkin, 2005, pp. 33-34]. "A trouble with this approach is that in most cases clustering is applied to phenomena of which nothing is known" [Mirkin, 2005, p. 34]. Here, cluster analysis is regarded as the process of generating a classification based on empirical data in a situation in which clear theoretical concepts and definitions are absent and the patterns and laws governing the situation are unknown (see [Mirkin, 2005, p. 36]). The concept of every application (available as open-source code in the R language [R Development Core Team, 2008]) used throughout this thesis is based on this idea.

The goal of this work is to provide an open-source framework for cluster analysis that is founded on a swarm-based projection method and uses a human-understandable visualization approach based on a topographic map of high-dimensional data structures, with the option of 3D printing (see [Thrun et al., 2016a]). This framework should be sufficiently stable while remaining adaptive and exhibiting sufficient plasticity to permit the creation of clusters of various shapes. It should include only a very few non-sensitive parameters that can be visually deduced by a non-professional data miner without any need to understand the theory behind them.

To achieve this goal, expertise on various topics from various areas of research will be required. It is the author's experience that experts in different fields rarely share or exchange practical approaches, and almost nobody is interested in providing and willing to provide easily available and human-understandable solutions to domain experts.

Here, the main hope is to be able to provide reproducible cluster analysis solutions for nonprofessional data miners and to deliver human-understandable concepts of high-dimensional data structures that are simultaneously able to be processed by machines. In the context of the Databionic swarm (DBS) approach, the author attempts to build, use and explain connections

among various fields of research; to be precise, the author will illustrate connections between cluster analysis [Hennig et al., 2015; Jain/Dubes, 1988], the imitation of collective behavior [Beni/Wang, 1993; Bonabeau et al., 1999; Reynolds, 1987], the visualization of information [Venna et al., 2010] and its evaluation, machine learning applications [Herrmann/Ultsch, 2008c], game theory [Nash, 1951], symmetry considerations in physics [Feynman et al., 2007, pp. 147-153, 745] and emergence [Ultsch, 2007]. Undoubtedly, making connections between different schools of thought sometimes requires simplifications. For example, with regard to the collective behavior of bees, the fact that bees have a queen who influences their behavior remains unaddressed in this work. Such simplifications are necessary for analytical modeling and applications of cluster analysis.

Chapter 2 addresses most of the necessary definitions and lays the groundwork for all of the mathematical notation used throughout the thesis. The literature reviewed in chapter 3 shows how common clustering methods tend to implicitly assume the patterns or structures sought in data. The reviewed clustering methods are grouped based on their definitions of generalized neighborhoods.

Chapter 4 introduces and classifies common methods of projecting high-dimensional data into two dimensions. Such projections are necessary to cope with the pitfalls of higher dimensions (see, e.g., [Bouveyron/Brunet-Saumard, 2014, pp. 55-57; Verleysen et al., 2003]). Two- or three-dimensional projections will always result in errors; however, gaining a spatial understanding of more than three dimensions is typically an excessively complex task for humans.

Chapter 5 presents examples to depict the typical errors encountered and describes efforts to manage these errors by means of the U-matrix visualization approach [Ultsch, 2003a]. By contrast, chapter 6 demonstrates a more stringent mathematical approach based on quality measures (QMs) presented in the literature. The evaluation of 19 QMs yields a grouping of the QMs based on their implied characterization of structures of high-dimensional data using the definition of neighborhoods introduced in this thesis. Consequently, it is not possible to generalize any of the QMs. If it were possible, the corresponding optimization approaches would not imply any prior assumptions about the structures of high-dimensional data and, consequently, would outperform any other projection methods.

Chapter 7 discusses a nature-inspired and behavior-based system of data science with the goal of using emergence, instead of the optimization of an objective function, for data visualization and clustering.

Building on the insights gained in chapter 7, chapter 8 introduces the DBS concept. Because it relies on the self-organization of data and emergence, DBS does not imply any particular structure that is sought in data. In the context of the projection, visualization and clustering of artificial or high-dimensional data, chapters 10-12 compare DBS with various common methods and apply the DBS framework both to reproduce known insights and to gain new knowledge about various types of data, e.g., multivariate time series or genetic data.

Readers may skip certain chapters depending on their interests. However, the contents of some chapters are based on insights from previous chapters, as indicated by arrows in Figure 1.1, which outlines the organization of this work. Please note, that due to technical limitations the figures and equations are numbered chapter wise.

Figure 1.1: Dependency graph of the chapters. BBS: behavior based systems; QAV: Quality Assessments of Visualizations; DBS: Databionic swarm. The underlying concept of DBS is based on insights from chapters 3, 5 and 7 (orange). The evaluation of DBS is performed in three steps (green): general validation in chapter 10, the reproduction of known knowledge in chapter 11, and the generation of new knowledge, as validated by domain experts, in chapter 12.

License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made. **Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **2 Fundamentals**

The first section of this chapter familiarizes the reader with the definitions of the basic notation and terminology used in this thesis. Concepts of graph theory are introduced in the next section. They give rise to a new concept of neighborhoods, which is utilized in several chapters. The last section explains a possible approach to knowledge discovery, which is applied in chapters 11 and 12.

# **2.1 Basic Definitions**

# **Hilbert space**

Let be a vector space above a field Κ with the following properties for every pair of elements :Κ∈ߙ and ∋ ሻݖ ,ݕ ,ݔሺ

1.) 〈.,.〉: x → Κ is a non-degenerate symmetric bilinear form:

 0 〈ݔ ,ݔ〉: ∋ ݔ ∀ .a

$$\text{b. } \langle \mathbf{x}, \mathbf{y} \rangle\_{\mathcal{H}} = 0, \forall \text{ } \mathbf{y} \in \mathcal{H} = \succ \mathbf{x} = \mathbf{0}$$

Thus, is a Hilbert space (for further details, see [Bronstein et al., 2005, pp. 635-636; Nolting, 2001, p. 22]).

# **Bra-ket notation**

Bra-ket notation 〈. |.〉 is used in physics to describe functions or vectors in a Hilbert space when the coordinate system of the vectors is irrelevant. The left part is called the bra (ۦ.(| , and the right part is the ket (|. ۧ). This notation is used to describe physical states (it is also called Dirac notation, as described in [Dirac, 1981, pp. 15-22]; for a formal introduction, see [Nolting, 2001, pp. 147-148]).

# **Operator**

An operator ܣመ is an unambiguous mapping of each element |ߙۧ of the subset ܦఈ ⊆ to an element |ߚۧ ∋ ܹ ⊇ such that หߚۧ ൌ ܣመ หߙۧ ൌ |ܣመ ߙۧ, where ܦఈ is the definition range of ܣመ and the set of all |ߚۧ is the domain of ܣመ, as defined in [Nolting, 2001, p. 153]; see also [Bronstein et al., 2005, pp. 49,639-640]. An "operator is considered to be completely defined when a result of its application to every ket vector [|ߙۧ [is given" [Dirac, 1981, p. 23].

# **Observation**

An observation *f* is a set of measured values for the properties of a phenomenon. It is described in the bra-ket notation as the change from one physical state ݕۦ |to another physical state |ݔۧ that results from the measurement of the operator ݂መ, as denoted by ݂ ൌ 〈ݕ|݂መ|ݔ) 〈see [Feynman et al., 2006, pp. 145, 147]). Such an observation *f* is a measurement of a physical process.

# **Feature**

Each individually measurable property *r* of a phenomenon being observed can be mapped to an operator ݎ̂ that can be applied to a physical state |ݔۧ] Stöcker et al., 2007, p. 744]. Such an individually measurable property is called a **feature**, **attribute** or **observable**. Here, an approximately continuous distribution of values in the vector space Թௗ is additionally assumed for a **variable** (see the definition of the **distribution of a variable**)**.**

# **Data**

A batch of data is defined as a matrix 〈݅หܣመ ห݆〉 ൌ ܣ, in which **facts1** about a physical state are summarized based on observations of the form 〈ݕหܣመ หݔ 〈ൌ ∑ 〈ݕ|݅〉〈݅หܣመ ห݆〉〈݆|ݔ 〈of a phenomenon in a Hilbert space, where ۦ݅|, ۦ݆|,| ݅ۧ and |݆ۧ are the basic states relevant to the phenomenon (for further discussion, see [Feynman et al., 2006, pp. 147-150]).

## **Distribution of a variable**

A formal distribution ݂݀ is defined as the probability density of a feature ݎ:

݂݀ሺݎሻ ൌ lim → 〈௫,| ௫〉 √ሺሻ [Nolting, 2001, p. 150]. If the feature *r* is continuous, then it is called a **variable** ݖ ∋ Թௗ, and *df* is called its probability density function (**pdf**) (see [Goodfellow et al., 2016, p. 58]). Here, when it describes how the relative probability of a variable ݖ takes on a given value, such a distribution is a pdf that is assumed to be normalized as follows [Walck, 2007, p. 15]: ݂݀ሺݖሻ݀ݖ ൌ 1 <sup>ஶ</sup> ିஶ .

"Statisticians often use the distribution function or as physicists more often call it the cumulative function which is defined as ݂ܿ݀ሺݖሻ ൌ ݂݀ሺݖሻ݀ݖ <sup>௫</sup> ିஶ " [Walck, 2007, p. 15].

If not elaborated further, here, the distribution of a variable *z* is regarded as an approximation of its pdf; for further details, see, for example, [Bock, 1974, p. 250; G. Ritter, 2014, p. 275 ff], and for types of pdfs, see [Walck, 2007].

# **Dirac delta function**

The Dirac delta function ߜ is a function with the following properties [Jackson, 1999, p. 31]: 1.) ߜሺݖെܽሻ ൌ 0 iff ݖ ് ܽ 2.) ߜሺݖെܽሻ ൌ ൜1, if ݖ ൌ ܽ lies in the integration area under the curve 0, otherwise

# **Density of data**

Let *dn* be the number of observations in an **elementary volume** (see [Bronstein et al., 2005, p. 491]) ݀ௗ෨ ݒ ൌ ݀ݒଵ ∗ ݀ݒଶ ∗. . ݀ݒௗ෨ ൌ ݀ݒ <sup>Ԧ</sup>of the Hilbert space Թௗ෨ (henceforth, Թௗ); then, the density of the data is defined as ߩሺݒ Ԧሻ ൌ ௗ ௗ௩ሬԦ , where ߩ: Թௗ → Թ is the density field function. Here, ߩ is subject to the condition that N is the number of data points defined by

ܰ ൌ ߩሺݒԦሻ <sup>Թ</sup> <sup>d</sup>ݒ Ԧൌ <sup>∑</sup> ߜሺݒԦെݒԦሻdݒ Ԧே <sup>Թ</sup> ୀଵ , in analogy to [Jackson, 1999, p. 33], where ߜ is the Dirac delta function and ߩሺݒԦሻ ൌ ∑ ݍߜሺݒԦെݒԦሻ ே ୀଵ is the charge density of point charges. Then, the **homogeneity** of the data is defined as

ܰ ൌ ߩሺݒԦሻdݒ ԦԹ ൌ ߩdݒ ԦԹ ൌ ߩ dݒ ԦԹ , where ߩ ൌ const.

 1 See [Fayyad et al., 1996, p. 6].

#### **Pattern**

A "[p]attern is an expression *E* in a language *L* describing facts [*F*] in a subset ܨா of *F*. *E* is called a pattern if it is simpler than the enumeration of all facts in ܨா" [Fayyad et al., 1996, p. 7]. Here, the expression *E* is "simpler" if it describes a group of similar (see the definitions of **metric space** and **distance** below) or homogeneous observations.

In graph theory, a pattern may be described by a neighborhood *H* (see the graph theory section for details). If the observations are not directly comprehensible, such a pattern is called a *hidden pattern*.

#### **Discontinuity in data**

A set of data can exhibit discontinuity if

$$\int\_{\mathbb{R}^d} \rho(\vec{v}) \mathrm{d}\vec{v} \neq \rho\_0 \int\_{\mathbb{R}^d} \mathrm{d}\vec{v},$$

which means that the density of data ߩ depends on its location ݒ Ԧin the Hilbert space Թௗ; Discontinuities can occur when interruptions or distortions exist in the homogeneity of the data, or in the continuity of the distribution of the data, in Թௗ. Thus, there are elementary volumes dݒ Ԧ with high density and elementary volumes dݒ Ԧwith low density or even empty elementary volumes. In the one-dimensional case, such a discontinuity can be mathematically defined as an essential or jump discontinuity. In two or three dimensions, a discontinuity may manifest as a spatial separation (see, e.g., Figure 2.1 or chapter 5 and 9, the Hepta data set).

In a higher-dimensional case, a discontinuity represents a change in the characteristics of facts, resulting in multiple patterns (see, for example, the leukemia data set, chapter 3, Figure 3.7 and chapter 9).

Figure 2.1: Spatial separation of data, after [Handl et al., 2005].

#### **Metric space and distance**

Let a metric space be represented by an ordered pair *(M, d)*, where *M* is an arbitrary set and *d* is a metric on *M*, i.e., a function

$$\mathsf{d}: \mathsf{M} \times \mathsf{M} \to \mathbb{R}$$

such that for any *l*, *j*, *m* ∈ *M* ,

$$\begin{aligned} d(l,j) &= d(j,l) \\ d(l,j) &\ge 0 \\ d(l,j) &= 0, \text{iff } l = j \end{aligned}$$

and the triangle inequality is satisfied as follows:

$$d(l, f) + d(j, m) \ge d(l, m)$$

Then, the metric *d* is also called a **distance** (see [Bronstein et al., 2005, pp. 624-625]). By contrast, for a **dissimilarity**, denoted by ݀መ, the triangle inequality may not apply ([Bock, 1974, pp. 25-26]. The distance between two **similar** points ݈, ݆ ∈ ܯ is small, whereas that between two **dissimilar** points ݈, ݆ ∈ ܯ is large. Transformations exist between a dissimilarity ݀መ and a distance *d* (e.g., [Bock, 1974, pp. 77-79]).

If the distance is defined in an output space O, it is denoted by *d(l, j)*, whereas a distance defined in an input space *I* is denoted by *D(l, j).* An example of a metric space is a Hilbert space that is a real-numbered vector space Թୢ of *d* dimensions. If the distances in a space are defined as Euclidean distances, then the corresponding space is called a Euclidean space.

#### **Data set**

A data set consists of a finite set of observations ݂∈F⊂ୢ෩ of ݀ሚ observed features.

In this work, observations ݂ are assumed to be vectors *l* in a metric space *M*, and features are assumed to be variables, if not stated otherwise.

#### **Input space**

An input space ܫ⊃Թௗ is the d-dimensional space consisting of ݀≤ d෨ variables in a data set that have been selected for a given task and contains *n* data points: ܫ ൌ ሼ݈ଵ,…,݈, ݊ ∈ Գሽ. The properties of an input space are as follows (see [Lee/Verleysen, 2007, p. 243]):

#### **Data point**

A data point ݈∈ܫ is a numeric vector consisting of one observation for each of the *d* variables in the input space, where a vector is an array of numbers arranged in a specific order defined with respect to the d variables.

<sup>2</sup> Note that, in general, the number of data points has greatly increased over time [Goodfellow et al., 2016, p. 21 , Fig. 1.8] and therefore the precise number may change with time

#### **Object**

When the data of interest are a set of facts F consisting of numerical, ordinal or nominal scaled entries, each fact ݂ ∈ ܨ, such that f ∉ Թୢ, is called an object or **case**.

An object can be regarded as a generalization of a data point. If an object can be interpreted (has a meaning within itself), then it contains **information** ([Ultsch, 2016c]; see also [Ultsch, 1994, p. 2]).

#### **Output space**

An output space ܱ⊂Թ is the m-dimensional space such that *m<d* in which, for each point ݆∈ܱ, a mapping to a data point ݈ of the input space I⊂ Թௗ exists.

#### **Machine learning**

The field of machine learning concerns computer programs that can imitate learning behavior [Natarajan, 2014] (see also [Goodfellow et al., 2016, p. 99]). Machine learning comes in two general forms<sup>3</sup> (see [Murphy, 2012, p. 2]). *Unsupervised learning* refers to the task of finding patterns in unlabeled data. Since the data are unlabeled, no reward function exists that can be used to evaluate potential results. If the data set is labeled, then *supervised learning* is possible. A typical supervised learning task is classification or regression. A typical unsupervised learning task is cluster analysis.

#### **Label**

A label is a tag ݃ ∈ ሼ1, . . . , ݇ሽ ⊂ Գ attached to an object ݂∈ܨ that identifies the object via a mapping ݂: ሼ1, . . . , ݇ሽ → ܨ. The labels of such a set of objects range from *one* to *k* [Hennig et al., 2015, p. 2], where k is the number of groups of objects. Here, it is assumed that a label exists for every object.

#### **Classification**

A classification ܥ ൌ ሼܩଵ, ܩଶ,…ሽ is a system of subsets [Bock, 1974, p. 22] such that ܥ⊃ୢ෩ . A subset ܩ ൌ ሼ݈ଵ, ,…݈ሽ݅ ∈ Գ, , is a set of k observations. In an exclusive classification, the subsets are disjunct, denoted by ܩଵ ∩ ܩଶ ൌ ∅; in a non-exclusive classification, elements that overlap between two subsets may exist, denoted by ܩ ∩ ܩ ്∅ . However, overlapping classification is not considered here (for various types of classification, see Figure 2.2 or [Hennig et al., 2015, p. 45]). Supervised and unsupervised classifications are defined as in the context of machine learning.

<sup>3</sup> Reinforcement learning is not considered in this context; semi-supervised learning (e.g. active learning) uses labeled data as well as unlabeled data.

Figure 2.2: Tree of classification types, after [Jain/Dubes, 1988, p. 56]. This work concentrates on unsupervised classification (see unsupervised machine learning).

#### **Classifier**

A classifier is an algorithm that constructs a function ܥ݈ݏ: ܨ → ሼ1, … , ݇ሽ ⊂ Գ that maps objects ݂∈ܨ to class labels ݃ ∈ Գ.

In terms of understandability, a distinction can be drawn between symbolic and sub-symbolic classifiers [Ultsch/Korus, 1993]. Symbolic classifiers are able to acquire knowledge (for a detailed description, see the last section of this chapter). By contrast, sub-symbolic classifiers (e.g., KNN classifiers) are only able to integrate knowledge [Ultsch, 1994], because a characteristic property of a sub-symbolic representation of data is that a single object alone does not contain information (see [Ultsch, 1994, p. 2]).

#### **Projected point**

A projected point ݆ሺݔଵ,..,ݔሻ ൌ ଔԦ is a vector of *m* scalars ݔ in the output space ܱ ⊂ Թ, where a vector is an array of numbers arranged in a specific order such that each individual number can be identified by its index.

#### **Projection**

Let ݆∈ܫ denote data points in the input space *I*⊂ Թௗ, and let ݈∈ܱ denote projected points in the output space O⊂ Թ. Then, a mapping proj: I → O, j ↦ l is called a projection iff ݉ ൌ ݀. ≫ ݉ ∧ ݐݏ݊ܿ

Note that unlike for a projection method, for a manifold learning method, the dimensionality of the output space ݉ depends on the data set (see, e.g., [Lee/Verleysen, 2007, pp. 14-15]).

#### **2.2 Concepts of Graph Theory Applied to Patterns**

This section uses graph theory to describe patterns found in data.

#### **Graph**

"A graph [Γሿ is a pair [Γ ൌ ሺV, Eሻ] consisting of a finite set ്ܸ∅ and a set E of two-element subsets of V. The elements of V are called vertices. An element e = (a, b) of E is called an edge with end vertices a and b. […] [In such a case,] *a* and *b* are adjacent or neighbors of each other" [Jungnickel, 2013, p. 2].

A graph ߁ is called undirected if, for every edge ݁ሺܽ, ܾሻ in E, the edge ݁ሺܾ, ܽሻ is also in E. A graph is called a weighted graph if a number (weight) is assigned to each edge.

#### **Directed graph**

A "directed graph or, for short, a *digraph* is a pair ߁ ൌ ሺܸ, ܧሻ consisting of a finite set *V* and a set *E* of ordered pairs *(a, b),* where ܽ ് ܾ are elements of V" [Jungnickel, 2013, pp. 25-26].

#### **Direct adjacency**

Let Γ be a graph, and let j be a point in a metric space M; then,

 ൟܧ ∋ ሻݒ ,ݒሺ݁ ∧ ∃ ܸ ∈ ݒ |ܯ ∋ ݈൛ ൌ ሻܯ ,߁ ݆,ሺ

is the set of points that are directly adjacent to j. The direct adjacency is defined by the specified graph.

#### **Adjacency matrix**

A digraph Γ with a vertex set ሼ1, . . . , ݊ሽ is specified by an ݊ ൈ ݊ matrix ܣ ൌ ሺܽ ሻ, where ܽ ൌ 1 if and only if ሺ݅, ݆ሻ is an edge of Γ, and ܽ ൌ 0 otherwise. *A* is called the adjacency matrix of Γ [Jungnickel, 2013, p. 40].

#### **Path**

Let ሺ݁ଵ,...,݁ሻ be a sequence of edges in a graph ߁. If there exist vertices ݒ,...,ݒ such that ݁ ൌ ݒିଵݒ for ݅ ൌ 1, . . . , ݊, then the sequence is called a walk; if ݒ ൌ ݒ, one speaks of a closed walk (Figure 2.3). A walk for which the ݁ are distinct is called a trail (Figure 2.3), and a closed walk with distinct edges is a closed trail. If, in addition, the ݒ are distinct, then the trail is a path [Jungnickel, 2013, p. 5].

Figure 2.3: Examples of trails, walks and paths [Jungnickel, 2013, p. 6 Fig. 1.5]: (a, b, c, v, b, c) is a walk but not a trail, and (a, b, c, v, b, u) is a trail but not a path [Jungnickel, 2013, p. 5].

### **Connected Graph**

Two vertices *a* and *b* of a graph *Γ* are called connected vertices if a walk exists with start vertex *a* and end vertex *b*. If all pairs of vertices of *Γ* are connected, then *Γ* itself is called a connected graph. For any vertex a, we consider *a* to be a trivial walk of length 0, such that any vertex is connected with itself. Thus, connectedness is an equivalence relation on the vertex set of *Γ*. The equivalence classes of this relation are called the connected components of *Γ*. Thus, *Γ* is connected if and only if its vertex set *V* is its unique connected component [Jungnickel, 2013, p. 6].

#### **Lattice**

A connected graph Γ with a particular well-defined two-dimensional tiling (tessellation) is defined as a lattice. A ݊ݔ݉ lattice has n vertices on the x-axis and m vertices on the y- axis. If the tiling is rectangular (every vertex has exactly four perpendicular edges) it will be called a **lattice**  (tiling) in this work, if the tiling is hexagonal (every vertex has exactly three edges) this will be called a **grid** (tiling) in this work.

#### **Shortest path**

For a connected graph Γ, there exists a distance *D(a, b*) between two vertices a and b that can be defined as the shortest path between these vertices [Jungnickel, 2013, pp. 65-66] as follows: For each path ܲ ൌ ሺ݁ଵ,…,݁ሻ, let the length of P be ሺܲሻ:ൌሺ݁ଵሻ ⋯ ሺ݁ሻ; then, the distance between two vertices a and b in *(Γ, p)* is defined by

ܩሺܽ, ܾ, ߁ሻ ൌ ൜ ∞, if ܾ is not accessible from ܽ minሼሺܲሻ: ܲ is a path from ܽ to ܾ ݅݊ ߁ሽ, otherwise

Let the vertices be denoted by points ݈, ݆ ∈ ܯ in the metric space *M*; then, G(l, j, ߁ (is the notation if the points ݈ and ݆ lie in the input space *I*, and gሺl, j, ߁ (is the notation if they lie in the output space *O*.

Note that ݀ሺܽ, ܽሻ ൌ 0 always holds because an empty sum is considered to have a value of 0, as usual. If no explicit length function is given, then the shortest paths and distances in a graph are defined using a length function that assigns a length of ሺ݁ሻ ൌ 1 to each edge *e* [Jungnickel, 2013, p. 66]. An algorithm for calculating the shortest paths in a graph is described in [Jungnickel, 2013, pp. 83-87]. The authors Lee and Verleyson have claimed that graph distances outperform the traditional Euclidean metric in terms of dimensionality reduction [Lee/Verleysen, 2007, p. 227].

#### **Acyclic graph**

Let ሺܯ,≽ ሻ be a partially ordered set (a poset, for short), which consists of the set *M* together with a reflexive, antisymmetric and transitive relation ≼, and let *M* correspond to a digraph ߁ with the vertex set *M* and with edges defined by pairs *(a, b)* such that ܽ ≺ ܾ; then, because of the transitive property, ߁ is acyclic [Jungnickel, 2013, p. 49].

#### **Tree**

A tree is a graph ߁ that satisfies the following three conditions [Jungnickel, 2013, pp. 7-8]:

The vertices in a tree are often called nodes. If ሺܽ, ܾሻ is an edge in a tree, then *a* is called the parent of *b*, and *b* is a child of *a*. If a path exists from *a* to *b* (്ܾܽ), then a is a proper ancestor of *b* and b is a proper descendant of *a* [Safavian/ Landgrebe, 1990, p. 2]. If a node has no descendant, it is called a leaf; if a node has no ancestor, it is called a root.

#### **Directed acyclic graph (DAG)**

A DAG is a directed tree (see above) that contains no cycles and one vertex, defined as the root, into which no edges enter. There is a unique path from the root to every vertex [Safavian/Landgrebe, 1990, p. 3]. Every vertex has a descendant called a child, except for the leaf vertices, which do not.

#### **Decision tree**

Let <sup>ܩ</sup> be a subset of a classification ܥ ൌ ሼܩଵ,…,ܩ,...ሽ⊆ୢ෩ ; then, a decision tree is a tree with the following properties:

#### **Decision tree learning**

Decision tree learning refers to a type of supervised machine learning in which decision trees are used (see [Safavian/Landgrebe, 1990]).

#### **Binary tree**

A binary tree is an ordered tree such that [Safavian/Landgrebe, 1990, p. 3] (see also the definition of a DAG)

#### **Lemma 1**

Let ߁ ൌ ሺܸ, ܧሻ be a connected graph with a positive length function p. Then, *(V, D)* is a finite metric space, where the distance function is defined as ܦൌܩሺܽ, ܾሻ [Jungnickel, 2013, p. 68].

#### **Proposition 1**

Any finite metric space can be represented by a pair (߁, ) (network) with a positive length function *p* [Jungnickel, 2013, p. 68].

#### **Ultrametric space**

Note that a metric space can be represented by a tree if and only if the following condition holds for any four vertices *x*, *y*, *z*, and *t* of the given metric space [Jungnickel, 2013, p. 69]:

ሻሻݖ ,ݕሺ݀ ሻݐ ,ݔሺ݀ ,ሻݐ ,ݕሺ݀ ሻݖ ,ݔሺ݀ሺݔ݉ܽ ሻݐ ,ݖሺ݀ ሻݕ ,ݔሺ݀

Changing the triangle inequality to this condition implies an ultrametric space.

#### *2.2.1 Patterns Defined as a Generalization of Neighbourhoods*

Here, it is argued that by using shortest paths and direct adjacency, the patterns that exist in data can be generalized to neighborhoods H of an extent *k*.

Let k∈Գ, *k>0*, let Γ be a connected graph, let ݆ be a point in a metric space ܯ, and let ܩሺ݆, ݈, ߁ሻ be the shortest path between ݆∈ܯ and an arbitrary point ݈∈ܯ ;then (1),

$$H\_f(k, \Gamma, M) = \{ l \in M \mid G(l, j, \Gamma) \le k \} \tag{1}$$

is the neighborhood set of the point *j* and *k* the neighborhood extent. The neighborhood *H* can define a pattern in the input space4 .

The easiest example is a neighborhood defined by distances in a Euclidean graph. In the context of graph theory, a Euclidean graph is an undirected weighted graph of the highest order with respect to all other graphs discussed here, because every vertex is connected to every other vertex. Note that the weights of the vertices in a Euclidean graph need not necessarily be defined by the Euclidean metric. Another representation of a neighborhood H is a Delaunay graph ࣞሺܸ, ܧሻ, which is a subgraph of a Euclidean graph. A Delaunay graph ࣞሺܸ, ܧሻ is based on Voronoi cells [Toussaint, 1980]. Each cell is assigned to one data point, and the size of a cell is characterized in terms of the nearest data points surrounding the point assigned to that cell. Within the borders of one Voronoi cell, there is no position that is nearer to any outer data point than to the data point within the cell. Thus, a neighborhood of data points is defined in terms of direct links between borders of Voronoi cells that induce an edge E in the corresponding Delaunay graph [Delaunay, 1934]. In short, a Delaunay graph represents a graph for a neighborhood ܪሺ1, ࣞ, ܯሻ. A neighborhood H can also be represented by a Gabriel graph ܩሺܸ, ܧሻ [Gabriel/Sokal, 1969], which is a subgraph of a Delaunay graph ࣞሺܸ, ܧሻ in which two points are connected if the line segment between the two points is the diameter of a closed disc that contains no other points within it (empty ball condition). A Gabriel graph represents a graph for a neighborhood ܪሺ1, ܩ, ܯሻ. Another case that is often considered is that of a neighborhood ܪ*)knn,* ܭ*, M),* where the number of nearest neighbors of a point *j* is defined by the number of vertices connected to this point in the *K*-nearest-neighbor graph (KNN graph), e.g., [Brito et al., 1997]. Here, we will use the shorter notation *H(knn, M).*

Figure 2.4: Four points and their Voronoi cells: *D(l, k)>D(l, m)* illustrate the different types of neighborhoods: unidirectional versus direction-based.

<sup>4</sup> Such neighborhoods *H* will prove useful for various evaluation steps, as summarized in Fig. 2.5.

Neighborhoods of points can be divided into two types, namely, *unidirectional* and *directionbased* neighborhoods. Consider the four points shown in Figure 2.4. The points *l*, *k*, *j*, and *l* are in the same neighborhood ܪሺ1, ࣞ, ܯ (in the corresponding Delaunay graph, but the points *l* and *m* are never neighbors in this graph, even if the distance *D(l, m)* is smaller than *D(l, k).* Thus, in this neighborhood definition, the direction information is more important than the real arrangement of the points in space as characterized by the distances *D*.

However, if a neighborhood is defined in terms of a KNN graph, then the points *l* and *m* could be in the same neighborhood ܪሺ݇݊݊, ܭ, ܯ(, and the points *l* and *k* could be in different neighborhoods, depending on the value of ݇݊݊ and on the ranking of the distances between these points. Therefore, this type of neighborhood is called unidirectional. In other words, it can be said that the points *l*, *j*, and *m* are more *dense* with respect to each other than they are with respect to *k*. Thus, unidirectional neighborhoods defined in terms of KNN graphs or unit disk graphs [Clark et al., 1990] can be used to define neighborhoods based on density.

#### **2.3 Overview of Knowledge Discovery**

 *"The term knowledge discovery in databases […] was coined in 1989 to refer to the general process of finding knowledge in data and to emphasize the 'high-level´ application of particular data mining methods" [Fayyad et al., 1996, p. 3].*

In 1996, Fayyad et al. used this term in his introduction to "From Data Mining to Knowledge Discovery" as follows:

*"Knowledge discovery in databases is the non-trivial process of identifying valid, novel, potentially useful, and ultimately understandable patterns in data" [Fayyad et al., 1996, p. 6].*

Dropping the suffix *in databases*, the term *knowledge discovery* was extensively discussed in [Mörchen, 2006, pp. 6-7]. According to the definition used in that work, *knowledge discovery* is "data mining with the goal of finding knowledge, i.e., novel useful, interesting, understandable, and automatically interpretable patterns" [Mörchen, 2006, p. 7]. The definition of *data mining* as given in [Mörchen, 2006, p. 7] is

*"The process of finding hidden information or structure in a data […] [set.] This includes extraction, selection, preprocessing, and transformation of features describing different aspects of the data".*

The following overview in Figure 2.5 presents a possible approach to knowledge discovery, as applied in chapters 11 and 12. It is not claimed here that this view is the only approach available in this research field. The remainder of this chapter will describe the various tasks involved in knowledge discovery which are shown in Figure 2.5.

Figure 2.5: The step-wise process of knowledge discovery, as inspired by [Fayyad et al., 1996, p. 10; Ultsch, 2000b]. The systematic process may contain loops between any steps [Behnisch/Ultsch, 2015, p. 52]. This work focuses on Clustering analysis which will be separately discussed in the next chapter, but in general applying Machine learning algorithms would be the 4th step.

#### *2.3.1 Feature Selection*

In the first step, the "features must be properly selected so as to encode as much information as possible concerning the task of interest. […] minimum information redundancy among the features is a major goal" [Theodoridis/Koutroumbas, 2009, pp. 596-597] (see also [Lee/Verleysen, 2007, p. 230]). Redundancy refers to a case in which certain features of a data set are not independent of each other [Lee/Verleysen, 2007, pp. 1-2]. For example, if the two variables ݈ and ݆ are correlated, then ܦሺ݈, ݆ሻ ൌ ඥ∑ ݈ െ ݆ is no longer a Euclidean distance [Cormack, 1971, p. 326].

#### *2.3.2 Preprocessing*

 *"Preprocessing the data to be mined is utterly important for a successful outcome of the analysis. If the data is not cleansed and normalized, there is a high danger of getting spurious and meaningless results. Cleansing includes the removal of outliers, i.e., data objects with extreme values, replacement of missing values, or the removal of erroneous corresponding data sets" [Mörchen, 2006, pp. 7-8].*

Sometimes, this first step is already referred to as feature extraction [Bishop, 2006, p. 2]. Many data mining methods rely on the concept of (dis-)similarity between pieces of information encoded in data. For example, for Euclidean distances, "normalization of the data needs to be considered to avoid undesired emphasis of features with large ranges and variances" [Mörchen, 2006, p. 8] (see also [Jain/Dubes, 1988, p. 38]). This process of creating such "syntetic" data features that retain the most important information of a pattern in question is here called feature extraction (consistent with [Mirkin, 2005, p. 208]).

#### *2.3.3 Feature Extraction*

The first step of feature extraction is to determine the distribution of each individual variable.

 *<sup>&</sup>quot;Important tools for this inspection are the quantile-quantile plot (QQ-plot) and kernel estimators for the probability density function (pdf). Here we use the PDE method for pdf estimation [Ultsch, 2003b] as it is specially designed to uncover subsets in the variables" [Behnisch/Ultsch, 2015, p. 54].*

A QQ-plot makes it possible to compare the given distribution of a variable to standard distributions. Additionally, box-whisker diagrams (boxplots) may be used to visualize the quartiles of a variable.

## *<IP_ADDRESS> Transformations*

*"Real valued data often comes from domains where variables have greatly varying variances because of different scales. Variables with large variances are likely to dominate the obtained distance structure, e.g. when using Minkowski metrics. To overcome this problem, each variable is linearly transformed (standardized) such that the estimated variance is the same on all variables. The Z-score scheme transforms a variable's values* ݔ ← ሺݔ െ ݉ሻ/ߪ *with mean m and standard deviation σ"* [Herrmann, 2011, p. 28].

If a variable can be non-linearly transformed to a normal distribution, the Box-Cox algorithm (see [Asar et al., 2014]) is often used to estimate the factor of the transformation. With an approximation of the factor obtained from the ladder of powers [Tukey, 1977], an "understandable" transformation, e.g., "log" or "sqrt," can be applied that is as near as possible to the factor of the Box-Cox algorithm. "These allow for hypotheses on why the distribution is shaped in a particular way" [Behnisch/Ultsch, 2015, p. 56].

For non-normally distributed variables (e.g., a variable with a multimodal distribution), a meaningful variance ߪଶ may be difficult to estimate. "Instead, a (robust) min/max-standardization transforms a variable's values ݔ ← ௫ିሺ௫ሻ ௫ሺ௫ሻିሺ௫ሻ with robust estimates ݉݅݊ሺݔሻ, ݉ܽݔሺݔሻ for minimum and maximum values. There is empirical evidence by Milligan and Cooper [Milligan/Cooper, 1988] that min/max standardization is to be preferred over Z-score, especially if variances of underlying distributions is [sic] hard to estimate" [Herrmann, 2011, p. 28]. In this context, ݉ܽݔሺݔሻ and ݉݅݊ሺݔሻ are estimated as the 95th and 5th percentiles, respectively, of the distribution [Herrmann, 2011, p. 127].

### *<IP_ADDRESS> Dimensionality Reduction*

A common approach to feature extraction is dimensionality reduction (DR). To cope with the "curse of high dimensionality" (for further details, see [Verleysen et al., 2003]), dimensionality reduction reduces an input space I⊂ Թௗ to an output space O⊂ Թ such that ݉൏݀ [Lee/Verleysen, 2007].

*"All difficulties that occur when dealing with high-dimensional data are often referred to as the 'curse of dimensionality´. When data dimensionality grows, the good and well-known properties of the usual 2D or 3D Euclidean spaces make way for strange and annoying phenomena" [Lee/Verleysen, 2007, p. 3].*

The various phenomena related to this concept are explained in [Lee/Verleysen, 2007, pp. 4-9] (see also [Bellman, 1957]). A DR method is usually either a manifold learning method or a projection method. DR methods such as autoencoders [Hinton/Salakhutdinov, 2006], Isomap [Tenenbaum et al., 2000] or local linear embedding (LLE) [Roweis/Saul, 2000] that are designed to find a manifold5 that represents a given set of high-dimensional data6 are called *manifold learning* methods. Such methods are disregarded here because these manifolds usually have more than two dimensions. DR methods of the type known as projection methods are

 5 "A manifold is a connected region. Mathematically, it is a set manifold of points, associated with a neighborhood around each point. From any given point, the manifold locally appears to be a Euclidean space." [Goodfellow et al., 2016, p. 160] 6

Often described using the term *intrinsic dimension* (e.g., [Lee/Verleysen, 2007, pp. 18-24, 41, 47ff]).

separately introduced in chapter 4. There, the focus is placed on methods that attempt to visualize information by means of projections that are restricted to visualizing high-dimensional data in a two-dimensional space while preserving their structure (for details, see chapter 5). The quality of a projection critically depends on the concept of dissimilarity that is chosen to be applied to the input space *I*. This concept could be a definition based on either distance or local proximity. An index used to evaluate the quality of a projection is called a quality measure (QM), and 19 QMs are introduced in chapter 6.

#### *2.3.4 Cluster Analysis*

Many data mining methods rely on some concept of the dissimilarity between pieces of information encoded in the data of interest. These methods are used for cluster analysis, and common approaches will be described in the next chapter. Cluster analysis is the task of unsupervised classification that results in a clustering. Given a data set *I* that contains *n* data points, the objective of cluster analysis is to group the data points into *K* disjoint subsets of *I*, denoted by ܿଵ,…,ܿ [Hennig et al., 2015, p. 2]. "A clustering is […] the partition obtained" with

Κ ൌ ሼܿଵ,…ܿሽ. If a data point *l* belongs to a cluster ܿ, then it has the class label ݃∈Գ. In the literature, this process is often called hard clustering to distinguish it from methods such as fuzzy clustering, in which a fractional degree of membership is assigned to each ݈∈ܫ] Jain et al., 1999].

#### **Cluster**

No generally accepted definition of clusters exists in the literature [Hennig et al., 2015, p. 705]. When describing clusters, the term *pattern* is often used (e.g., [Theodoridis/Koutroumbas, 2009]).

Here, consistent with Bouveyron et al., it is assumed that a cluster is a group of similar objects [Bouveyron et al., 2012]. Chapter 3 will elaborate on this statement while presenting the definition of *natural* clusters.

#### **Intracluster Distance**

Let ܿ ⊂ ܫ be a cluster such that ∀ܿ ⊂ ܫ, where , ݍ ∋ ሼ1, … , ݇ሽ and ്ݍ, ܿ∩ ܿ ൌ ሼ ሽ; then, the distance ܫ݊ݎݐܽሺܿሻ ≔ ܦሺ݈, ݆ሻ between two data points ݆, ݈ ∈ ܿ, is called an intracluster distance.

#### **Intercluster Distance**

Let ܿ ⊂ ܫ and ܿ ⊂ ܫ be two clusters such that , ݍ ∋ ሼ1, … , ݇ሽ , ܿ ∩ܿ ൌ ሼ ሽ, and ്ݍ ; then, the distance ܫ݊ݐ݁ݎ൫ܿ, ܿ൯ ൌD(݆, ݈) between two data points ݆ and ݈ in the two clusters, ݆∈ܿ and ݈∈ܿ, is called an intercluster distance.

#### **Compact Structures**

Compact structures in a data set are mainly defined by distances *d* if discontinuity in data exist such that the intracluster distances are small and the intercluster distances are large. Note, that the distance distribution is often bimodal if the data structures are compact. This type of structures leads to natural clusters (see chapter 3).

#### **Connected Structures**

Connected structures in a data set are mainly defined by density ߩሺݒ Ԧሻ if discontinuity in data exist. If a connected graph Γ is chosen appropriately regarding the data set, these data structures are based on neighborhoods ܪ ሺ݇, Γ, ܯሻ. This type of structures leads to natural clusters (see chapter 3).

#### *2.3.5 An Approach to Knowledge Acquisition*

If, for a given data set, there exist labels defined by a clustering or a domain expert, the next step may be to determine what each cluster means [Behnisch/Ultsch, 2015, p. 65] or what kind of knowledge can be acquired from it7 .

*"Under knowledge we understand a symbolic representation of objects, facts and rules for an interpreter with symbol processing capability, e.g. a human8 . In particular, knowledge is communicable by word or writing" [Ultsch, 1994, p. 1] (see also [Ultsch, 1987, p. 22]).*

Knowledge has the properties of being valid, comprehensible, nontrivial, potentially innovative and useful in practice [Behnisch/ Ultsch, 2015, p. 52]. It can be stored in a knowledge base, which "is an organized collection of knowledge together with operations for accessing and manipulating knowledge" [Ultsch, 1987, p. 22]. One example of a representation of knowledge is a rule [Ultsch, 2016c], which is defined as a prescription regarding how to generate, interpret and manipulate facts [Ultsch, 1987, p. 22].

In the context of knowledge discovery, knowledge acquisition can be defined "as the encoding of knowledge into the formal representation scheme of a knowledge-based system [KBS]" [Ultsch, 1987, p. 23]); here, a KBS is defined as "a computer program that contains an explicit, formal representation of knowledge in a knowledge base and is capable of [drawing conclusions9 ]" [Ultsch, 1987, p. 23]. In another context, researchers may interview domain experts "to become educated about the domain and to elicit the required knowledge, in a process called knowledge acquisition" [Russell et al., 2003, p. 217]. In short, knowledge acquisition can be described as a process that leads to a formal representation of knowledge (see [Aikins, 1983]), for example, a process leading to the generation of rules required for a computer program, e.g., DENDRAL [Russell et al., 2003, p. 22] or MYCIN [Aikins, 1983]. One possible approach to knowledge acquisition is to use machine learning [Russell et al., 2003, p. 687]. With regard to understandability, the machine learning methods used for this purpose can be classified as either symbolic or sub-symbolic methods [Ultsch/Korus, 1993].

*"Sub-symbolic methods model the structure of data using many numerical parameters. They are usually aimed at prediction or classification. The output of sub-symbolic methods often depends on the values and interactions of most or all model parameters. They fail to explain the prediction or classification. There are certainly areas of data mining where it is sufficient to build such black-box models that can approximately reproduce a classification or predict future data. An important requirement for knowledge discovery is the interpretability of the results. In many domains the expert wants to know why a decision was made or what a […] pattern describes. Comprehensible descriptions of the models are crucial for success in this case"* [Mörchen, 2006, p. 120].

For the acquisition of knowledge through cluster analysis, symbolic methods are preferable, as described in chapters 11 and 12 (see also [Ultsch, 1994]). In chapter 12, decision tree learning

<sup>7</sup> In another context one would like to explain a prediction done by a machine learning algorithm.

<sup>8</sup> For humans 7±2 rules appear to be the optimum [Miller 1956].

<sup>9</sup> Formally defined as *inference* in [Ultsch, 1987, p. 22].

is used in a knowledge acquisition approach called Classification And Regression Tree (CART) analysis [Breiman et al., 1984]). This method relies on a binary tree in which the splitting criteria (decisions) for the vertices are expressed in terms of the Gini index (for further details, see [Safavian/Landgrebe, 1990, p. 15]).

"A class is described by a number of conditions" [Ultsch/Korus, 1993, p. 3] that lead to the generation of a subset ܩ ⊃ ܥ defined by a previously identified clustering. Additionally, for each class, a unique class label ݃∈Գ exists for all ∋ ܩ. Every observation ∋ܩ can be unambiguously described by one or more properties that are shared among all observations of ܩ. Here, the conclusion that an observation can be correctly assigned to a class ܩ is reached based on the conditions defining a path (rule) from the corresponding leaf to the root of the binary tree, and this conclusion is called the decision to place in ܩ. Therefore, the class ܩ has a semantic characterization because it is characterized by the rules governing the decision tree, which allow this class to be distinguished from other classes. Here, it is assumed that the last step in the evaluation of a clustering is to ask domain experts to validate the identified classes.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made. **Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International

# **3 Approaches to Cluster Analysis**

Many data mining methods rely on some concept of the similarity between pieces of information encoded in the data of interest. Various names have been applied to these clustering methods, depending largely on the field of application in data science. For example, in biology the term "numerical taxonomy" is used [Thorel et al., 1990], in psychology the term Q analysis is sometimes employed, market researchers often talk about "segmentation" [Arimond/Elfessi, 2001] and in the artificial intelligence literature, unsupervised pattern recognition is the favored label [Everitt et al., 2001, p. 4]. The corresponding methods can be either data-driven or needdriven. The latter, called also constraint clustering [Tung et al., 2001] aims at organizing the true structure to meet certain application requirements such as energy aware sensor networks, privacy preservation, and market segmentation [Ge et al., 2007, p. 320]. An overview of constrained clustering algorithms can be found in [Basu et al., 2008].

Here, however, the focus is placed on data-driven10 methods, in which patterns present in the data are used to identify homogeneous groups of objects [Arabie et al., 1996, p. 8 ff.]. Consequently, the term *cluster analysis* is used to refer to a step in the knowledge discovery process (chapter 2, Figure 2.5.). Let it be assumed that in Figure 3.1 (top left), the first data set (I) contains two variables11. The division of this homogeneous data set into different patterns would be called dissection [Everitt et al., 2001, p. 7]. By contrast, *natural clusters* do not require dissection; instead, they are clearly separated in the data [Duda et al., 2001, p. 539; Theodoridis/Koutroumbas, 2009, pp. 579, 600], as shown in the second data set (II) in Figure 3.1 (top right).

No generally accepted definition of clusters exists in the literature [Hennig et al., 2015, p. 705]. Additionally, Kleinberg showed for a set of three simple properties (scale-invariance, consistency and richness), that there is no clustering function12 satisfying all three [Kleinberg, 2003]. By concentrating on distance and density based *structures13*, this work restricts clusters to "natural" clusters (see section 2) and therefore omits the axiom of richness where all partitions should be achievable. Consequently, only natural clusters, in which objects are similar within clusters and dissimilar between clusters [Bouveyron et al., 2012], are considered here. For example, the distance distribution in the input space can be bimodal, indicating a distinction between the inter- versus intracluster distances: in data set I in Figure 3.1 (bottom left), no large intercluster distances exist and the distribution of the distances is unimodal, whereas in data set II in Figure 3.1 (bottom right), the distribution of the distances is bimodal because data set II contains two natural clusters with a large intercluster distance. Another example is the case in which the number of data points in one *elementary volume* (݀ݒ (Ԧof the input space is higher than that in another elementary volume ݀ݒԦ, which can be estimated using a nonparametric technique for density estimation (e.g., kernel density estimation). In a third example, local proximities can be defined as structures based on neighborhoods ܪሺ݇, ߁,ܯሻ (see chapter 2.2.1).

M. C. Thrun, *Projection-Based Clustering through Self-Organization and Swarm Intelligence*, https://doi.org/10.1007/978-3-658-20540-9\_3

<sup>10</sup> The progress in an "algorithmic activity" is enforced by data w.r.t. patterns (as opposite to intuition or personal experience, e.g. through the setting of parameters). 11 In fact, this figure shows a CCA projection of the leukemia data set (see chapter 9). 12 "[A]ny function f that takes a set S of n points with pairwise distances between them, and returns a partition of

S" [Kleinberg, 2003, p 2]. 13 They can be described as patterns identified based on discontinuity.

Figure 3.1: Data set I is an approximately homogeneous data set with patterns that form no natural clusters (left, top). The distance distribution in this case is not bimodal (left, bottom). Data set II contains two natural clusters with a large intercluster distance (right, top). The distance distribution is bimodal here (right, bottom). See Figure 12.2 or supplement B for a high-dimensional example. Distance distributions was generated using the AdaptGauss CRAN package [Thrun/Ultsch, 2015; Ultsch et al., 2015].

#### **3.1 Common Clustering Methods**

Clustering methods can be broadly divided into two groups: hierarchical and partitional methods [Jain, 2010]. Partitional clustering methods simultaneously divide a set of data points into subsets. Because we are concentrating on *natural clusters,* overlapping clustering is not considered here. It should be remarked that the choice of the clustering algorithm to be used is more important than the choice of the distance calculation [Jain/Dubes, 1988, p. 140].

A prominent example of a partitional clustering method is the well-known *k-means* method of [MacQueen, 1967] (originally from [Steinhaus, 1956]). It proceeds as follows: Once the number of clusters has been chosen, a random initialization of cluster centers, called centroids, is performed in the input space. Then, the nearest data points to each centroid are assigned to that centroid. After the mapping of the data points, the centroids are moved such that the distances from the assigned points to their corresponding centroids are minimized. This process is performed repeatedly. Figure 3.2 illustrates four iterations of the process. In summary, k-means centroids are average points rather than individual data points. Details about the algorithm can be found in [Hennig et al., 2015, p. 68ff].

By contrast, the clustering method called partitioning around medoids (PAM), introduced in [L. Kaufman/Rousseeuw, 1990], minimizes the sum of the distances from the data points within a cluster to one chosen data point in the same cluster, called the medoid [Mirkin, 2005, p. 181]. In other words, the average distance between a medoid and a subset of data points in the same cluster is minimized. Aside from the change from centroids to medoids, the algorithm can be formulated analogously to k-means [Mirkin, 2005, p. 182].

Hierarchical clustering algorithms are based on the "representation of data as a hierarchy of clusters nested over set-theoretic inclusion" [Mirkin, 2005, p. 112]. In the agglomerative approach, such an algorithm begins with each data point in its own cluster and successively merges the most similar pairs of clusters to form a cluster hierarchy14.

A typical visual representation of this process is called a dendrogram (Figure 3.3). A dendrogram is a tree showing a hierarchical structure of distance-based connections between subsets of points. The similarity between points or groups of points depends on the algorithm. [Bock, 1974] demonstrated (see chapter 2 for details) that for every dendrogram, an ultrametric space can be constructed in which the triangle inequality is redefined as

ܦሺ݈, ݆ሻ ݉ܽݔ ሺܦሺ݈,݉ሻ, ܦሺ݉, ݆ሻሻ.

Figure 3.2: Steps of iteration using the k-means algorithm. After a random initialization of three centroids the nearest data points are assigned to each centroid. Then the centroids are moved to minimalize the distances.

<sup>14</sup> The divisive approach is not considered here (see [Mirkin, 2005, p. 113 ff] for details).

Figure 3.3: Dendrogram of the Hepta data set based on the Ward algorithm. Large changes in fusion levels of the ultrametric portion of the Euclidean distance in the Ward algorithm (y-axis) indicate the best cut. Seven clusters are indicated by red boxes at the y-axis value of 10. If only small changes in the fusion levels exist, it indicates that the algorithm is not able to find a cluster structure.

One of the most common hierarchical clustering algorithms is called *single linkage* (SL) [Florek et al., 1951; Sokal/Sneath, 1963], in which the clustering process is agglomerative [Jain et al., 1999]. In SL, the similarity between two subsets of data points is defined as the minimum distance between data points in these subsets [Duda et al., 2001, p. 553].

Let ܦ ෩be the distance between two clusters ܿଵ ⊂ ܫ and ܿଶ ⊂ ܫ, and let ܦሺ݈, ݆ሻ be the distance between two data points in the input space I; then, SL is defined based on (see [Hennig et al., 2015, p. 9]) ܦ෩ሺܿଵ, ܿଶሻ ൌ min ∈ భ,∈మ ܦሺ݈, ݆ሻ.

In graph theory terminology, this process generates a tree [Duda et al., 2001, p. 553]. If it is allowed to continue until all subsets of points are linked, the result is a (minimal) spanning tree (MST) [Duda et al., 2001, pp. 553, 554; Jain/Dubes, 1988, p. 70]. Of all common algorithms developed before 1968, only SL satisfies all conditions of a "theoretically valid" clustering (see [Jardine/Sibson, 1968] for details).

Another hierarchical clustering algorithm that will be used here is called the *Ward* algorithm [Ward Jr, 1963]. In the Ward algorithm, the similarity between two subsets of points is based on an optimal value of an objective function, which commonly is the sum of squared errors (*SE)*.

Let ܿ ⊂ ܫ and ܿ ⊂ ܫ be two clusters such that ݎ, ݍ ∋ ሼ1, … , ݇ሽ and ܿ ∩ܿ ൌ ሼ ሽ for ݎ്ݍ, and let the data points in the clusters be denoted by ݆ ∈ ܿ and ݈ ∈ ܿ, with the cardinality of the sets being ݇ ൌ หܿห and ൌ |ܿ| and with

 $m(c\_q) = \frac{1}{k} \sum\_{l=1}^{k} j\_l$  and  $m(c\_r) = \frac{1}{p} \sum\_{l=1}^{p} l\_l$ .

then, the *SE* is defined as (see [Theodoridis/Koutroumbas, 2009, pp. 661-663])

$$SE = \frac{k \ast p}{k + p} \sum\_{l=1}^{n} \left( m(c\_k) - m(c\_p) \right)^2$$

In Figure 3.3, the ultrametric property of the Ward algorithm is represented in a dendrogram (for further details, see [Duda et al., 2001, p. 557; Everitt et al., 2001, p. 68ff; Jain/Dubes, 1988]). If the values on the y axis "for the levels are roughly evenly distributed throughout the range of possible values, then there is no principled argument that any particular number of clusters is better or more natural than another" [Duda et al., 2001, p. 551]. "Large changes in fusion levels are taken to indicate the best cut" [Everitt et al., 2001, p. 76]. The cut depicted in Figure 3.3 generates a clustering consisting of seven clusters of roughly equal size. The next clustering method used in this work is called spectral clustering.

*"[It] is a class of graph-based techniques that unravel the structure properties of a graph using information conveyed by the spectral decomposition [eigendecomposition [see [Goodfellow et al., 2016, pp. 42-44]]] of an associated [Laplacian] matrix. The elements of this matrix code the underlying similarities among nodes [data points] of the graph" [Theodoridis/Koutroumbas, 2009, p. 772].*

 *"The K principal eigenvectors of the Laplacian matrix provide a mapping of the objects into K dimensions. To obtain clusters, the resulting K-dimensional vectors are clustered by standard methods, usually K-means. There are various interpretations of this. […]. For these [Euclidean] data, spectral clustering acts as a remarkably robust linkage method." [Hennig et al., 2015, p. 10].*

There is a close resemblance between spectral clustering and manifold learning methods [Theodoridis/Koutroumbas, 2009, p. 779]. Here, the clustering algorithm of [Ng et al., 2002] is used to take advantage of the open-source implementation of this method that is available in the R language [R Development Core Team, 2008].

 "Clustering via mixtures of parametric probability models is sometimes in the literature referred to as 'model-based clustering´" [Hennig et al., 2015, p. 10]. With the clustering algorithm of [Fraley/Raftery, 2006] in mind, here, this clustering method is called the *mixture of Gaussians* (MoG) method. The MoG method uses the *expectation maximization* (EM) algorithm (for further details on the EM algorithm, see [Bishop, 2006]).

*The EM algorithm is "an algorithm of alternating maximization applied to the likelihood function for a mixture of distributions model. At each iteration, EM is performed according to the following steps: (1) Expectation: Given parameters of the mixture* ܲ *and individual density functions* ܽ*, find posterior probabilities for observations to belong to individual clusters* ݃ *[…]. (2) Maximization: given posterior probabilities* ݃*, find parameters* ܲ*,* ܽ *maximizing the likelihood function" [Mirkin, 2005, p. 178].*

The MoG method suffers "from the well-known curse of dimensionality [Bellman, 1957], which is mainly due to the fact that model-based clustering methods are over-parametrized in high-dimensional spaces" [Bouveyron/Brunet-Saumard, 2014, p. 53]. To solve this problem, "for model based clustering, variable selection can be tackled within a Bayesian framework" [Bouveyron et al., 2012]. In the case of the MoG clustering method, the optimal model can be calculated according to the Bayesian information criterion [Aho et al., 2014] for parameterized Gaussian mixtures that are EM initialized using hierarchical agglomeration [Fraley/Raftery, 2002, pp. 10-12].

*"In each hierarchical agglomeration, each stage of merging corresponds to a unique number of clusters, and a unique partition of data. A given partition can be transformed into indicator variables […] which can then be used as conditional probabilities in an M-step of EM for parameter estimation, initializing an EM iteration" [Fraley/Raftery, 2002, p. 11]. Here, the R package mclust is used [Fraley/Raftery, 2006].*

# **3.2 Structure of Natural Clusters**

*"Clusters can be of arbitrary shapes (structures) and sizes in a multidimensional pattern space. Each clustering criterion imposes a certain structure on the data, and if the data happen to conform to the requirements of a particular criterion, the true clusters are recovered. Only a small number of independent clustering criteria can be understood both mathematically and intuitively. Thus the hundreds of criterion functions proposed in the literature are related and the same criterion appears in several disguises"* [Jain/Dubes, 1988, p. 91].

This section analyzes common clustering algorithms from the perspective of structures, whereas in various other sources, the clustering criterion or objective function has been understood only intuitively. Here, it is argued that the main argument of Jain and Dubes has received overall consent from the clustering community: Different clustering methods tend to implicitly assume different structures of clusters [Duda et al., 2001, pp. 537, 542, 551; Everitt et al., 2001, pp. 61, 177; Handl et al., 2005; Theodoridis/Koutroumbas, 2009, pp. 862, 896; Ultsch/Lötsch, 2016].

# *3.2.1 Types of Structures Sought by Clustering Algorithms*

The argument of Handl et al. is partially adopted here, in which natural clusters are considered to exhibit two types of structures, called compact and connected structures [Handl et al., 2005], as depicted in Figure 3.4. Clusters with compact structures show small variations in their intracluster distances; connected structures are based on the idea of neighborhoods of data points [Handl et al., 2005]. Here, a compact structure is considered to be mainly defined by interversus intracluster distances, whereas a connected structure is primarily defined by neighborhoods *H* of data. Using the definitions presented in section 2.2.1, neighborhoods can be identified based on graph theory. This can result in connected structures consisting of either unidirectional or direction-based neighborhoods.

Figure 3.4: Two types of cluster structures, compact (left) and connected (right), taken from [Handl et al., 2005]. Here, a compact structure is considered to be mainly defined by intra- versus intercluster distances, whereas a connected structure is primarily defined based on neighborhoods ܪ ሺ݇, Γ, ܯሻ and the density of the data.

An example of an algorithm that seeks compact clusters is the k-means clustering algorithm, which imposes a spherical cluster structure [Duda et al., 2001, p. 542; Handl et al., 2005, p. 3202; Hennig et al., 2015, p. 61; Mirkin, 2005, p. 108; Theodoridis/Koutroumbas, 2009, p. 742] such that the clusters cannot be too elongated [L. R. Kaufman/Rousseeuw, 2005, p. 117]. This cluster structure can be found in a data set if "the data points are actually normally distributed" (…) because "the sample mean tends to fall in the region where the samples are most densely concentrated" [Duda et al., 2001, p. 537]. The k-means algorithm is sensitive to noise and outliers [Theodoridis/Koutroumbas, 2009, p. 744]. "This drawback […] gave rise to the k-medoids algorithms […]." The PAM algorithm is less sensitive to outliers. Because of its strong similarity to the k-means algorithm, it is assumed here that PAM also yields a compact spherical cluster structure.

Examples of algorithms that seek connected clusters include density-based methods such as DBscan [Ester et al., 1996] and SL [Handl et al., 2005]. Because SL searches for nearest neighbors [Cormack, 1971, p. 331], it tends to produce connected and chain-like structures [Duda et al., 2001, p. 554; Everitt et al., 2001, p. 67; Hartigan, 1981; Jain/Dubes, 1988, pp. 64-65; Theodoridis/Koutroumbas, 2009, p. 660]. A nearest neighbor is also a Delaunay neighbor (Figure 3.4), leading to a direction-based connected structure of clusters. Spectral clustering is based on graph theory and consequently searches for connected structures [Ng et al., 2002, p. 5] of clusters with "chain-like or other intricate structures" [Duda et al., 2001, p. 582]. This indicates that such an algorithm also searches for direction-based connected clusters (see also [Hennig et al., 2015, p. 10]). "They [spectral clustering methods] are well-suited for the detection of arbitrarily shaped clusters, but can lack robustness when there is little spatial separation between the clusters" [Handl et al., 2005, p. 3202].

The Ward algorithm is sensitive to outliers and tends to find compact clusters of equal size [Everitt et al., 2001, p. 61, Tab. 1] that are ellipsoidal in structure [Ultsch/Lötsch, 2016]. The MoG method uses a mixture-of-distributions approach, which leads to connected clusters. Contrary to [Handl et al., 2005], it is argued here that the MoG method should be able to separate clusters that are non-linear separable (e.g., Chainlink [Ultsch/Vetter, 1995]). Jains and Dubes report that "fitting a mixture density model to patterns" creates clusters with hyper-ellipsoidal shapes [Jain/Dubes, 1988, p. 92]. [Handl et al.] report that the MoG method is very effective for well-separated clusters [Handl et al., 2005, p. 3202].

In the case of self-organizing mapping (SOM)15, the structures have been reported to be of "very general shapes" [Duda et al., 2001, p. 582; Ultsch/Lötsch, 2016]. Similarly to the emergent SOM (ESOM)/U-matrix clustering method [Ultsch et al., 2016a], the Databionic swarm (DBS) method that is discussed later in this work also uses the concept of emergence16, through which novel properties can arise in a system. Emergence leads to clusters whose structures are not predefined.

To summarize, the cluster structures that are theoretically sought by various methods are visualized in Figure 3.5. It should be noted that clustering methods that search for clusters with connected structures should also be able to find compact clusters as long as the distance between

<sup>15</sup> However, for k-means-SOM of the batch type, spherical or well-separated structures have been reported [Handl et al., 2005, p. 3202] (see the SOM section in chapter 4 for the differences between ESOM and k-means-SOM). 16 Definition, see chapter 7.3, p. 81-82

Figure 3.5: Overview of the cluster structures that common clustering algorithms tend to find. It is based on the literature, except for the MoG algorithm17, for which an educated guess is made. The subgroup of DBscan clustering is characterized based on arguments presented in section 3.2.1, for the definition of emergent see chapter 7.3.

clusters is large or the density between clusters is very low (see also [Handl et al., 2005, p. 3202]); e.g., "single-linkage clusters detect high-density clusters if there is a low enough valley separating them" [Hartigan, 1981]. However, methods that search for compact and spherical structures cannot be expected to find connected structures.

#### *3.2.2 Quality of Clustering*

*"[The quality of clustering is measured using a] "procedure for validating a cluster structure […]. This can be based on an internal index, an external index or resampling. An internal index scores the degree of correspondence between the data and the cluster structure. An external index compares the cluster structure with a structure given externally. A resampling is used to see whether the cluster structure is stable with respect to data change" [Mirkin, 2005, p. 205] (see also [Jain/Dubes, 1988, p. 161ff]).*

Internal and external indices are also often called *intrinsic* or *extrinsic* indices, respectively; here, they are referred to as *supervised* or *unsupervised* indices, respectively. The simplest example of a supervised index is the accuracy, which is defined as follows:

$$\text{\textbf{шири}} \text{ ол } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{ о } \textbf{u} \text{$$

$$\text{Accuracy} \left[ \% \right] = \frac{\left[ \text{No. of true positives} \right]}{\left[ \text{No. of cases} \right]} \tag{3.1}$$

In Eq. 3.1, the number of true positives is the number of labeled data points for which the label defined by a prior classification is identical to the label defined after the clustering process. To determine either the number of clusters or the clustering quality, two approaches are generally possible. Covariance matrices can be calculated, or the intra- versus intercluster distances can be compared to evaluate the homogeneity versus heterogeneity of the clusters. In the literature, a sufficient overview of 15-30 indices has already been provided [Charrad et al., 2012; Dimitriadou et al., 2002], and these indices will not be further discussed here. A special type of unsupervised indices, referred to as quality measures for projection methods, will be separately

<sup>17</sup> Also known as model-based clustering.

introduced in chapter 6. Two unsupervised indices and corresponding visualizations are presented in the following sections.

#### *<IP_ADDRESS> Heatmaps*

A heatmap is an example of an unsupervised index. For the ordering of the data points in heatmaps, dendrograms are often used. They enable the visualization of high-dimensional information and dissimilarity matrices without projecting them into a lower-dimensional space. Their use strongly depends on the sequence of the observations. For cluster validation, it is desirable to plot observations that are in the same cluster together [Hennig et al., 2015].

*"[A heatmap] consists of a rectangular tiling, with each tile shaded on a color scale to represent the value of the corresponding element of the data set. The rows (columns) of the tiling are ordered such that similar rows (columns) [in the sense that they are in the same cluster] are near each other" [Wilkinson/Friendly, 2012]. "The cluster heat map is a rectangular tiling of a data matrix with cluster trees appended to its margins. Within a relatively compact display area, it facilitates inspection of joint cluster structure" [Wilkinson/Friendly, 2009].*

Unlike in [Wilkinson/Friendly, 2009; Fig. 1], in Figure 3.7, the dendrogram between the variables is disregarded and only the ݊ݔ݊ heat map of the distance matrix is shown.

#### *<IP_ADDRESS> Silhouette plots*

The Silhouette plot is a common unsupervised index for visual evaluation of a clustering [L. R. Kaufman/Rousseeuw, 2005].

*"A score function s:* ܺ → ሾെ1, 1ሿ *evaluates the positioning of data objects inside their assigned cluster. Let a(x) denote the average distance between x and all other objects of the same cluster, and b(x) denotes the smallest average distance between x and all objects of another cluster. The silhouette score follows as* ሺݔሻ ൌ ሺ௫ሻିሺ௫ሻ ௫ሼሺ௫ሻ,ሺ௫ሻሽ *. Silhouette scores similar to 1 indicate objects that have been assigned to an appropriate cluster, whereas −1 indicates objects that have been badly classified. Silhouette scores similar to 0 indicate objects that lie in between clusters. Each cluster is represented by one silhouette, showing which objects lie within the cluster and which objects merely hold an intermediate position. The entire clustering is displayed by plotting all silhouettes into a single diagram, from which the quality of the clusters can be compared" [Herrmann, 2011, pp. 91-92].*

A reasonable clustering is characterized by a silhouette width of greater than 0.5, and an average width below 0.2 should be interpreted as indicating a lack of any substantial cluster structure [Everitt et al., 2001, p. 105]. However, it is evident that silhouette scores assume clusters that are spherical or Gaussian in shape [Herrmann, 2011, pp. 91-92].

#### **3.3 Problems with Clustering Methods**

To illustrate several problems encountered when using common clustering methods, a domain expert measured genetic data for subjects who were known either to be healthy or to have one of 3 subtypes of leukemia. Here, a typical knowledge discovery task could be to identify patterns in the cancer subtypes based on the four diagnoses leading to the prior classification**.**

 *"[I]t is a common practice among researchers to employ a variety of different clustering techniques to analyse a dataset, and to use visual inspection18 and prior biological knowledge to select what is considered the most 'appropriate' result" [Handl et al., 2005, pp. 3202-3203].*

 Consequently, the first step would be to confirm that the structure defined by the classification distinguishing the healthy patients from the non-healthy ones does indeed exist in this data set.

<sup>18</sup> The application of visual inspection will be reported in chapter 6, Fig. 1, resulting in arbitrary projections.

The data set used as an example to illustrate the general problem described above contains data representing 7747 variables for 554 subjects (see chapter 9 for details). Of the subjects, 109 are healthy, 15 have acute promyelocytic leukemia (APL), 266 have chronic lymphocytic leukemia (CLL), and 164 have acute myeloid leukemia (AML). There is a possibility that some subjects might be misclassified, but a future publication will address this diagnostic.

The heatmap and the silhouette plot presented in Figure 3.7 and 3.6 show that this data set is defined by discontinuities because the intracluster distances are small and the intercluster distances large. Hence, the leukemia data set is a high-dimensional data set with natural clusters that are specified by the illness status and defined by discontinuities19.

Table 3.1 shows the accuracies of common clustering algorithms computed by comparing the clustering results with the prior classification made available by the domain expert. The default settings were used for all algorithms, and the number of clusters was assumed to be four. The MoG algorithm cannot be applied without first using dimensionality reduction methods because the dimensionality of the data set is too high. Only one algorithm (Ward) is able to fully reproduce the prior classification. However, a classification should typically be reproduced using more than one algorithm, and the reproduction of a classification with 100% accuracy is unusual.

This example illustrates that "Clustering algorithms will create clusters whether the data are naturally clustered or purely random" [Jain/Dubes, 1988, p. 201] and "By imposing a predefined shape on the clusters, classical algorithms occasionally suggest a cluster structure in homogenously distributed data or assign points to incorrect clusters" [Ultsch/Lötsch, 2016].

To summarize, the unsupervised indices, namely, the heatmap and the silhouette plot, agree with the prior classification provided by the domain expert, whereas the external index of accuracy and the projections of the data5 disagree with the domain expert. The question arises whether this data set contains natural clusters and, if so, how the structure of these natural clusters can be correctly identified or how the optimal clustering (or projection) algorithm can be chosen for the knowledge discovery task. This work will propose approaches and solutions to these problems.

Figure 3.6: Silhouette plot of the leukemia data set indicates a cluster structure.

<sup>19</sup> It should be remarked that common data-driven methods as well as the heatmap and Silhouette plot do not reproduce the (sub) classification(s) of AML (like FAB subtypes) or CLL of research in this area, e.g. [Bene et al., 1995; Bennett et al., 1985; Vardiman et al., 2009; Haferlach et al., 2010], for CLL [Rosenwald et al., 2001].

Table 3.1: Accuracy results for common clustering algorithms. No result could be calculated for the MoG algorithm (also known as model-based clustering).

Figure 3.7: The heatmap of the leukemia data set with at least one outlier (red line). The intracluster distances are distinctively smaller than the intercluster distances. Cls1 =APL, Cls2= healthy, Cls3=CLL, Cls4=AML.

License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made. **Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **4 Methods of Projection**

Dimensionality reduction techniques reduce the dimensions of the input space to facilitate the exploration of structures in high-dimensional data. Two general dimensionality reduction approaches exist: manifold learning and projection. Manifold-learning methods attempt to find a sub-space in which the high-dimensional distances can be preserved. These sub-spaces may have a dimensionality of greater than two. However, only two- or three-dimensional representations of high-dimensional data are easily graspable for to the human observer.

The goal of this chapter is the visualization of structures in high-dimensional data. Venna et al. argued that "manifold learning methods are not necessarily good for […] visualization […] since they have been designed to find a manifold, not compress it into a lower dimensionality" [Venna et al., 2010, p. 452], and it has been shown by van der Maaten et al that they do not outperform classical principal component analysis (PCA) for real-world tasks [L. J. van der Maaten et al., 2009].

Therefore, this chapter focuses on common projection methods. Many projection methods are characterized by an objective function that is optimized using gradient descent or a corresponding learning algorithm. The quality of the projection and, consequently, of the visualization will critically depend on the similarity concept chosen as the basis of the objective function, which may be based on either distance or local proximity; thus, the methods will be categorized on this basis. This chapter will attempt to relate the various projection approaches to the compact and connected structure types introduced in the previous chapter.

# **4.1 Common Approaches**

Here, projection is used as a method for visualizing high-dimensional data in a two-dimensional space such that the discontinuities in the data are captured. Thus, the quality of a projection critically depends on the chosen similarity concept. This concept may be defined based on either distance or local proximity. The former type of similarity describes the arrangement of all given points in space and is sometimes called topography; the latter compares local neighborhoods and is sometimes called topology. Here, projections are called focusing if they are constructed using an iterative learning process that first adapts to the global intercluster distances and then focuses on more local intracluster distances.

# *4.1.1 Principal Component Analysis (PCA)*

PCA assumes that the directions in the input space that show the highest variance contain the most information about the data set [Hotelling, 1933]. The coordinate system of the input space is replaced with a (principal) coordinate system in which the variance of the data is maximized. This is achieved by finding a set of weighted linear combinations of the original variables, where the weights are found through eigendecomposition (for a definition, see [Goodfellow et al., 2016, pp. 42-44]).

Pearson proposed an equivalent definition based on an objective function in which the average projection cost is minimized [Pearson, 1901]. The projection cost is defined in terms of the mean squared distances between the points ݈∈ܫ and the projected points ݆∈ܱ:

$$E = \frac{1}{n} \sum\_{l \in I}^{n} D(l, \hat{\jmath}) \qquad \qquad \text{(4.1)}$$

where ଔ̂ൌ݆ ∑ ܾ ∗ ݑෝప ୀାଵ ൌ ∑ ܾ ∗ ݑ ୀଵ ∑ ܾ ∗ ݑෝప ୀାଵ has the same dimension as ݈ ∈ ܫ.  The minimization of *J* is achieved by choosing the basis vectors to be eigenvectors of the covariance matrix constrained by the orthonormality conditions [Duda et al., 2001, pp. 114-117]:

$$Cov(l, j) = \frac{1}{n} \sum\_{l \in I} (l - mean\_l) \left( j - mean\_l \right) \tag{4.2}$$

$$Cov(l, j) \* u\_l = \lambda\_l u\_l \tag{4.3}$$

Now, the objective function E can be redefined in (4.4) in terms of the eigenvalues ߣ in (4.1) as

$$E = \sum\_{\ell=m+1}^{n} \lambda\_{\ell} \qquad \qquad \text{(4.4)}$$

where *n* is the dimension of the input space I and m is the dimension of the output space O. The largest eigenvalues correspond to the 1, … , ݉ dimensions with the largest variance. Dimensions of the input space with small variances are discarded. Thus, PCA is an orthogonal projection of the data into a lower-dimensional space. It should be noted that "PCA remains a rather basic method and suffers from many shortcomings" [Lee/Verleysen, 2007, p. 226].

#### *4.1.2 Independent Component Analysis (ICA)*

*"Independent component analysis (ICA) is a method for finding underlying factors or components from multivariate (multi-dimensional) statistical data. What distinguishes ICA from other methods is that it looks for components that are both statistically independent, and nonGaussian" [Hyvärinen et al., 2004].*

Let ܫ ൌ ሺ݈ଵ,…,݈ሻ be defined as the matrix of the data in the input space. ICA assumes that *I* is a linear combination of non-Gaussian independent components *S* as follows:

$$I = \mathcal{S} \* A \qquad \text{ (4.5)}$$

where *A* is a linear mixing matrix and ܵ ൌ ሺ݆ଵ,…,݆ሻ,݆ ∈ܱ. ICA unmixes *I* by estimating a matrix ܹൌܣିଵ such that

$$I\*\mathcal{W} = \mathcal{S} \tag{4.6}$$

With the goal of estimating *W*, the central limit theorem and matrix search can be used to maximize the non-Gaussianity. In the fastICA algorithm [Hyvärinen, 1997], the non-Gaussianity is defined as the negentropy *F*, and it is approximately maximized by maximizing the objective function in (4.7)

$$E(j) \approx \left[ F\{ G(j) \} - F\{ G\{ N(m=0, \mathbf{s}=1) \} \} \right]^2 \tag{4.7}$$

where *N* is a Gaussian and G is a contrast function, e.g., ܩሺݑሻ ൌ െ݁ݔሺെ ௨<sup>మ</sup> <sup>ଶ</sup> <sup>ሻ</sup>.

Constraints on the estimated contrast function *G* include pre-whitening and the centering of the data in the input space [Hyvärinen et al., 2004].

#### *4.1.3 Non-linear metric multidimensional scaling (MDS) techniques*

Multidimensional scaling (MDS) was originally proposed by [Torgerson, 1952]. MDS techniques attempt to preserve the pairwise distances *D(l, j)* of the input space in the output space to the greatest possible extent. Therefore, MDS techniques minimize an objective (error) function *E* that is, as given in [Kruskal, 1964b], defined as

$$E(D,d) = \sum\_{j,l=1, j$$

where *f(D(l, j))* is a *non-metric*, monotonic transformation of the distances in the input space [Kruskal, 1964a, p. 7]. *E* is often called the stress, and *E* is minimized in an attempt to reproduce the general rank ordering of the distances. This minimization is usually performed via gradient descent.

However, the objective function *E* depends on the scale on which the distances are measured. It is preferable to normalize the objective *E* to reduce it to the same units in which the distances are expressed (Eq.4.9). Sammon mapping [Sammon] is one type of MDS technique and uses the error function

$$E(D,d) = \frac{1}{\Sigma\_{j,l=1,j<1}^n D(l,j)} \sum\_{j,l=1,j$$

#### *4.1.4 Curvilinear Component Analysis (CCA)*

When a non-linear structure is being analyzed, MDS cannot reproduce all distances. Therefore, [Demartines/Hérault] proposed a projection method that favors local neighborhoods. Curvilinear component analysis (CCA) attempts to reproduce short distances before reproducing long distances [Demartines/Hérault, 1995]. The objective function is defined in (4.10) as

$$E(D,d) = \sum\_{j,l=1, j$$

where ݄: Թ → ሾ0,1ሿ is a neighborhood function that depends on a radius R as follows:

$$h(D(l,j),R) = \begin{cases} 1, & ifD(l,j) \le R \\ 0, & otherwise \end{cases} \tag{11}$$

#### *4.1.5 t-Distributed Stochastic Neighbor Embedding (t-SNE)*

The t-distributed stochastic neighbor embedding (t-SNE) technique is an enhanced version of SNE [Hinton/Roweis, 2002] in which the Kullback-Leibler divergence (KLD) is symmetrized and the crowding problem solved. The latter is achieved by redefining the conditional probabilities in the output space O through the application of Student's t-distribution with

$$p(l|j) = \begin{cases} \frac{\{1 + \mathbf{d}(l,j)^2\}^{-1}}{\sum\_{\mathbf{l}, \mathbf{j} \in \mathbf{l}} (1 + \mathbf{d}(l,j)^2)^{-1}}, & l \neq j \\\ 0, & l = j \end{cases} \tag{4.12}$$

In [Van der Maaten/Hinton], the distance between two data points is redefined as the conditional probability that j would pick l, where ݈, ݆ ∈ ܫ, as follows:

$$P(l|j) = \begin{cases} \frac{\exp\left(-\frac{D\{l,j\}^2}{2\sigma\{l\}^2}\right)}{\sum\_{l,j\in I} \exp\left(-\frac{D\{l,j\}^2}{2\sigma\{l\}^2}\right)}, & l \neq j\\ & 0, \qquad l = j \end{cases} \tag{4.13}$$

where ߪሺ݈ሻ is the variance of a Gaussian that is centered on data point *j*. If the projection is correct, then the conditional probabilities will be equal [Van der Maaten/Hinton]. Therefore, the objective function is defined using the symmetric KLD in (14) as

$$E = \sum\_{l} \sum\_{j} \frac{P(l|j) + P(j|jl)}{2n} \* \log\left(\frac{\frac{P(l|j) + P(j|jl)}{2n}}{p(l|j)}\right) \tag{4.14}$$

#### *4.1.6 Neighborhood Retrieval Visualizer (NeRV)*

[Venna et al., 2010] reintroduced the idea of misses used by [Ultsch/Herrmann, 2005], where misses are similar data points ሺ݈ூ, ݆ூ) ∈ ݅ that are mapped onto far separated points ሺ݈ை, ݆ைሻ ∈ ܱ [Ultsch/Herrmann, 2005]. Conversely, if a pair of closely neighboring positions ሺ݈ை, ݆ைሻ represents a pair of distant data points, then this pair is called a false positive. From the information retrieval perspective, this approach allows one to define the precision ܨ and the recall ܨோ for the case in which the neighborhoods are simply binary. However, [Venna et al., 2010] goes a step further by replacing such binary neighborhoods with probabilistic ones, which are loosely inspired by the SNE approach [Hinton/Roweis, 2002]. The neighborhood of the point l is defined in terms of the relevance of the ݆∈ܫ points around l:

$$p\_l(j) = \frac{\exp\left(-\frac{D\{l, j\}^2}{\sigma\_l^2}\right)}{\Sigma\_{k \neq j} \exp(-\frac{D\{l, k\}^2}{\sigma\_l^2})}\tag{4.15}$$

where ߪ is set to the value for which the entropy of ሺ݆ሻ is equal to log(knn) and knn is a rough upper limit on the number of relevant neighbors that is set by the user [Venna et al., 2010]. The authors propose a default value of 20 effective nearest neighbors. Similarly, the corresponding neighborhood in the output space is defined as

$$q\_l(j) = \frac{\exp\left(-\frac{d(l,f)^2}{\sigma\_l^2}\right)}{\Sigma\_{k \neq f} \exp(-\frac{d(l,k)^2}{\sigma\_l^2})}\tag{4.16}$$

These neighborhoods are compared based on the mean of the KLD, which is used to define the precision ܨ and recall ܨோ:

$$F\_R = -\frac{1}{N} \Sigma\_l^N \Sigma\_{f\neq l} p\_f(l) \ast \log(\frac{p\_f(l)}{q\_f(l)}) \tag{4.17}$$

$$F\_P = -\frac{1}{N} \sum\_{l}^{N} \sum\_{j \neq l} q\_f(l) \ast \log(\frac{q\_f(l)}{p\_f(l)}) \tag{4.18}$$

The objective function is then defined in (19) as

$$E = \lambda \sum\_{\text{lj}} p\_f(\text{l}) \* \log\left(\frac{p\_f(\text{l})}{q\_f(\text{l})}\right) + (1 - \lambda) \sum\_{\text{lj}} q\_f(\text{l}) \* \log\left(\frac{q\_f(\text{l})}{p\_f(\text{l})}\right) \tag{4.19}$$

The objective function E is non-linearly optimized via conjugate gradient descent. In the absence of prior knowledge, the neighborhoods *p* are defined as symmetric Gaussians or heavytailed distributions. The weighting between precision and recall must be set by the user using the parameter ߣ. Weighting precision over recall means that if points are similar to each other in the output space, then they will also be similar to each other in the input space, whereas weighting recall over precision means that if points are similar in the input space, then they will also be similar in the output space. Note that the KLD and the symmetric KLD do not follow the triangle inequality for metric spaces.

The projection approach used in the Neighborhood Retrieval Visualizer (NeRV) method is randomly initialized by default, resulting in stochastic projections (see Figure 4.1). However, there exists an option to use PCA projection for initialization.

#### **4.2 Emergent Self-Organizing Map (ESOM)**

Self-organizing (feature) map (SOM) was invented by [Kohonen, 1982a, 1982b] and is a type of unsupervised neural learning algorithm. In contrast to other neural network models20 a SOM consists of an ordered two-dimensional layer of neurons called units. Neurons are interconnected nerve cells in the human neocortex [H. Ritter et al., 1992, p. 22], and the SOM approach was inspired by somatosensory maps (e.g. see [Hennig et al., 2015, p. 421] cites [Haykin, 1994], see also [Kandel, 2012, p. 335]). There are two types of SOM algorithms: online and batch [Fort et al., 2001]. The first is stochastic, whereas the second is deterministic, which means that it yields reproducible results for a given parameter setting. However, Fort et al. have argued "that randomness could lead to better performances" [Fort et al., 2001, p. 12].

The main differences between batch-SOM [Kohonen/Somervuo, 2002] and online-SOM [Kohonen, 1995] lie in the updating and averaging of the input data. In batch-SOM, prototypes (see Eq. 4.20 below) are assigned to the data points and the influences of all associated data points are calculated simultaneously, in contrast to online-SOM, in which sequential training of the neurons is applied (as described in detail below). The batch-SOM method has been shown to produce topographic mappings of varying quality depending on the pre-defined parametrization [Fort et al., 2001], and "the representation of clusters in the data space on maps trained with batch learning is poor compared to sequential training" [Nöcker et al., 2006]. An important comparison between the batch-SOM approach and ant-based clustering was presented by [Herrmann/Ultsch, 2008c] and will be elaborated upon in chapter 7. No objective function is used in online-SOM [Lee/Verleysen, 2007, p. 241], and SOM remains a reference tool for twodimensional visualization [Lee/Verleysen, 2007, p. 244].

In one common approach to applying the SOM concept, the algorithm acts as an extension of the k-means algorithm [Cottrell et al., 2016] or is a partitioning method of the k-means type [Murtagh/Hernández-Pajares, 1995]. In such a case, only a few units are used in the SOM algorithm to represent the data [Reutterer, 1998], which results in direct clustering of the data. Here, each neuron can be considered to represent a cluster. For example, Cottrell and de Bodt

<sup>20</sup> For an overview, see [H. Ritter et al., 1992], for deep learning see [Goodfellow et al., 2016].

used 4x4 units to represent the 150 data points in the Iris data set ([Ultsch et al., 2016a] cites [Cottrell, 1996]). Therefore, the conventional SOM algorithm is called k-means-SOM here. This SOM algorithm also has two common extensions called Heskes-SOM [Heskes, 1999] and Cheng-SOM; these two extensions include objective functions [Cheng, 1997] and are not discussed further in this thesis**.** The optimization of objective functions in general will be discussed in chapter 6, where it will be argued that it is not useful for the goal of this thesis. Chapter 7 will show that objective functions are incompatible with self-organization.

The other approach to applying SOM is to exploit its emergent phenomena through self-organization, in which case it is necessary to use a large number of neurons (>4000) [Ultsch, 1999]. This enhancement of the online-SOM approach is called emergent SOM (ESOM). In such a case, the neurons serve as a projection of the high-dimensional input space instead of a clustering, as is the case in k-means-SOM.

Let ܯ ൌ ሼ݉ଵ,… ,݉ ሽ be the positions of neurons on a two dimensional lattice21 (feature map) and ܹ ൌ ሼݓሺ݉ሻ ൌ ݓ| ݅ ൌ 1, … ݊ሽ the corresponding set of weights or prototypes of neurons, then, the SOM training algorithm constructs a non-linear and topology-preserving mapping of the input space by finding the best matching unit (ܯܤܷ (for each ݈∈ܫ:

$$bbmu(l) = \underset{m\_l \in \mathcal{M}}{\operatorname{argmin}} \{ D(l, \mathcal{w}\_l) \}, \quad i \in \{1, \dots, n\} \tag{4.20}$$

if in Eq. 4.20 a distance in the input space I between the point *l* and the prototype ݓ is denoted. In each step, SOM learning is achieved by modifying the prototypes (weights) in a neighborhood as follows:

$$
\Delta\mathfrak{w}(R) = \mathfrak{v}(R) \* h(bm\mathfrak{u}(l), m\_l, R) \* (l - \mathfrak{w}(m\_l))\tag{4.21}
$$

The cooling scheme is defined by the neighborhood function ݄: ܯ ൈ ܯ ൈ Թା → ሾെ1,1ሿ and the learning rate ߟ: Թା → ሾ0,1ሿ, where the radius R decreases until ܴൌ1 in accordance with the definition of the maximum number of epochs. In contrast to all previously introduced projection methods, no objective function is used in the ESOM algorithm. Instead, ESOM uses the concept of self-organization (see chapter 6 for further details) to find the underlying structures in data. The structure of a (feature) map is **toroidal**; i.e., the borders of the map are cyclically connected [Ultsch, 1999], which allows the problem of neurons on borders and, consequently, boundary effects to be avoided. The positions ݉∈ܯ of the BMUs exhibit no structure in the input space [Ultsch, 1999]. The structure of the input data emerges only when a SOM visualization technique called U-matrix is exploited [Ultsch/Siemon, 1990].

Let ܰሺ݆ሻ be the eight immediate neighbors of ݉ ∈ ܯ, let ݓ ∋ ܹ be the corresponding prototype to ݉, then the average of all distances between prototypes ݓ

$$u(j) = \frac{1}{n} \sum\_{l \in N(j)} D(\mathbf{w}(m\_l), \mathbf{w}(m\_l)), n = |N(j)|\tag{4.22}$$

A display of all U-heights in Eq. 4.22 is called a U-matrix [Ultsch/Siemon, 1990].

<sup>21</sup> In general this work uses the term grid if the resulting tiling is hexagonal and lattice if the resulting tiling is rectangular (see connected graph). In the context here the distinction is not important, therefore we use the term (feature) map.

*"By formalizing the displayed structures, [Lötsch/Ultsch, 2014] showed that the U-matrix is an approximation of the Voronoi borders of the high-dimensional points in the output space:*

*Let bmu(l) and bmu(j) be the BMUs of data points l and j, where bmu(j) and bmu(l) have bordering Voronoi cells. On the borderline, there is a vertical plane (AU-height), which is the distance D(l, j) > 0 between the data points in the input space. In sum, the abstract U-matrix (AU-matrix) is the Delaunay graph of the BMUs weighted by the corresponding Euclidean distances in the input space" [Thrun et al., 2016a, p. 9].*

#### *4.2.1 Visualizations of SOMs*

This section is reproduced in its entirety from [Thrun et al., 2016a]. The result of every Kohonen SOM algorithm is a set of neurons located on a map where a set *W* of prototypes corresponds to a set *M* of positions. In general, the positions on *M* are restricted to a grid/lattice, but a few approaches exist that change the positions in *M*, like Adaptive Coordinates [Merkl/Rauber, 1997]. Because these approaches are not grid/lattice based, they are not considered any further. BMUs define the locations of input points on the map. However, they exhibit no structure of the input space for a SOM [Ultsch, 1999]. However, the goal is to grasp the high-dimensional data structure and possibly even visualize cluster boundaries. Therefore, post-processing of the neurons is required for an informative representation of high-dimensional data. Three standard approaches are found in the literature:

The first approach projects the set *W* of prototypes with MDS [Torgerson, 1952] or some of its variants to a two-dimensional space [Kaski et al., 2000; Sarlin/Rönnqvist, 2013]. The result is mapped into the CIELab color space [Colorimetry, 2004]. In this uniform color space, perceptual differences in colors correspond to Euclidean distances in the map space as precisely as possible [Kaski et al., 2000]. The next two approaches visualize either the distances or density of the prototypes.

The second approach defines receptive fields around each position in *M*. The unified distance matrix (U-matrix), [Ultsch/Siemon, 1990] or one of its variants [Häkkinen/Koikkalainen, 1997; Hamel/Brown, 2011; Kraaijveld et al., 1995] , represents distances of prototypes (see equations above) by using proportional intensities of gray shades, color hues, shape or size. In [Kraaijveld et al., 1995], every neuron corresponds to a pixel. The gray value of each pixel is determined by the maximum unit distance from the neuron to its four neighbors (up, down, left, right). The larger the distance is, the lighter the gray value is. In [Häkkinen/Koikkalainen, 1997], additional unit distance visualization approaches are explained. The shapes and sizes of the receptive fields describe the dissimilarity of corresponding neurons. Apart from the U-matrix, visualizations of receptive fields in three dimensions or specific components of prototypes with receptive fields in two dimensions have been attempted [Vesanto, 1999]. It is also possible to add SOM quality measures to the receptive fields in a third dimension, e.g., [Vesanto et al., 1998].

The third approach connects the positions *M* by way of a specific scheme. In [Hamel/Brown, 2011], in addition to a U-matrix approach, neurons are connected with lines along the maximum gradient. The authors claim that clusters are the always-connected components of the graph defined by the U-matrix. [Merkl/Rauber, 1997] omitted the receptive fields approach, merely connecting map positions with lines, where the connection intensities reflect the similarity of the underlying prototypes. [K. Tasdemir/Merenyi, 2009] proposed the CONNvis technique, which visualizes the feature map by connecting neurons whose corresponding prototypes are adjacent in an input space with a dimensionality equal to that of the high-dimensional data. The

width of each connection line is proportional to the strength of the connection [K. Tasdemir/Merenyi, 2009].

In sum, all above described visualizations of large SOMs require an expert in the field for interpretation. To the best of the present author's knowledge, there are no 3D visualizations of ESOMs based on a 2D feature map currently in use22.

#### *4.2.2 Clustering with ESOM*

Combining ESOM with the U\*-matrix approach enables an application of [Ultsch et al., 2016a]:

 *"A single wall of AU-matrix represents the true distance information between two points in the data space. Valid density information at the midpoints between a BMU and a second BMU is calculated for [the] P-matrix, since the same volumes, i.e. spheres of a predefined radius, are used. The AU\*matrix therefore represents the true distance information between two points weighted by the true density at the midpoint. The representation is such that high densities shorten the distance and low densities stretch this distance. Using transitive closure for these weighted distances allows classical clustering algorithms (AU\*clustering) to actually perform distance- and density-based clustering, taking into account the complex structure of partially entwined clusters within the data."*

In contrast to the Databionic swarm approach, in which the shortest paths between AU-distances are calculated23, this clustering approach uses only the direct neighborhood of the projected points. A computation of the abstract P-matrix is necessary because ESOM itself does not consider density. Overlaying a political map on the U\*-matrix map reveals errors made by the ESOM algorithm during the annealing process. The political map shows the Voronoi areas of each cluster, where the color of each cluster area corresponds to the cluster label. The clustering is solid if every cluster consists of only one connected area, of which the borders are mountain ranges. The clustering process is sensitive to the parcel window parameter that is required for estimating the density of the high-dimensional data, and the clustering process is mostly conducted through an interactive approach requiring human intervention24.

#### **4.3 Types of Projection Methods**

In the previous section, it was shown that projection methods such as CCA, MDS and NeRV are characterized by an objective function that is optimized using gradient descent or a corresponding learning algorithm, whereas others, such as ESOM, are not. However, the first obvious difference between types of projection methods is that between linear projection methods such as PCA or ICA and non-linear projection methods. Linear projection methods are only able to rotate the high-dimensional data space and choose the most interesting dimensions, such as the dimensions with the highest variance, as is the case for PCA.

In contrast to this approach, non-linear projection methods are able to disentangle structures, e.g., represent the Chainlink data set25 in such a way that the two clusters are separated in the output space. The next major distinction between projection methods is the deterministic versus the stochastic approach. Some projection methods will always produce the same projection in the output space if all parameters remain unchanged. However, for many projection methods, such as t-SNE, their projections in the output space will drastically change with different trials

<sup>22</sup> Standard ESOM visualizations using the U-matrix are shown in supplementary D. 23 See chapter 7 for details.

<sup>24</sup> For this reason, the ESOM/U-matrix clustering approach cannot be compared with other approaches in chapter 10. 25 See the next chapter for details.

even when all settings of the projection method remain unchanged (see also examples in chapter 5, Figure 5.2). Hence, the results of deterministic methods are always reproducible, whereas stochastic methods may yield irreproducible results and require a statistical approach to assess their quality. Similarly to MDS techniques, deterministic projection methods are often based on Lyapunov functions (for further details, see [Lyapunov, 1992]). Here, it is assumed that linear and MDS techniques should only be able to visualize compact structures, which are based on the intra- versus intercluster distances of natural clusters (see the previous chapter for details).

Stochastic methods are mainly characterized by either a focusing approach or a self-organizing approach. Let k be the neighborhood extent, and let ߁ be a graph; then, a projection method is of the focusing type if the result is constructed through an iterative learning process that adapts first to global neighborhoods ܪሺ݇ଵ 1, ߁, ܫሻ and later to local neighborhoods ܪሺ݇ଶ, ߁, ܫሻ, where ݇ଵ ݇ଶ. Therefore, such methods should be capable of visualizing connected structures (see the previous chapter for details) if the annealing process is correctly chosen.

*Self-organization* is defined as spontaneous pattern formation by a system itself, without any central control26 [Kelso, 1997, p. 8 ff.]. By means of self-organization, some projection methods, such as ESOM or Pswarm, are able to project data without requiring an objective function. Thus, self-organizing methods do not implicitly predefine the structures that are sought in the data of interest. The Pswarm projection method will be introduced in chapter 8 as part of the Databionic swarm clustering approach. An overview of the various types of projection methods is shown in Figure 4.1.

Assumptions regarding the types of structures that the projection methods in Figure 4.1 are able to visualize will be either disproven or verified in chapter 10 based on 100 trials per projection method (with the exception of ICA due to technical difficulties) of five artificial three-dimensional data sets.

<sup>26</sup> Further explained in chapter 7, p.79 ff.

Figure 4.1: Overview of different types of projection methods. Here, it is argued that linear methods and MDS techniques are only able to visualize compact structures (shaded with the first pattern), whereas focusing projection methods should be able to visualize connected structures (shaded with the second pattern) if the annealing scheme is correctly chosen. For self-organizing methods, the structures that are sought in the data are not implicitly predefined. The ellipses indicate that this overview includes only common projection methods. Pswarm will be introduced in chapter 8 as a new approach based on swarm intelligence.

License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made. **Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **5 Visualizing the Output Space**

Projection methods are a common approach to dimensionality reduction with the aim of transforming high-dimensional data into a low-dimensional space. For data visualization purposes, projections into two dimensions are considered here. However, when the output space is limited to two dimensions, the low-dimensional similarities cannot completely represent the high-dimensional distances, which can result in a misleading interpretation of the underlying structures.

Nonetheless, visualization techniques based on scatter plots produced using a projection method (usually principal component analysis (PCA)) remain the state of the art in cluster analysis (e.g., [Everitt et al., 2001, pp. 31-32; Hennig et al., 2015, pp. 119-120, 683-684; Mirkin, 2005, p. 25; G. Ritter, 2014, p. 223]). Even if one disregards that "PCA remains a rather basic method and suffers from many shortcomings" [Lee/Verleysen, 2007, p. 226], visualization based on such a scatter plot is questionable in principle. Several two-dimensional scatter plots of elementary three-dimensional data sets and one high-dimensional data set (see also Figure 6.1 in the next chapter) will be presented to illustrate this claim.

Thereafter, structure preservation will be defined in this chapter to serve as the basis for a new method of visualization. This new concept with regard to the visualization of projected points in a two-dimensional output space is called the generalized U-matrix approach. In the generalized U-matrix approach, similarities between high-dimensional data are represented as valleys, and dissimilarities are represented as mountains or ridges. For the computation of the generalized U-matrix, the generation of the topographic map (see chapter 5.3) and island visualization the CRAN R package GeneralizedUmatrix was used [Thrun/Ultsch, 2017b].

# **5.1 Examples**

In Figure 5.1, the Hepta data set is shown. The Hepta data set [Moutarde/Ultsch, 2005] consists of 7 clusters that are clearly separated by distance, which means that the intracluster distances are small and the intercluster distances are large (for details, see chapter 9). This gives rise to structures that are clearly defined by discontinuity and consequently can be characterized as natural clusters.

Projections of the Hepta data set obtained by applying three of the projection methods introduced in the previous chapter are shown in Figure 5.2: PCA, curvilinear component analysis (CCA) and t-distributed stochastic neighbor embedding (t-SNE). In total, four projections are evaluated, including two t-SNE projections, denoted by *t-SNE (1)* and *t-SNE (2).* PCA yields the best representation of the clusters. With the default parameters, CCA adds excessive gaps around three points. In t-SNE (1), generated using the default parameter settings of t-SNE, the density of the data is overestimated, and wide gaps are added between two points and their corresponding cluster. When one parameter of the t-SNE algorithm is changed, resulting in t-SNE (2), the data clusters are not preserved because many random gaps are added.

Figure 5.1: The three-dimensional Hepta data set consists of 7 clusters that are clearly separated by distance. One cluster (green) has a higher density. Every cluster is ball-like in shape.

Figure 5.2: Visualizations of four cases of the projection of the Hepta data set into a two-dimensional space generated with [Thrun et al, 2017b].

> **Top left**: PCA projects the data without disrupting any clusters. This is the best-case scenario for a projection method. **Top right**: CCA disrupts two clusters by falsely projecting 3 points. This is the standard-case scenario.

> **Bottom left**: t-SNE does not correctly visualize the density of the data set at all, and one cluster is disrupted through the false projection of two points. Projection methods are often unable to correctly capture the density of data**. Bottom right**: When one parameter of the t-SNE algorithm is chosen incorrectly, all clusters are completely disrupted. This is the worst-case scenario for a projection method.

Figure 5.3: Chainlink data set and PCA projection generated with [Thrun et al., 2017]. The projection suffers from local backward projection error (BPE) and forward projection error (FPE) only in two small areas around a low number of points, but the visualization still shows low structure preservation.

The Chainlink data set [Ultsch, 2005c] consists of two clusters in Թଷ. Together, both clusters form intricate links in a chain and therefore cannot be separated by linear decision boundaries. Both rings are intertwined in Թଷ and have the same average distance and density (Figure5.3 left). The data lie on two well-separated manifolds; however, the global proximities contradict the local ones in the sense that the center of each ring is closer to some elements of the other class than it is to elements of its own class (for details, see chapter 9). PCA projection completely fails to preserve the structures in this data set because PCA merely rotates the data set and the discontinuities are not linearly separable.

#### **5.2 Structure Preservation**

Let ݇0, ݇∈Գ, let Γ be a connected graph, and let j be a point in a metric space M; then,

$$H\_l(k, \Gamma, M) = \{ l \in M \mid G(l, j, \Gamma) \le k \} \tag{5.1}$$

is the neighborhood set of j with k as the neighborhood extent, where ܩሺ݈, ݆, Γሻ is the minimum distance among all possible path distances (for details, see chapter 2, Eq. 1).