John Aldo Lee

Cyril De Bodt

Ludovic Journaux

Lucile Sautot

Proximities in dimensionality reduction

Keywords: medical imaging, radiation oncology, and medical physics dimensionality reduction, machine learning, visualization, clustering, optimization

and data mining.

-Ludovic Journaux was born in 1976 in Dijon, France. He received the Ph.D. degree in computer sciences from the

Introduction

Data collection becomes easier and cheaper every day. One purpose of accumulating large amounts of data (Big Data) and making it publicly available (Open Data) is essentially the hope that analysts will be able to extract some useful knowledge or information out of it. Information means here that data is hoped not to be just random, that some structure or patterns can be spotted and possibly interpreted. Beyond retrospective analysis of data, information is expected to be valuable in prospective applicative cases, where the same patterns could be found back or even predicted in new data. Extracting information can thus be seen as summarizing large amounts of data into exploitable, more concise knowledge.

Both structure and pattern are terms that connote some visual or geometrical perception in space. Let us leave aside text data, which are mostly linear, to focus rather on multivariate data. Then, a spatial representation is more or less straightforward. Each datum or observation can be thought of as a point (or individual sample) in a multidimensional space, with one axis (or dimension) associated with each variable (or attribute). Structure in data means that observations are not occupying space in a dense, entropic, and trivial way, like uniform or normal distributions that would span all dimensions indistinctively. Instead, some regions of space are expected to be (nearly) empty. Abstract structures that live in a data space are for instances subspaces, clusters, or outliers.

Subspaces can be linear, like hyperplanes. A straight line or a plane in our three-dimensional space are typical examples. More complex, curved hypersurfaces, are said to be nonlinear and termed manifolds. Any curve, a paraboloid, or an ellipsoid are all one-or two-dimensional manifolds embedded in our three-dimensional space (Figure 1a).

Clusters are densely populated regions of space, separated by zones where data is much more sparse (Figure 1b). To some extent, modes of a distribution can be thought of as clusters. Outliers are then rare and isolated observations lying in between or outside clusters.

Abstract structures like subspaces, clusters, and outliers are high-level objects that are loosely defined and difficult to identify in data. The task of finding them is facilitated by considering lowerlevel indicators that characterise structure. Such indicators aim for instance at quantifying possible linear and nonlinear dependences between variables, like covariance and mutual information. Other indicators deal with pairwise proximity relationships between data, like affinities, distances, or even scalar products.

Clustering and outlier detection are the task of identifying possible clusters and outliers in data, usually by labeling instances. Projection methods and manifold learning aim to spot linear and nonlinear subspaces, typically by providing a lower-dimensional representation of data. Such representation, often called projection or embedding, should ideally keep just as many dimensions as necessary. For instance, one dimension is enough to represent a curve, a plane or curved surface requires two dimensions, but more complex structures might need more dimensions to be correctly represented. Looped structures like a circle or spherical shell might need more than one or two dimensions to be rendered faithfully.

Dimensionality reduction (DR) is a task even more general than projecting or embedding [START_REF] Carreira-Perpiñán | A review of dimension reduction techniques[END_REF][START_REF] Lee | Nonlinear dimensionality reduction[END_REF]. It encompasses them, but could also achieve other goals, like forcing a twodimensional representation even though the underlying structure in data might span more than two dimensions. A typical application of DR is exploratory data visualisation. A two-dimensional representation turns out to be visually convenient, whereas the actual intrinsic dimensionality or Figure 1: Figure 1a shows a surface in shades of grey that depicts an hourglass-shaped Swiss roll, as an instance of a two-dimensional subspace embedded in our usual three-dimensional space. Figure 1b presents eight dense clusters, separated by large empty regions in the three-dimensional space.

even the true nature of the underlying structure is yet to be discovered. For instance, DR can yield a representation where not only manifolds, but also clusters and outliers can be visually detected. Beyond visualisation, DR can also be used to compress data, by keeping their principal components of variability and getting rid of the minor ones, like those only affected by statistical noise.

How does DR work intuitively? The general idea is to capture the abstract, high-level structure of data (manifolds, clusters, etc.) by evaluating appropriate low-level structure indicators (scalar products, distances, affinities, neighbourhoods, etc.) in the multivariate data space. Next, a lowerdimensional representation of data (projection, embedding, etc.) is devised such that the same indicators, now computed on the lower-dimensional representation space, yield nearly identical values. In that sense, the result of DR might be deemed to represent data faithfully. For example, clusters might be faithfully rendered if small intra-cluster and long inter-cluster distances are preserved in their low-dimensional representation. Distance preservation has been a very popular principle of DR for many decades, although it can be flawed in certain cases.

Over the course of its history, DR has moved from rather global indicators like covariance between dimensions to more and more local ones, like pairwise distances, similarities, or probabilistic neighbourhoods. Most modern DR methods focus now on preserving local proximity from the high-dimensional data space to the lower-dimensional representation. Their objective could be briefly put as: represent dissimilar data far apart and similar data close to each other. Their quality assessment follows the same principle.

Figure 2 illustrates with a toy example how modern DR methods can capture the proximities and neighbourhoods in high-dimensional data. Markers in Fig. 2a represent three-dimensional data points. Obviously, these have a complex and particular distribution in the three-dimensional space. Namely, the circles are sampled from a two-dimensional surface (half a cylinder), the triangles are distributed on a one-dimensional subspace (a spiral), the crosses define a smaller cluster, and the single square seems to be an outlier, as it is lying alone and far from other data points. In real applications, noise can affect measurements, explaining why the circles and triangles slightly deviate from their respective two-and one-dimensional subspaces.

The application of a DR method to data shown in Figure 2a should render the four types of identified structures. This goal is achieved in Figure 2b, where a two-dimensional embedding is obtained by running the Neighborhood Retrieval Visualizer (NeRV) [VPN + 10]. In particular, circles are now distributed on a piece of plane, which means that the half cylinder were nicely unrolled. The spiral is correctly represented as a one-dimensional string. The crosses still lie in a small cluster, while the square remains an outlier. Hence, as previously underlined, DR attempts to represent dissimilar data far apart in the embedding and, conversely, similar data close to each other. Here, shades of grey are used only as a visual cue to allow seeing corresponding points on figures 2a and 2b.

In a practical cases, data dimensionality can be much higher, making direct visualization very difficult. Therefore, a representation such as in Fig. 2a is not available. This highlights the interest of low-dimensional embeddings provided by DR methods, as they often give users some insights about the data structure through visualization, like in Figure 2b.

(a) (b)

Figure 2: Illustration of some concepts of dimensionality reduction. Figure 2a: Three-dimensional data space. Figure 2b: Two-dimensional embedding obtained after applying a modern DR method.

The rest of this chapter is organized as follows. Section 2 briefly reviews the history of DR and its main trends over time. Section 3 presents the different ways to categorize DR methods, considering properties of their respective models and optimization techniques. Section 4 defines the necessary symbols and concepts that are used further. After a few definitions, section 5 describes a few emblematic DR methods in more details, to illustrate their main categories and principles. Section 6 briefly discusses how to assess quality of DR, whereas Section 7 illustrates and compares DR methods with data. Section 8 deals with alternative metrics and distances. Section 9 draws the conclusions and sketches some open issues and perspectives.

Brief history of DR

The history of DR spans more than a century and it has benefitted from trends and progress in other fields like social sciences (psychometry, econometrics, and statistics in the broad sense), artificial intelligence (machine intelligence), computer science (graph theory), and applied mathematics (optimization theory).

The history of DR begins with principal component analysis, a linear projection method that has been discovered independently several times in different contexts [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF][START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF][START_REF] Karhunen | Zur Spektraltheorie stochastischer Prozesse[END_REF][START_REF] Loeve | Fonctions aléatoire du second ordre[END_REF] between 1901 and 1948. Starting from coordinates of points in some space, PCA looks for the projection that both preserves most of the observed variance and minimizes the residues between the initial points and their projection onto the linear subspace [START_REF] Jolliffe | Principal Component Analysis[END_REF]. PCA is described in Section 5.1.1.

Similar to PCA, classical metric multidimensional scaling (MDS) aims to project points that are known through pairwise relationships like similarities or dissimilarities [START_REF] Young | Discussion of a set of points in terms of their mutual distances[END_REF][START_REF] Torgerson | Multidimensional scaling, (i): Theory and method[END_REF][START_REF] Gower | Some distance properties of latent root and vector methods used in multivariate analysis[END_REF]. In the particular cases where these relationships are scalar products or Euclidean distance, classical metric MDS turns out to yield the same solution as PCA. Section 5.2.1 deals with this early version of MDS.

Both PCA and classical metric MDS are considered to project data in a usual, linear way, if distances are Euclidean. Starting from the late 1960s, metric MDS, also known as stress-based MDS, generalizes classical MDS and allows for nonlinear mappings of data, thanks to explicit formalization of a very generic and tunable criterion of distance preservation from one space to the other [START_REF] Shepard | The analysis of proximities: Multidimensional scaling with an unknown distance function (parts 1 and 2)[END_REF][START_REF] Sammon | A nonlinear mapping algorithm for data structure analysis[END_REF][START_REF] Demartines | Curvilinear component analysis: A self-organizing neural network for nonlinear mapping of data sets[END_REF]. This criterion is an objective function that is often called stress, which assesses how pairwise distances d ij computed from the low-dimensional coordinates match the corresponding distances δ ij in the high-dimensional data space. In particular, this framework allows for putting more weight on the preservation of short distances as compared to longer ones. This has been the first attempt to implement the intuition that proximities [START_REF] Shepard | The analysis of proximities: Multidimensional scaling with an unknown distance function (parts 1 and 2)[END_REF] are more important than more distant, looser relationships [Sam69, DH97].

In the sixties and seventies, yet another generalization of MDS was investigated, with the early intuition that distances measured in two spaces of largely differing dimensionalities are not necessarily comparable [START_REF] Donoho | High-Dimensional Data Analysis: The Curse and Blessings of Dimensionality. Aide-mémoire for a lecture for the American Math[END_REF][START_REF] François | The concentration of fractional distances[END_REF]. Nonmetric MDS relaxes strict distance preservation (δ ij = d ij for all i, j) into ordinal distance preservation (δ ij ≤ δ ik ⇒ d ij ≤ d ik for all i, j, k) [START_REF]Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis[END_REF][START_REF] Takane | Nonmetric individual differences multidimensional scaling: an alternating least squares method with optimal scaling features[END_REF]. Section 5.2.2 describes stress-based (metric and nonmetric) MDS.

In the eighties, artificial neural networks raised much interest, with a peak when backpropagation was popularized in 1986 [START_REF] Rumelhart | Learning representations by backpropagating errors[END_REF]. A few years earlier, a neural model known as self-organizing feature maps (SOMs) was published [START_REF] Kohonen | Self-organization of topologically correct feature maps[END_REF]. This model performs both vector quantization (or clustering, like K-means) and a particular form of nonlinear DR, particularly suited for data visualization. The neural trend continued in the nineties with the advent of auto-associative networks, composed of an encoder and decoder in sequence [Kra91, Oja91, UNN91, DC93]. Difficulties to train such networks with many layers of neurons stacked in the encoder and decoder were addressed only recently in the early years of the deep learning trend [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. Section 5.3 investigates these neural approaches of DR based on vector quantization or auto-association.

In the late nineties and early 2000s, developments in DR were influenced mainly by the interest in kernel techniques in the booming field machine learning. Briefly put, the so-called kernel trick theorizes the fact that some particular functions termed kernels induce a scalar product in an unknown, possibly very high-dimensional space, without dealing explicitly with coordinates (or features) in that space. The pioneering work in that direction was kernel PCA, which 'kernelizes' classical metric MDS into a method of nonlinear mapping [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]: MDS is applied in the feature space induced by a nonlinear kernel, instead of the raw data space. Other methods followed, investigating numerous kernels and various ways to build them, like graph-approximations of geodesic distances [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], commute-time distance [YVW + 05], (heat) diffusion kernels [CLL + 05], or even optmized kernels [START_REF] Weinberger | Unsupervised learning of image manifolds by semidefinite programming[END_REF]. All these methods are closely related to graph theory [BN02, RS00, DG03] and spectral clustering [SFYD04, BVP + 03]. While these methods are appealing theoretically, they entail usually some arbitrariness in the choice of the kernel or difficulties in its optimization. Section 5.4 deals with spectral methods of DR, with a focus on kernel PCA and Laplacian eigenmaps.

At the time of writing, the latest trend in DR investigates methods of stochastic neighbour embedding (SNE), which started in 2002 [START_REF] Hinton | Stochastic neighbor embedding[END_REF]. The next milestone occurred a few years later, in 2008, with t-distributed SNE, a variant of SNE that gained popularity thanks to its relative simplicitly and disruptively good experimental results [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF]. Stochastic neighbour embedding aims at preserving probabilistic neighbourhoods from one space to the other: if some point ξ j is deemed to be a neighbour of ξ i with good probability in the HD space, then its LD counterpart x j should also be a neighbour of x i with approximately the same probability. The properties of these stochastic neighbourhoods have been investigated [START_REF] Hinton | Stochastic neighbor embedding[END_REF][START_REF] Lee | Shift-invariant similarities circumvent distance concentration in stochastic neighbor embedding and variants[END_REF][START_REF] Vladymyrov | Entropic affinities: Properties and efficient numerical computation[END_REF][START_REF] Lee | Two key properties of dimensionality reduction methods[END_REF], as well as objective functions to assess their HD-LD discrepancy [VPN + 10, BHBV12, LRB + 13]. Variants with simplified neighbourhoods [START_REF] Carreira-Perpiñán | The elastic embedding algorithm for dimensionality reduction[END_REF], multi-scale neighbourhoods [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF] and missing data robustness [START_REF] De Bodt | Nonlinear Dimensionality Reduction with Missing Data using Parametric Multiple Imputations[END_REF] exist, as well as scalable approximations of SNE with tree structure [START_REF] Hinton | Stochastic neighbor embedding[END_REF][START_REF] Van Der Maaten | Accelerating t-SNE using tree-based algorithms[END_REF][START_REF] Yang | Scalable optimization of neighbor embedding for visualization[END_REF][START_REF] De Bodt | Extensive assessment of Barnes-Hut t-SNE[END_REF]. Section 5.5 presents an overview of probabilistic neighbourhood preservations, with details on SNE and its variants.

In parallel to the development of new DR methods, the issue of DR quality assessment (QA) has been investigated at some occasions. It is obviously important but difficult, DR being an unsupervised learning task. The most obvious way to evaluate quality of DR results is to look at how well the criterion of the considered DR method has been optimised. Popularity of stressbased MDS has led to the common adoption of stress functions, like Sammon's stress, as quality indicators [START_REF] Sammon | A nonlinear mapping algorithm for data structure analysis[END_REF]. This approach is, however, intrinsically biased towards methods working with the very same objective. For instance, specific quality criteria have been devised for SOMs [VDHM97, BHV99, GS96, GS97, BP92]. Indirect ways to assess DR quality are also widespread, like classification errors, provided data are class-labeled. In the recent years, generic QA frameworks have been proposed, most of them quantifying neighbourhood preservation, as this principle makes only minimal assumptions about the data and method at hand [VK01, AC06, FC07, VK07, VPN + 10, LV09, LRB + 13, LPOnV15].

Categories

Dimensionality reduction is generic objective that can be particularized in many ways. Also there are different models and techniques to achieve these objectives, leading to methods with specific features.

Linear versus nonlinear model

Dimensionality reduction amounts to finding a subspace that can 'explain' observed data and dependencies between dimensions or variables. These dependencies can be either simply linear or more complicated, namely, nonlinear. In other words, data are distributed on some linear or nonlinear subspace (or close to these if there is some noise). Technically, a nonlinear subspace is called a manifold.

Some DR methods rely on a linear model: they can only identify a linear subspace and project data onto it. More complex DR methods can disentangle nonlinear dependencies between variables and are thus capable of identifying manifolds. Nonlinear DR methods are more powerful than linear ones, but also more prone to overfitting. Nonlinear DR can somehow unfold manifolds before projecting them in a lower-dimensional space; this is sometimes abusively called nonlinear projection.

The assumption of an underlying manifold, also known as the manifold hypothesis, is sometimes deemed central in some DR methods, called manifold learners. Most of the time, however, these methods are applied like others to all sort of data, be there a manifold, or clusters.

The manifold hypothesis calls for another concept, called the intrinsic dimensionality of the manifold [START_REF] Fukunaga | Intrinsic dimensionality extraction[END_REF]. If data is provided as coordinates in some space, the dimensionality of this space is just the number of features or variables. On the other hand, the dimensionality of the manifold is a priori unknown. For instance, data distributed on a sphere has an ambient dimensionality of three and an intrinsic dimensionality of two (the sphere is a surface in the physical space). There are various methods of estimating the intrinsic dimensionality of data [START_REF] Camastra | Estimating the intrinsic dimension of data with a fractal-based method[END_REF][START_REF] Fukunaga | Intrinsic dimensionality extraction[END_REF].

Parametric versus nonparametric model

The distinction between these two types of models concerns mostly how the mapping from high-to low-dimensional coordinates is formalized and how it can be reused. In non-parametric methods, the model reduces to its simplest form: the DR method gives each available datum low-dimensional coordinates. This means that if new data come afterwards, the model will not be able to compute their associated low-dimensional. A complementary method, called an out-of-sample extension [BPV + 04], can possibly carry out this additional task.

On the other hand, parametric models usually assume some analytical function that formalizes the mapping from high-to low-dimensional coordinates. Typically this function is made as generic as possible with the use of many tunable parameters and DR then just consists in adjusting those. Parametric models can be directly reused to new, out-of-sample data. The price to pay is that no matter how generic they are designed, parametric models always rely on some assumption about the data mapping, that can be restricted to be linear like in PCA [START_REF] Jolliffe | Principal Component Analysis[END_REF] or constrained to be continuous like in auto-encoders [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF].

Local versus global DR

Very often, DR is categorized into local and global methods [START_REF] Silva | Global versus local methods in nonlinear dimensionality reduction[END_REF][START_REF] Roweis | Global coordination of local linear models[END_REF]. However, these opposed categories are rather unclear and fuzzy. Most of the time, local methods are those who consider only short-range distances or neighborhoods, with respect to the total number of data, whereas global ones take into account all distances and neighborhoods indifferently. Not considering explicitly long distances (or non-neighbors) does not necessarily mean that the model or method does not implicitly give them some default value, like infinity. In spectral methods like Laplacian eigenmaps [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF] or locally linear embedding [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], usually seen as local, reformulating them as classical metric MDS applied to commute-time distances transforms them into equivalent methods that would typically be considered global. Moreover, the local/global distinction appears mostly informal and can difficultly account for hybrid method like MDS and Sammon's mapping, where the importance of short and long distances can be modulated. A better way of determining whether a DR method or embedding is local or global would be to look at how it performs for small and large neighborhoods with the quality assessment tools detailed further below in Section 6. [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF].

Discriminative versus generative DR

This distinction deals with the direction the model or method maps data. A discriminative method typically maps high-dimensional data to a lower-dimensional embedding. Most methods belong to this category. A generative method works the other way around and starts from the lowdimensional space: a candidate embedding and the mapping to the high-dimensional space are refined in order to be able reproduce the observed data.

Discriminative DR can be non-parametric, whereas generative DR is typically either parametric or non-parametric but with a user-fixed low-dimensional space, like in SOMs. Some methods like auto-encoders provide parametric mappings in both directions (which are not necessarily the inverse of one another).

Mild versus about-right versus harsh DR

Dimensionality reduction can process data with largely varying numbers of features or coordinates, between three and thousands, and aim typically at spaces with one, two, three, or a few more dimensions. Hence the dimensionality gap between the initial and final spaces can vary a lot [START_REF] Carreira-Perpiñán | A review of dimension reduction techniques[END_REF].

About-right DR would target an embedding dimensionality that is guided by the estimated intrinsic dimensionality of data [START_REF] Takens | On the numerical determination of the dimension of an attractor, in dynamical systems and bifurcations[END_REF], under the manifold hypothesis detailed in Section 3.1 above. If there are for instance five latent parameters or degrees of freedom in the data generation process, then at first glance it would be counterproductive to reduce dimensionality below five, at the expense of squeezing data and possibly damaging part of its structure. Nevertheless, five dimensions remain difficult to explore visually and harshly crushing data down to three or two dimensions, whatever the actual intrinsic dimensionality is, might still be a useful tradeoff and provides some insights. Eventually, the converse is possible as well: if visualization is not intended, if preserving the data structure is essential, or if DR is just some form of preprocessing to compress data, then mild DR can target a dimensionality that is slightly higher than the intrinsic one [START_REF] Takens | On the numerical determination of the dimension of an attractor, in dynamical systems and bifurcations[END_REF].

Spectral versus nonspectral DR

Dimensionality reduction methods can also be categorized according to the optimization techniques they involves. Two main categories can be distinguished. The first one includes all methods that solve a convex optimization problem, typically with an spectral decomposition into pairs of eigenvalues and eigenvectors, like in PCA, Laplacian eigenmaps, LLE, and many other so-called spectral embedding methods [BH03, XSB06, SWS + 06].

The main advantage of solving a convex problem is the theoretical guarantee to get the global optimum. In the particular case of eigendecomposition for DR, the embedding dimensionality can be adapted very easily without reoptimizing: eigenvectors and eigenvalues are associated with dimensions and can be easily considered or discarded. On the other hand, the main shortcoming of convex optimization is that it usually forces working within a tight frame with little flexibility. Problem statement can be much more flexible when the convexity constraint is relaxed. The price to pay is that non-convex optimization can get stuck in (or close to) a local optimum. The hope is then that the local optimum of a problem stated without convexity constraint proves to be better than the global optimum of a problem where convexity is key in the statement. This remains an open and somewhat subjective debate but, currently, non-spectral (an thus non-convex) nonlinear DR methods seem to outperform spectral embedding [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF]. Non-spectral DR methods typically involve generic optimization techniques that are based on derivatives, like (stochastic) gradient descent.

Paradigms of DR

Yet another even more obvious way to categorize DR methods consists in distiguishing various principles or paradigms that guide their conception. For instance, PCA is often considered as a method that best preserves variance from the data space to the projection space. But PCA can also be seen as minimizing the error in a reconstruction problem, where data dimensionality is first reduced and then recovered in a sort of encoding and decoding sequence. The same idea drives auto-associative networks (also known as auto-encoders), except that they substitute linear projection and reconstruction operators with nonlinear mappings.

Another important paradigm in DR is the preservation of pairwise proximities, typically in terms dissimilarities or similarities. These can be scalar products, like in classical MDS [START_REF] Torgerson | Multidimensional scaling, (i): Theory and method[END_REF], or distances, such as in nonlinear variants of MDS [START_REF]Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis[END_REF]. Distance preservation can aim at perfect isometry for all distances, focus chiefly on small distances thanks a weighting scheme [Sam69, DH97, WS06], or keep track of ordinal distance relationships, like in nonmetric MDS [START_REF]Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis[END_REF]. The measured distance can be Euclidean, geodesic (or a graph approximation thereof), some implicitly defined distance defined by a kernel [SSM98, BN02, WS06], like the commute-time distance, or even some other metric.

To put the emphasis on local data structure, the use of similarities seems however more natural than distances. Distances grow and have significantly large values when points get far apart, whereas similarities are negligible in that case and become large when points get close to each other. Most modern successful DR methods, be they spectral or non-convex, rely rather on similarities than dissimilarities, like Laplacian eigenmaps, LLE, SNE, and their variants. Preservation of proximities, particularly in the form of (soft) neighborhoods, is thus the principle that proves to be the most efficient in DR currently.

Definitions

This sections defines some mathematical notation and a few concepts that will be used afterwards.

Scalars, vectors, and matrices

In this chapter, symbols in italic denote scalar numbers, like variable x, index i, parameter α, or constant N . Bold lower-case symbols denote column vectors, like x = [x k ] 1≤k≤M ∈ IR M , where M is the dimensionality of x, namely, its number of components or entries. The transpose of x is a row vector noted x T . Bold upper-case symbols denote matrices, like

A = [a ij ] 1≤i≤M,1≤j≤N ∈ IR M×N . Matrix A can be seen as collection of indexed column vectors, as in A = [a i ] 1≤i≤N . The transpose of A is A T = [a ji ] 1≤j≤N,1≤i≤M
. For two equal-size vectors u and v, their scalar product is computed as

u T v = v T u = M k=1 u k v k .
Similarly, the matrix-vector and matrix-matrix products are computed as

y = Ax = K k=1 a ik x k 1≤i≤M and C = AB = K k=1 a ik b kj 1≤i≤M,1≤j≤N
.

Data sets and coordinates in space

Let Ξ = [ξ i ] 1≤i≤N denote a set of N points in some M -dimensional data space, with M potentially very high. Each ξ i is a column vector from IR M holding the coordinates. Similarly, let X = [x i ] 1≤i≤N be the data set representation in a P -dimensional space, with P ≤ M . All along this chapter, Greek symbols usually refer to quantities that are relative to the high-dimensional (HD) space, whereas usual latin symbols denote their counterpart in the low-dimensional (LD) representation space.

Objective function and optimisation

The goal of DR is to arrange points x i ∈ IR P , with P potentially much smaller than M , such that X is representative of Ξ, namely, its preserves most of salient or relevant structural features. Representativeness of X with respect to Ξ is quantified formally by some mathematical function E(X; Ξ, θ), termed objective function, cost function, energy function, or stress function. Vector θ denotes some additional parameters that could modulate the function behavior. Usually, representativeness is made optimal by minimizing the value of E. This is achieved by a clever trialand-error procedure that progressively updates X, like (stochatistic) gradient descent or spectral decomposition.

Covariance and Gram matrices

Assuming that the variables in Ξ and X are centered, the sample covariance matrices in the HD and LD spaces are C Ξ = 1 N ΞΞ T and C X = 1 N XX T , with sizes M -by-M and P -by-P , respectively. Superscript T denotes the transpose of a matrix. In product

C = A T B, element c ij from C is computed as c ij = k a ki b kj .
The N -by-N Gram matrices of all pairwise scalar products in the HD and LD spaces are G Ξ = Ξ T Ξ and G X = X T X, respectively. Covariance assesses the linear relationships between two rows of Ξ or X, i.e., two variables, whereas the Gram matrix assesses pairwise relationships between columns of Ξ or X, i.e., observations.

Spectral decomposition

Square symmetric matrix A can be factorized into A = VΛV T , where V is orthogonal (V T V = I) and Λ is diagonal. The columns v k of V are the eigenvectors of A, while λ kk are their associated eigenvalues.

Let us consider some known, square, symmetric matrix Y that results from the products XX T or X T X, where X is unknown but has the same size as Y. For such products, all eigenvalues are guaranteed to be non-negative (0 ≤ λ kk ). The spectral decomposition would allow us to retrieve X by computing X = VΛ 1/2 or X = Λ 1/2 V T , respectively, where Λ is diagonal with entries λ -1/2 kk . To some extent, X can thus be interpreted as the square root of Y.

Norms and distances

The norm of some vector u is its length, noted u . In particular, the Euclidean norm is defined as

u 2 = √ u T u = k u 2 k .
The Frobenius norm extends the Euclidean norm to matrices, namely,

A F = ij a 2 ij .
The distances between the ith and jth points are defined as δ ij = ξ i -ξ j and d ij = x i -x j in the HD and LD spaces, respectively. Most of the time the distances are Euclidean, namely,

δ ij = ξ i -ξ j 2 = M m=1 (ξ mi -ξ mj ) 2 and d ij = x i -x j 2 = P p=1 (x pi -x pj ) 2 .
(1)

Distances are obvious ways to assess proximity between two vectors. For matrices A and B, the Frobenius distance A -B F can be used.

DR paradigms and emblematic methods

Variance preservation

The main representent in this category is principal component analysis (PCA), detailed here below. Many variants of PCA exist, like minor component analysis (focusing on axes with the least variance) or linear discriminant analysis (favoring axes that best separate classes in labeled data).

Principal component analysis

The intuition behind PCA can be traced back to a multivariate regression problem. Regression is a supervised learning task in which a model is trained to predict outputs y from inputs x. While scalar outputs y i are often considered, nothing prevents us to have multiple outputs like here in vector y. The model can be denoted by some parametric function ŷ = f θ (x), where ŷ is the predicted output. Performance is typically measured in terms of squared residues (or least squares) through

y i -ŷi 2 2 . For a complete data set (X = [x i ] 1≤i≤N , Y = [y i ] 1≤i≤N ), the mean squared error is then E(θ; X, Y) = 1 N N i=1 y i -ŷi 2 2 = Y -Ŷ 2 F . (2)
Let us now consider an apparently trivial problem where only input coordinates are available. Instead of having both X and y, we only have Ξ and no counterpart for y. In this case, we would like to approximate Ξ itself, which is expressed as

E(θ; Ξ) = 1 N N i=1 ξ i -ξi 2 2 = Ξ -Ξ 2 F , (3)
where ξ = f θ (ξ). A trivial solution to this reconstruction problem consists in choosing the identity function for f θ . To avoid triviality, let us assume that f θ consists actually of two successive functions e θ and d θ , namely ξ = d θ (e θ (ξ)), with x = e θ (ξ) and ξ = d θ (x), with the constraint that dimensionality of x is lower than that of ξ. particular coefficient (not shown). When several straight lines join together at the left side of some squares, it means that the corresponding component is obtained by adding up the components linked by the lines, weighted by their coefficients. Hence x 1 and x 2 are linear combinations of ξ 1 , ξ 2 and ξ 3 , and ξ1 to ξ3 are linear combination of x 1 and x 2 . The set of straight lines joining ξ (resp. x) to x (resp. ξ) encode e θ (resp. d θ ) as a linear transformation. The whole network tries to reproduce the identity function, under the constraint of the bottleneck x. The dimensionality reduction is performed by projecting ξ to x.

In PCA, functions e θ and d θ are thus chosen to be simple linear transformations, which somehow correspond to changes of coordinate systems. More precisely, only shifts, rotations, and mirroring operations are considered. This amounts to considering some linear subspace of the data space, like in Fig. 4. The axes [x 1 , x 2 ] define a 2D subspace of the 3D data space [ξ 1 , ξ 2 , ξ 3 ]. They are positioned such that projection residues are minimal in (3). In other words, PCA finds a linear projection subspace that maximizes proximity between original data and their reconstruction after projection. If rotations (and mirroring) are centered on the data set average, then the sum of squared residues can also be interpreted as a variance, computed perpendicularly to the x subspace, following the residue lines. Let us call it lost variance. Under the same assumption, variance can of course also be computed along the x axes, within the considered subspace. Let it be the preserved variance. Because the projection and backprojection then consist of rotations only, scales remain unchanged in all directions and the total variance, which is computed as the sum of the lost and preserved variances, remains constant. Rotating the subspace only transfer variance from lost to preserved or the other way round, like communicating vessels. PCA chooses the subspace with both minimal reconstruction residues and maximum projected variances.

In practice, PCA can be computed in several ways, with a first common step, which is always data centering. Subtracting the data set mean can indeed be proved necessary to minimize residues. The mean can always be added back after reconstruction. To compute the optimal rotation, one can either use the singular value decomposition of the centered data coordinates or the EVD of the covariance matrix. As a reminder, the sample covariance is written as 1 N ΞΞ T if Ξ is centered (i.e., Ξ1 T = 0). The SVD of the centered data set is Ξ = VΣU T . The P -dimensional PCA projection is then given by X = U T P , namely, the first P columns of U T . The EVD of the sample covariance is 1 N ΞΞ T = VΛV T . The P -dimensional PCA projection is then given by X = V T P Ξ, where V T P gathers the eigenvectors associated with the largest P eigenvalues. Reconstruction is carried out by Ξ = V P X. PCA's main limitation is the linearity of the projection and back-projection it involves, which is suited only to linear dependencies among data features. Nonlinear extensions can rely on a discrete, nonlinear subspace, like the self-organizing maps in Section 5.3.1 or nonlinear encoders and decoders implemented with artificial neural networks, like auto-encoder networks in Section 5.3.2

Distance preservation

Hereafter are presented various flavors of distance-preserving methods, after a short detour by classical metric multidimensional scaling. This method makes the transition between PCA and

ξ 1 ξ 2 ξ 3 x 1 x 2 (a) x 1 x 2 (b)
Figure 4: Illustration of PCA as a dimensionality reduction method finding the linear subspace preserving the most variance of the original data. Figure 4a: Three-dimensional data space [ξ 1 , ξ 2 , ξ 3 ]. The data points are illustrated by dots. x 1 and x 2 are the first two principal components of the data set. This means that the grey shaded plane, defined by axes x 1 and x 2 , is the linear subspace of dimension two preserving the most variance of the data points in three dimensions, once orthogonally projected on the plane. Figure 4b: Two-dimensional embedding obtained, which simply consists in the grey shaded plane of Figure 4a, on which the data points were projected. distance preservation. It relies on pairwise scalar products, which are mathematically related to both covariances and Euclidean distances.

Classical metric multidimensional scaling

Compared to PCA, the idea behind classical metric multidimensional scaling (CMDS) is that data should be described by quantities characterizing the proximity relationship for pairs of data, rather than the proximity between pairs of variables, like covariances.

The simplest pairwise quantity that can be defined is the scalar product (or dot product, or inner product). For data ξ i and ξ j , the usual scalar product can be written as ξ T i ξ j = ξ T j ξ i = k ξ ki ξ kj . The interpretation of scalar product is not very intuitive and relies mainly on angles. If the scalar product between ξ i and ξ j is null, then they are 90 degrees from each other with respect to the origin 0. If the scalar product is positive (resp. negative), then the angle is greater (resp. lower) than 90 degrees. For N data vectors in Ξ, all pairwise scalar products can be gathered in a symmetric N -by-N matrix called Gram matrix:

G Ξ = Ξ T Ξ.
As in PCA, the data set Ξ can be assumed to be centered. If it is not, then its mean 1 N Ξ1 = 1 N N j=1 ξ j can be subtracted from each ξ i . After centering, we have

N i=1 ξ i = 0.
The idea of CMDS is to reproduce as well as possible in the LD space the dot products that are observed in the HD space. Formally, this can be written as min

X G Ξ -G X 2 F = min X G Ξ - X T X 2
F , which is usually called the strain function. The solution relies on positive semidefiniteness of Gram matrices and the eigenvalue decomposition G Ξ = VΛV T , where V T V = I and Λ is diagonal. We have min X VΛV T -X T X 2 F . For target dimensionality P , the solution is then X = Λ 1/2 P V T P . Subscript P indicates that that we keep only the first P columns of V and the first P rows and columns of Λ.

What can be done if Ξ is not centered? Then, the centering matrix

C = I -1 N 11 T can be used. It is easy to see that XC = X(I -1 N 11 T ) = X -( 1 N X1)1 T is centered, since the data set mean 1
N X1 is explicity subtracted from all columns. The main interest of the centering matrix is that it can be applied directly to the Gram matrix. Double centering consists in computing CG Ξ C = C T Ξ T ΞC (since C is symmetric and invariant under transpose).

Centering is essential here since dot products are seen as a form of angular proximity measure with respect to the origin 0. If the data set is not centered, the origin might lie far away and then all angles would be small. On the other hand, if data is centered, then angles are likely to vary between 0 and 180 degrees, leading therefore to a more discriminative similarity. Centering also reinforces the duality between PCA and CMDS. For centered data, one might demonstrate that the eigenvalues and eigenvectors of either the covariance matrix 1 N ΞΞ T in PCA or the Gram matrix Ξ T Ξ in CMDS are closely related; the key to this duality is to look at the SVD of Ξ after centering.

Eventually, it is worth noticing that CMDS can also be applied to distance matrices as depicted in Figure 5, in addition to Gram matrices or raw centered coordinates. Let ∆ 2 = [δ 2 ij ] 1≤i,j≤N denote a matrix of squared pairwise distances. Then, double centering on ∆ 2 , namely, C∆ 2 C, can be shown to yield a valid Gram matrix (provided δ ij has all properties of a genuine distance). This can be easily seen in the case of the squared Euclidean distance, precisely defined with respect to scalar products:

δ 2 ij = ξ i -ξ j 2 2 = ξ T i ξ i -2ξ T i ξ j +ξ T j ξ j .
(Notice that distances are shift invariant and thus remain the same whether data are centered or not.) The effect of double centering is then to substract the left and right terms, leaving only -2ξ T i ξ j . Multiplying further by -1 2 yields the elements of the Gram matrix. Figure 5: Illustration of classical metric multidimensional scaling. Figure 5a: Three-dimensional space, in which a two-dimensional subspace corresponding to a half cylinder is embedded. Three data points (black dots) are lying on this subspace. The two solid black lines depict the Euclidean distances between two dots and the third one. The dashed line is the shortest path on the half cylinder between the two joined data points. Figure 5b: Two-dimensional embedding. The lowdimensional representation of the data points aims to preserve as much as possible the pairwise distances computed between their high-dimensional counterparts.

Stress-based metric multidimensional scaling

In its classical form, multidimensional scaling involves scalar products to characterize the (angular) proximities between data. Usual scalar products have a strong relationship with squared Euclidean distances, namely,

δ 2 ij = ξ i -ξ j 2 2 = ξ T i ξ i -2ξ T i ξ j + ξ T j ξ j .
From that perspective, CMDS can thus be seen not only as direct way to preserve scalar product (min X Ξ T Ξ -X T X F ), but also as an indirect, approximate way to preserve squared Euclidean distances.

If ∆ 2 = [δ 2 ij ] 1≤i,j≤N and D 2 = [d 2 ij ]
1≤i,j≤N denote the matrices of pairwise squared Euclidean distances, then squared distance preservation would be written as min X ∆ 2 -D 2 F . Notice that this problem can no longer be solved with spectral techniques like eigenvector decomposition. Instead, more generic optimisation techniques are needed, typically based on derivatives, like gradient descent. Another side effect of this change, perhaps even more important, is that the problem is no longer linear. One cannot guarantee anymore that the solution in X results from a linear transformation of Ξ as in PCA and CMDS.

Considering distances instead of angular proximity proves to be easier for the intuition, short distances translating immediately the notion of proximity. Starting from there, MDS could be extended to address the following cases:

• How raw (non-squared) Euclidean distances or even non-Euclidean distances can be preserved?

• How to put more emphasis on some pairwise relationships than on others?

The goals here are to broaden the applicability of MDS to other distances and to refine its objective, when for instance short distances (local structure) are deemed more relevant than longer ones (global structure). For this purpose, let us get rid of matrix notation and rewrite squared distance preservation as

min X ∆ 2 -D 2 F = N i,j=1 (δ 2 ij -d 2 ij ) 2
. This cost function is called sstress, standing for squared stress. This formula involves distances to the 4th power, meaning that large distances have much heavier weight than short ones. At the same time, short distances are those that describe local proximities and neighbourhoods. To make them gain some importance, a first workaround consists in using raw, non-squared distances. The resulting stress function can be written as min X i,j (δ ij -d ij ) 2 . Another workaround is to introduce explicit weights, as in min X i,j w ij (δ ij -d ij ) 2 , where 0 ≤ w ij . An emblematic choice for the weights is w ij = 1/δ ij , like in Sammon's nonlinear mapping [START_REF] Sammon | A nonlinear mapping algorithm for data structure analysis[END_REF]. Such weights explicitly grant short distances more importance. For example, twice longer a distance has half the weight. Notice that weights and importance are determined here by δ ij , the distance in the HD data space. The weighting scheme can be changed to take into account the distance in the LD representation space instead, leading to a method called curvilinear component analysis (CCA) [START_REF] Demartines | Curvilinear component analysis: A self-organizing neural network for nonlinear mapping of data sets[END_REF]. The weight is then typically written as

w ij = H(λ -d ij )
, where H(u) is Heaviside's step function (null if u < 0, one elsewhere) and λ is parameter controlling the radius of the neighbourhood around each x i . CCA is more flexible than Sammon's mapping, in the sense that if d ij gets larger than λ, then the weight gets null and the corresponding term in the cost function vanishes. In other words, CCA tolerates some distance stretching, whereas Sammon cannot. Distance stretching is useful to unfold a curved manifold or even to tear apart some piece of the manifold when necessary (typically to break loop structures [START_REF] Lee | How to project "circular" manifolds using geodesic distances[END_REF] like circles and sphere). A planisphere of the Earth is the typical example where this property comes in very handy.

Eventually, distance preservation has limitations when it comes to reduce very high-dimensional data (tens of dimensions) to very low-dimensional representations (typically, 2D), because distances have completely different statistical properties in those high-and low-dimensional spaces. Distances are said to concentrate in HD spaces. This concentration gap makes the direct comparison of distances, as involved in stress-based MDS, hardly relevant. This issue is addressed partly by non-metric MDS, hereafter, and more thoroughly by methods of stochastic neighbour embedding, as described further in Section 5.5.

Non-metric multidimensional scaling

An even more evolved form of MDS is called non-metric or ordinal MDS. In this variant, distance preservation is expressed as min X,f i,j (f (δ ij )-d ij ) 2 , where f is a growing function (or monotonic, with a positive derivative). In this problem, d ij is not told to match δ ij as such. Instead, d ij should reproduce δ ij up to transofrmation f . Since f is increasing, we have that δ ij ≤ δ kl implies that f (δ ij ) ≤ f (δ kl ). Introducing f in the cost function of MDS somehow relaxes distance preservation: instead of requiring strict one-to-one equalities d ij = δ ij and d ik = δ ik , nonmetric MDS merely enforces ordinal relationships in triplets (i, j, k), namely,

d ij = f (δ ij ), d ik = f (δ ik ), together with δ ij ≤ δ ik imply that d ij ≤ d ik .
In nonmetric MDS, both X and f need to be optimized and identified. In particular, f is obtained with a technique called isotonic regression.

Neural approaches of DR with data reconstruction

In this subsection, we present DR methods based on artificial neural networks.

This type of methods are inspired by biological neural networks. The main principles of neural networks are detailed as following.

An artificial neural network (ANN) is a set of nodes, which share connections. In classic ANN, a node is constituted with two functions:

• A combination function, which calculates a value with the input nodes and the weights of the connections.

• An activation function, which is applied to the value calculated with the combination function and transfered to the following layer of nodes.

A learning phase identifies the optimal weights of the synaptic connections with optimization methods. In a first part, we detail the self-organizing maps (Section 5.3.1). Secondly, we present the auto-associative neural networks (Section 5.3.2).

Self-organizing maps

Linear projection with PCA or classical MDS consists in fitting a plane through the data cloud and then projecting data onto it. If data is distributed near a nonlinear manifold, a possible extension of PCA consists in replacing the plane with a nonlinear structure. In the case of the SOM, the plane is usually restricted to a rectangular, articulated grid of neurons, which can deform more or less freely, like in Fig. 6. Each neuron has low-dimensional (LD) coordinates as a grid node

(a) (b)
Figure 6: Illustration of the concept of self-organizing maps. Figure 6a depicts a three-dimensional space in which a rectangular self-organizing map is folded according to a curved two-dimensional manifold. Figure 6b presents the corresponding low-dimensional grid.

and high-dimensional (HD) coordinates in the data space. The low-dimensional space is given by grid coordinates G = [g i ] 1≤i≤N and distances d(g i , g j ). The user fix them a priori, typically on a regular rectangular grid. In the high-dimensional space, the grid coordinates are denoted by Γ = [γ i ] 1≤i≤N ; they are adjusted as follows. The available data set is traversed several times, each traversal being called an epoch. At each epoch, data vectors ξ k are presented one by one. For each of those, the index of the closest grid node is denoted by j = arg min i ξ k -γ i 2 . The HD coordinates of all grid nodes are then updated according to

γ i ← γ i + αK d(gi,gj ) λ (ξ k -γ i )
, where K is a decreasing function from IR + to IR + . Function K can be a flipped step function (one below some threshold λ, zero beyond) or a Gaussian function with standard deviation λ. Parameter λ can be interpreted as a neighborhood radius in the grid space. Parameter α is a so-called learning rate and determines the neural plasticity, namely, the strength of the influence of any data vector ξ k on a neuron. In other words, high α gives low inertia to neurons, which are then easily attracted by data vectors, in a kind of Brownian motion in time. Parameter λ controls the influence of the LD grid. High λ means that the winning neuron together with its neighbors will be attracted by the data points. Reducing λ decreases this effect and neurons then become loosely connected and almost independent of each other, nearly overlooking their respective location in the grid. In practice, both α and λ decrease in time from epoch to epoch. The schedule of their respective decrease must be carefully adjusted for strong attractions and cohesions to occur in the beginning, to properly fit and unfold the SOM in the data space, while weaker attractions and cohesions allow for finer, smaller-scale adjustment later on.

Self-organizing maps have been developed mostly in a heuristic way and only a few attempts have been made to formulate their underlying cost function [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF]. The main intuition behind SOMs is to ensure consistency between neighborhoods in both spaces (LD grid and HD data), through the joint attraction of the closest neuron to a data vector and its neighbors. Such cohesion makes that grid neighbors of a neuron are also very likely to be its neighbors in the data space. In that sense, a SOM allows for visualizing proximities among data. The use of neurons, typically in smaller number than the data points, also relates SOMs to vector quantization and methods like competitive learning [START_REF] Rumelhart | Feature discovery by competitive learning[END_REF]. Neurons will distribute themselves within the data cloud so as to best represent them.

The use of neurons with a preset position in a grid makes the SOM an outlier in the field of dimensionality reduction. Instead of building a LD scatterplot of N points, representative of the N HD data vectors, the SOM goes the other way around, by fitting a predefined LD structure within the HD data distribution. From a visualization standpoint, the SOM proceeds by fetching from the HD data space some local properties of interest and displays it on the grid representation. Such local property can be a class or cluster label encountered near the neuron, the local density, or simply some component of ξ k . For an illustration of SOM application on the well known iris dataset, see [VHA + 99].

Auto-associative neural networks

Auto-associative networks [Kra91, Oja91, UNN91, DC93, HS06] are another neural approach to the problem of dimensionality reduction and, like the SOMs, they can be seen as a nonlinear extension of PCA. By fitting a plane through the data cloud, PCA seeks variance preservation and, hence, minimal loss of variance after projection on the plane. This also means minimal residues between the data points and their projected counterparts. PCA identifies the linear transformation (and its pseudo-inverse) that minimize these residues. Auto-associative networks replace the linear transformation with a nonlinear one. More specifically, artificial neural networks can approximate any function provided their architecture be adapted to the problem complexity. Complicated functions typically require more neurons and layers [START_REF] Bengio | Learning deep architectures for AI[END_REF][START_REF] Lecun | Deep learning[END_REF]. Figure 7 show a typical auto-associative network, with a succession of layers, among which the central narrower one forms a bottleneck and gives thus the whole network an hourglass shape. Auto-association means that input data (ξ) must be closely reproduced in the outputs ( ξ), with minimal residues, with the difficulty of passing through a central layer (x) with fewer dimensions than the inputs and outputs. Data flowing through the network must therefore be first compactly encoded and then decoded back to its initial form, hence the alternative name auto-encoder. Dots in the figure represent the neurons. They combine scalar products (like the linear projection in PCA), bias terms (extending mean subtraction in PCA), and a nonlinear function (unlike PCA). Nonlinearities are typically sigmoid functions or rectified linear units (ReLU), namely piecewise linear functions. The difficulty of using neural networks in an auto-associative setting is that it requires a 'deep' network, that is, with many intermediate layers, so that both halves, the encoder and decoder, are complex enough to approximate nonlinear functions. Parameters in those layers (synaptic weights) are adjusted by training the network, which is achieved by presenting small batches of input and output pairs. The reconstruction error associated with these batches is then computed, as well as its gradient with respect to the synaptic weights. Whereas inputs flow forward along the network to compute the corresponding outputs, gradient information goes back from output layers to input layers, in a process called gradient back-propagation, where synaptic weight are tuned to improve auto-association. Depth makes networks more difficult to train. In particular, gradient back-propagation becomes critical because the derivative chain-rule then involves many factors, with significant probability of observing small factors or, conversely, very large ones, leading to the dual phenomena of gradient vanishing and gradient explosion. New learning techniques and network architectures can partly mitigate these risks. Nevertheless, to date, auto-encoders are still outperformed by DR methods like stochastic neighbour embedding, at least in visualisation tasks where the representation dimension must be very low. Auto-encoders keep the advantage of being fully parametric methods, in both ways (dimension reduction and expansion back). Also, proximity is entailed in both the reconstruction error and neural network smoothness (similar HD inputs yield similar LD representations, like with the topographic mapping in SOMs).

Spectral methods

Kernel PCA

Until here, various generalisations of PCA take inspiration from geometry (scalar product or distance preservation in MDS, articulated grid instead of linear subspace in SOMs) or function approximation (neural network as universal approximators instead of linear projections). Kernel PCA extends PCA from a more theoretical standpoint, by generalizing the notion of scalar product. Strictly speaking, kernel PCA [SSM98, Wil01] is thus a nonlinear extension of classical MDS, the dual of PCA, acting on the Gram matrix of scalar products instead of the covariance matrix. Intuitively, kernel PCA replaces the usual scalar products across all pairs of raw data in the Gram matrix with a scalar product in a so-called feature space, that is, after a nonlinear transformation of data. The surprising trick here is that there is no need to carry out these two steps (mapping to the feature space, then scalar product computation) explicitly. Instead, a theoretical result allows one to compute the new Gram matrix at once, starting from the regular one, without specifying the nonlinear mapping to the feature space.

Function k from IR M × IR M to IR is said to be a kernel if it is symmetric (k(a, b) = k(b, a)) and positive semidefinite, that is, fulfilling Mercer's condition:

g(a)k(a, b)g(b) da db ≥ 0 , ( 4
)
where g is any square-integrable function. Then, kernel k has the following property:

k(ξ i , ξ j ) = φ(ξ i ) T φ(ξ j ) , (5)
where φ : IR M → F, ξ → φ(ξ) is an unknown mapping from the initial data space to some feature space F. The unknown mapped coordinates in Φ = [φ(ξ i )] 1≤i≤N are involved in matrix K = [k(ξ i , ξ j )] 1≤i,j≤N = Φ T Φ, which is then a symmetric and positive semidefinite Gram matrix by construction. Here, the specificity is that K is built directly, using kernel k, without computing the scalar products φ(ξ i ) T φ(ξ j ). In other words, we use the left-hand side of the kernel definition, without bothering about the right-hand side. The choice of k, beforehand, implicitly fixes φ, whereas the usual approach would have been to select φ first and then get the kernel afterwards.

As soon as a Gram matrix is available, classical MDS can be applied in order to get a Pdimensional projection, from the unknown feature space to some low-dimensional linear subspace. To carry out MDS, data must be centered. The apparent difficulty here is that centering must be performed on Φ in F, which is unknown. At this point, the magic of double centering enters into play. Computing C T KC, where C = I -1 N 11 T is the centering matrix, allows us to center Φ without knowing it, by acting on K instead of Φ. Eventually, eigenvalue decomposition C T KC = UΛU T yields Φ = Λ 1/2 U T . Keeping then only the eigenvectors associated with the largest P eigenvalues gives us a projection from the feature space F to a P -dimensional representation space. Although this projection is linear, the use of kernel function k induces a nonlinear mapping φ. Kernel PCA can then be seen as two-step procedure, with first a nonlinear mapping from the data space to the feature space and then a linear projection from the feature space to the LD representation.

In practice, the central question raised but unanswered by kernel PCA is the choice of a particular kernel. Theory determines which kernels are admissible, like polynomial functions of the raw scalar product or a radial Gaussian function, but it does not state which ones are useful in a specific task. Hence, the main theoretical advantage of kernel methods, namely, the fact that they induce a nonlinear mapping implicitly, can also be their main drawback, since the properties or behaviour of the induced mapping are difficult to investigate. For instance, some kernels induce a mapping towards an infinite-dimensional space, which turns out to be counter-productive in the context of dimensionality reduction, whereas it is useful in clustering or classification. This is the reason why kernel PCA is hardly used as an effective DR method. On the other hand, it has provided a theoretical background for the development of other spectral DR methods, where the kernel is purposefully chosen or even data-driven. Isomap [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] can be cast within this framework. Instead of considering Euclidean distances and the associated usual scalar product, it relies on the concept of geodesic distances, which somehow can be interpreted as a manifold-aware metric. The geodesic distance is measured not as the crow flies, but rather following the underlying manifold, as the ant would crawl. In practice, the geodesic distances are approximated with shortest paths in a graph of K-ary neighbourhoods. As a consequence, this is no longer a true kernel, the method works fine with minor workarounds. In contrast, the kernel property is carefully preserved in a closely related method, namely, maximum variance embedding (MVU, formerly known as semidefinite embedding, or SDE) [START_REF] Weinberger | Learning a kernel matrix for nonlinear dimensionality reduction[END_REF][START_REF] Weinberger | Unsupervised learning of image manifolds by semidefinite programming[END_REF]. In this method, the kernel is not fixed arbitrarily in advance, nor precomputed from a graph as in Isomap, it gets continuously updated and optimised under constraints with a technique called semidefinite embedding. This technique updates the Gram matrix and ensures that it remains valid. All these DR methods end by applying CMDS to a modified, 'kernelized' Gram matrix, whose dominant eigenvectors yield an embedding.

Laplacian eigenmaps

Proximity preservation, spectral embedding, and kernelised Gram matrix are the key elements of another DR method, called Laplacian eigenmaps. The name Laplacian comes from a matrix operator in the domain of graph theory. The idea is that edges or connections in a graph indicate vertices or data that are similar, whereas non-connected items would be more dissimilar. Following the intuition of previous DR methods, one might want to find a low-dimensional representation of the graph where similar, connected vertices are close to each other, whereas non-connected ones are lying farther away.

Formally, let us consider a non-directed graph, with N vertices, and connections among those. The connections can be encoded in an adjacency matrix

A = [a ij ] 1≤i,j≤N where a ij = 1 if j and i are connected 0 otherwise . ( 6
)
Usual ways to define the adjacency matrix when starting from data vectors ξ i in Ξ are to connect either the closest K neighbours of ξ i or all neighbours lying within an ball around ξ i .

If the graph edges are weighted, then a weight matrix can be defined, too. For instance, we could have

W = A or W = [w ij ] 1≤i,j≤N where w ij = exp(-δ 2 ij /(2σ 2 )) if j and i are connected 0 otherwise , ( 7
)
where δ ij is the distance between data vectors ξ i and ξ j , and σ is some bandwidth. The value of σ should be such that the Gaussian is already close to zero for the farthest connected neighbours. Having a weight matrix for the neighbours, an embedding X of Ξ, in which neighbours are represented close to each other, can be found by minimizing

E(X; W) = N i,j=1 w ij x i -x j 2 2 . ( 8
)
A trivial solution to this problem is to map all data to a single point, that is, to have x i = x j for all i and j, regardless whether they are connected or not. To take into account the fact that non-neighours should be represented separately, let us assume that X is centered (X1 = 0) and constrain its covariance matrix 1 N XX T to be diagonal and full rank. This way, dimensions in X are decorrelated with nonzero variance. In other words, some spread of the embedding gets guaranteed.

Another observation is that the above problem can be reformulated using the Laplacian matrix, defined as L = D -W, with D = diag(W1). The Laplacian matrix allows us to write

E(X; W) = N i,j=1 w ij x i -x j 2 2 = Tr(X T LX) . (9)
Then, the solution that minimizes the weighted distances between neighbours while preserving some spread is arg min X E(X; W) such that X1 = 0 and X T DX = I. The advantage of reformulating the problem in this way is that the Laplacian matrix is square, symmetric, and positive semidefinite.

Its eigenvalue decomposition can be written L = UΛU T . The solution X = UD -1/2 fulfills the constraint, that is, XDX T = UD -1/2 DD -1/2 U T = I, since D is diagonal and U is orthogonal (UU T = U T U = I). So far the solution gathers all N eigenvectors of L. To get a P -dimensional embedding, and in order to minimize the cost function, the trailing P eigenvectors must be kept, namely, those with the smallest nonzero eigenvalues. Notice that L is rank-deficient with at least one zero eigenvalue and a trivial uniform eigenvector. Taking the trailing eigenvectors instead of the dominant ones is major difference with kernel PCA. Also, the considered Laplacian matrix is sparse (assuming only a few edges per vertex in the graph), in contrast to the dense Gram matrices used in CMDS, and thus in KPCA, Isomap, MVU, etc. Nevertheless, a form of duality exists between these two categories of methods. In particular, the pseudo-inverse of the Laplacian (UΛ + U T , where λ + ii = 1/λ ii if λ ii = 0) can be considered to be a Gram matrix. The associated distance is the commute-time distance [START_REF] Saerens | The principal components analysis of a graph, and its relationships to spectral clustering[END_REF], involving random walks in the neighbourhood graph, which can be seen as a 'kernelised' distance like the geodesic distance in Isomap.

Other methods using sparse matrices are, for instance, locally linear embedding (LLE) [RS00] and Hessian LLE [START_REF] Donoho | Hessian eigenmaps: New locally linear techniques for high-dimensional data[END_REF]. Laplacian eigenmaps is also closely related to spectral clustering [NJW02, BVP + 03], a method that uses the Laplacian matrix mostly in the same way to better separate clusters by mapping them to a feature space before projecting them onto a lower-dimensional space. Whereas spectral clustering is well studied and quite successful in practice, spectral embedding as a DR tool suffers from the same limitation as kernel PCA, namely, the a priori mapping induced by the graph Laplacian can actually increase data dimensionality, instead of reducing it, which is detrimental to DR but beneficial to clustering.

Probabilistic neighbourhood preservation

As described in Section 5.2, many DR methods consider that the proximities between data points can be well captured by measuring their distances. As a consequence, the computation of lowdimensional coordinates of the data is naturally driven by a distance preservation criterion. The low-dimensional space should indeed well render the high-dimensional proximities between the data, as exemplified by Figure 2.

Many recent DR developments however question the last assumption that distances well summarize proximity relationships [START_REF] Lee | Shift-invariant similarities circumvent distance concentration in stochastic neighbor embedding and variants[END_REF]. Distances in spaces with varying dimensions indeed behave differently. Due to the norm concentration effect induced by the curse of dimensionality [START_REF] Lee | Nonlinear dimensionality reduction[END_REF], distances in high-dimensional spaces tend to concentrate toward much higher values than in lowdimensional ones. This misleads criteria aiming to exactly preserve distances between spaces with different dimensions.

The last observations motivated the development of proximity indicators being more robust when computed in different spaces. In particular, these indicators should have the ability to be invariant with respect to distance shifts. This would allow them to at least partly alleviate the undesirable consequences of the norm concentration.

In this context, the Stochastic Neighbor Embedding (SNE) method [START_REF] Hinton | Stochastic neighbor embedding[END_REF] constituted a real paradigm shift in DR. It is based on the definition of similarities between pairs of data points in both HD and LD spaces. These similarities more faithfully characterize the proximity relationships

HDS LDS K ξ i ξ j σ ij x i x j s ij

Minimize divergence

Figure 8: Illustration of similarity-based dimensionality reduction. In both high-and lowdimensional spaces (HDS and LDS), the data points are depicted using black dots. The HDS and LDS similarities between each pair of points are computed. This defines fields of influences around each data point in both spaces. The ones around ξ i and x i are represented by the grays shades with decreasing intensities as the distances to ξ i and x i increase. This leads to soft Kary neighbourhoods around each datum, as represented around ξ i . As the HDS configuration is fixed, the low-dimensional positions of the data points are computed by minimizing the divergence between the HDS and LDS similarities.

among the data, precisely thanks to their invariance with respect to distance shifts. From the pairwise distances, the latter similarities are defined for i = j in the high-and low-dimensional spaces as

σ ij = exp(-π i δ 2 ij /2) k,k =i exp(-π i δ 2 ik /2) and s ij = exp(-p i d 2 ij /2) k,k =i exp(-p i d 2 ik /2) . (10)
Since proximities are only relevant to evaluate between different data items, σ ii and s ii are set to 0. The precisions of the normalized Gaussian functions centered on ξ i and x i are denoted by π i and p i . σ ij (resp. s ij ) can be interpreted as the probability of ξ j (resp. x j ) to lie in a soft Gaussian neighbourhood around ξ i (resp. x i ). The precisions of the Gaussian functions allow to characterize these neighbourhoods in terms of their sizes. Smaller (resp. larger) precisions indeed flatten (resp. narrow) the Gaussian functions, defining larger (resp. smaller) neighbourhoods as more (resp. less) data points have a non-negligible probability to belong to them. Hence the precision value and the soft neighbourhood size uniquely determine each other. As it is more natural to reason in terms of sizes than precisions, one usually fix the typical number of neighbors around each data point, which is termed as the perplexity. The precisions of the Gaussian functions are then deduced to equalize the soft Gaussian neighbourhoods sizes and the perplexity.

More formally, the perplexity K i of the soft neighbourhood around ξ i is defined as

K i = exp(H i ) and H i = - N j=1 σ ij ln σ ij , (11)
where H i denotes the entropy of the distribution of σ ij . As intuitively explained above, the neighbourhood width, the size K i and the entropy H i increase as the precision π i decreases. In all (single-scale) variants of SNE, a unique perplexity value K is provided by the user in order to adjust each HD precision π i by solving (11) with

K i = K for 1 ≤ i ≤ N .
In order to identify the LD precisions p i , the same procedure can not be applied as the coordinates x i are unknown. To circumvent this issue, SNE fixes p i = 1 for all i.

Having defined the similarities σ ij and s ij , X can be computed in order to match σ ij and s ij for all 1 ≤ i, j ≤ N . As these similarities sum up to one ( j σ ij = j s ij = 1), σ i = [σ ij ] 1≤j≤N and s i = [s ij ] 1≤j≤N define discrete probability distributions. Divergences can therefore evaluate their mismatch. This process is illustrated by Figure 8.

The Kullback-Leibler divergence is used in SNE. It can be developed as

D KL (σ i s i ) = j σ ij ln(σ ij /s ij ) . (12)
The cost function of SNE [START_REF] Hinton | Stochastic neighbor embedding[END_REF] then consists in the sum of (12) for each i:

E(X; Ξ, p, π) = i D KL (σ i s i ) , (13)
where p = [p i ] 1≤i≤N and π = [π i ] 1≤i≤N . By optimizing this cost function, one can find the low-dimensional coordinates of the data points. The SNE method is hence fully defined. Over the years, enhancements of the basic SNE method were proposed in order to tackle problems with increasing difficulties. One can first cite the t-SNE method [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF]. It modifies SNE in two ways:

• The HD similarities are made symmetric in order to simplify the expression of the cost function gradient:

σ ij ← σ ij + σ ji 2N . (14)
• The LD similarities are defined using a Student t-distribution with one degree of freedom:

s ij = 1 + d 2 ij -1 k =l (1 + d 2 kl ) -1 . (15)
At the exception of these two modifications, t-SNE remains identical to SNE, in particular concerning its cost function. Although symmetrizing the HD similarities is not clear to have a significant effect on the performances [LRB + 13], modifying the LD similarities is crucial. Using a heavy-tailed distribution in the LD space has been shown to allow to partly alleviate the so-called crowding problem. When reducing the dimensionality of very high-dimensional spaces, one can not hope to preserve exact neighborhood relationships in the LD space as in the HD one as there is much more "place" around each datum in high dimensions. Hence data points lying at a moderate distance from each other must be pushed much farther away in the LD space in order to best preserve close neighborhoods. The use of a heavy-tailed distribution in the LD space therefore allows data points at a moderate distance in the HD space to lie far apart in the LD one while preserving a LD similarity close to their HD one. This avoids misleading the cost function. On the other side, the Student t-distribution in the LD space induces an exponential transformation between the HD and LD distances [LRB + 13], which prevents t-SNE to yield isometric embeddings of linear manifolds.

Another variant of SNE is the Neighborhood Retrieval Visualizer (NeRV) [VPN + 10]. It uses the exact same HD and LD similarities as SNE, at the exception that it fixes p i = π i . On the other hand, it modifies its cost function by blending two dual KL divergences in a so-called type 1 mixture [START_REF] Cichocki | Families of alpha-beta-and gamma-divergences: Flexible and robust measures of similarities[END_REF]:

D κ KLt1 (σ i s i ) = (1 -κ)D KL (σ i s i ) + κD KL (s i σ i ) , (16)
where parameter 0 ≤ κ ≤ 1 trades off both terms. This idea of mixing dual objective functions was previously investigated in the distance preservation paradigm [START_REF] Venna | Local multidimensional scaling with controlled tradeoff between trustworthiness and continuity[END_REF]. NeRV cost function is then written as E(X; Ξ, p, π, κ)

= i D κ KLt1 (σ i s i ) . (17)
This divergence mixture can be interpreted as a way to balance the precisions and recalls of the retrieved LD neighborhoods around the data points [VPN + 10].

Introducing the notation z i = κσ i + (1 -κ)s i , another KL divergences combination can be defined as

D κ KLt2 (σ i s i ) = κD KL (σ i z i ) + (1 -κ)D KL (s i z i ) , (18)
D 1/2
KLt2 (σ i s i ) is termed as the type 2 symmetric KL divergence [START_REF] Cichocki | Families of alpha-beta-and gamma-divergences: Flexible and robust measures of similarities[END_REF], or symmetric Jensen-Shannon divergence [START_REF] Burbea | On the convexity of some divergence measures based on entropy functions[END_REF][START_REF] Lin | Divergence measures based on the shannon entropy[END_REF]. The cost function becomes

E(X; Ξ, p, π, κ) = 1 κ(1 -κ) i D κ KLt2 (σ i s i ) , (19)
which is effectively used in the Jensen-Shannon embedding (JSE, or 'Jessie') method [LRB + 13], relying on the same HD and LD similarities as SNE. Other divergences can be defined and are discussed in [START_REF] Bunte | Stochastic neighbor embedding (sne) for dimension reduction and visualization using arbitrary divergences[END_REF].

As the obtained low-dimensional coordinates result from the optimization of the cost function, the success of similarity-based methods is crucially dependent on this last step. Yet cost functions of SNE-like methods are known to be difficult to optimize [HR03, vdMH08, LRB + 13, LPOnV15]. They indeed most probably present many local minima, in particular when K is small compared to N since the localness of the similarities is then amplified. The squared distances in the normalised Gaussian functions furthermore poorly scale the gradient. These observations suggested to rely on optimization schemes benefiting from second order information about the cost function [START_REF] Yang | Multiplicative updates for t-sne[END_REF][START_REF] Van Der Maaten | Fast optimization for t-sne[END_REF]. Limited-memory BFGS [START_REF] Ghahramani | Updating quasi-Newton matrices with limited storage[END_REF] therefore replaces gradient-descent in NeRV [VPN + 10] and JSE [LRB + 13]. Moreover by modifying the values of the HD similarities in the few first iterations, either directly [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF] or with an increased value of K [VPN + 10, LRB + 13], one can typically lessen the influence of local minima.

Until this point, all the detailed SNE variants considered a single, fixed perplexity value K in order to define the HD similarities. As detailed above, the latter perplexity implicitly determine the soft neighborhood size around each datum. The relationships between pairs of data points inside these neighborhoods will hence be considered while designing the low-dimensional embedding, but not the ones outside these neighborhoods. As K is usually small with respect to N , SNE-like methods thus best preserve close neighborhoods around each datum. However they usually do not capture the global structures of data sets as they only focus on local interactions. Large-scale neighborhoods are hence often poorly reproduced [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF]. Using only a single perplexity value K indeed enforces to choose between preserving local or global structures. As one typically wish to retrieve both, combining several single-scale similarities, defined using different perplexities, into multi-scale ones could allow to focus both on local and global interactions [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF].

With this idea in mind, an index h can be introduced in the single-scale similarities:

σ hij = exp(-π hi δ ij /2) k,k =i exp(-π hi δ ik /2)
and

s hij = exp(-p hi d ij /2) k,k =i exp(-p hi d ik /2) . ( 20
)
The HD and LD precisions of the ith soft neighbourhoods on scale h are denoted by π hi and p hi . Again, σ hii and s hii are set to 0. Multi-scale similarities can then be defined as

σ ij = 1 L Lmax h=Lmin σ hij and s ij = 1 L Lmax h=Lmin s hij , (21)
with 1 ≤ L min ≤ h ≤ L max and where the considered number of scales is L = L max -L min + 1.

Obviously if

N j=1 σ hij = 1 for all h and i, then N j=1 σ ij = 1 for all i. It should also be noted that the multi-scale similarities reduces to single-scale ones if L min = L max .

Multi-scale SNE, t-SNE, NeRV and JSE methods can hence be defined [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF][START_REF] De Bodt | Perplexity-free t-SNE and twice Student tt-SNE[END_REF], using multi-scale similarities instead of single-scale ones in their respective cost functions. In order to determine the precisions π hi at the various scales, perplexity values distributed over the whole range from 1 to N can be fixed without any user intervention. One could for example start from a small one such as K = 2, and increase it by powers of two:

K h = 2 h-1 K with 1 ≤ L min ≤ h ≤ L max = log 2 (N/K ) ,
• denoting the rounding operator. Defining L max in this way allows to consider large neighbourhoods with respect to N while preventing σ hij to become nearly uniform. Optionally, setting L min > 1 allows to more focus on large neighbourhoods.

On the other side, the determination of p hi is related to the estimation of the correlation dimension [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF]. Finally multi-scale optimisation can be carried out by introducing the smaller scales one at the time during the optimisation. This procedure allows to ease the multiscale methods tuning [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF].

Quality assessment

Dimensionality reduction is by nature an unsupervised learning task, namely, there is no ground truth in most practical cases and, therefore, no obvious performance assessment. It is the user's responsiblity to devise a performance criterion, which is often directly formalised into the objective function of the methods like those detailed above. This situation, however, makes it difficult to compare methods to each other.

For example, one could rely on variance preservation to assess DR quality, but doing so would bias evaluation favorably for PCA and CMDS, although nonlinear DR methods are designed and expected to outperform them in many cases. Similarly, stress functions have sometimes been deemed generic enough to be adopted as a universal quality measure. Again, this choice would not be fair as stress functions favor MDS-like methods that precisely optimize them. Moreover, distance preservation is known to be suboptimal in many DR problems where distances heavily concentrate in high-dimensional data spaces. Measures of data reconstruction, like those used in PCA and AEs, could be used as well and appear again as sufficiently generic. In this case, however, they apply only to parametric and inversible DR.

Yet another possibility to assess DR quality is to rely on some external ground truth that is not directly involved in the DR problem. The most typical example occurs when multivariate data is primarily collected to solve a supervised task. Then, data consists of both of inputs ξ i and associated desired outputs y i , like class labels or some regression output. The main issue with this approach is that it is indirect and requires such additional labels or variables to be systematically available even if they are not strictly necessary for DR itself, and possibly not relevant. For instance, class labels would typically be used to assess DR quality by measuring the classification performance in the embedding space. Better classification results would then mean better embedding. Notice also that DR quality is assessed only near class boundaries, since embedding errors within the bulk of a class would pass undetected.

With the recent cornucopian development of many new DR methods, sometimes very complicated, several attempts have been made to assess the quality of DR results in a generic and methodindependent way. Instead of variance or distance preservation, these tools focus on neighborhood preservation. Intuitively, this means that distances must not be preserved in value, but just in ranking. In other words, let us consider three points ξ i , ξ j , and ξ k and their embeddings x i , x j , and x k ; if ξ i is taken as vantage point, with δ ij < δ ik , then distance preservation aims at δ ij = d ij and δ ik = d ik , whereas neighborhood preservation would relax this into δ

ij ≤ δ ik ⇒ d ij ≤ d ik .
The main motivation is that distance rankings are expected to be relatively immune to norm concentration, contrarily to raw distances.

Formally, distance rankings call for the definition of ranks and K-ary neighborhoods. Points can be ranked in rising order of their distance with respect to a given point of interest. The rank of ξ j with respect to ξ i in the HD space is written as ρ ij = |{k : δ ik < δ ij or (δ ik = δ ij and 1 ≤ k < j ≤ N )}|, where |A| denotes the cardinality of set A. Similarly, the rank of x j with respect to x i in the LD space is r ij = |{k : d ik < d ij or (d ik = d ij and 1 ≤ k < j ≤ N )}|. Notice that reflexive ranks are set to zero (ρ ii = r ii = 0) and that ranks are unique, i.e., there are no equalvalued ranks: ρ ij = ρ ik for k = j, even if δ ij = δ ik . This means that nonreflexive ranks belong to {1, . . . , N -1}. The non-reflexive K-ary neighbourhoods of ξ i and x i are the sets defined by K i = {j : 1 ≤ ρ ij ≤ K} and n K i = {j : 1 ≤ r ij ≤ K}, respectively. With this notation, a measure of K-ary neighbourhood preservation [START_REF] Lee | Quality assessment of dimensionality reduction: Rankbased criteria[END_REF] can be written as

Q NX (K) = 1 KN N i=1 |ν K i ∩ n K i | , (22)
which is the average proportion of preserved K-ary neighbours, for 1 ≤ K ≤ N -1. It varies between 0 (empty intersection) and 1 (perfect agreement). Knowing that random coordinates in X lead on average to

Q NX (K) ≈ K/(N -1) [CB09], the useful range of Q NX (K) is N -1 -K, which depends on K.
Therefore, in order to fairly compare or combine values of Q NX (K) for different neighbourhood sizes, the criterion can be rescaled, like in [Kun07, LRB + 13], to get

R NX (K) = (N -1)Q NX (K) -K N -1 -K , (23)
for 1 ≤ K ≤ N -2. This modified criterion indicates the improvement over a random embedding and has the same useful range between 0 (random) and 1 (perfect) for all K.

Examples of DR results

In order to illustrate the results of various DR methods and their quality assessment, Brendan Frey's faces are used. This database is a collection of images of a single person's face, with various facial expressions, extracted as consecutive frames of a video clip. The data base includes N = 1965 images, with size 28 times 20 pixels, some of them being illustrated in Fig. 9. The images are simply vectorized, namely, the 28 rows of 20 pixels each are concatenated to form long 560-dimensional vectors ξ i . Starting from these high-dimensional vectors, dimensionality reduction can determine low-dimensional representations. If the target dimensionality is two, then the representations can be displayed as scatter plots, where each point corresponds to a facial expression. In the scatter plots, the same color code as in Fig. 9 is used.

Figure 10 shows the embeddings of Frey's faces with various DR methods. The considered ones are:

• A linear projection with classical MDS (CMDS) [START_REF] Young | Discussion of a set of points in terms of their mutual distances[END_REF], which is equivalent to PCA [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF].

• Variants of stress-based MDS, namely non-metric (NMDS) [She62, Kru64], Sammon's nonlinear mapping (NLM) [START_REF] Sammon | A nonlinear mapping algorithm for data structure analysis[END_REF], and curvilinear component analysis (CCA) [START_REF] Demartines | Curvilinear component analysis: A self-organizing neural network for nonlinear mapping of data sets[END_REF].

• Spectral embedding with Laplacian eigenmaps (LE).

• Probabilistic neighbour embedding, with SNE [START_REF] Hinton | Stochastic neighbor embedding[END_REF], NeRV [VPN + 10], JSE [LRB + 13], their multi-scale counterpart [START_REF] Lee | Multi-scale similarities in stochastic neighbour embedding: Reducing dimensionality while preserving both local and global structure[END_REF], and t-SNE [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF].

This short list covers many types of methods, from linear to nonlinear, from old to quite recent, and with various approaches. The focus is set on neighbour embedding techniques, which are currently the state of the art in the domain. For these methods, the perplexity parameter is set to 32 by default, except for the multi-scale extensions, which consider automatically several neighbourhood sizes. Figure 11 shows the quality assessment of the scatter plots in Fig. 10. In this plot, each colored curve corresponds to an embedding. The horizontal axis is a log-scale of neighbourhood size K. The vertical axis indicates R NX (K), that is, the relative improvement over an arbitrary embedding with random coordinates. A null value means no improvement, whereas 100% is synonymous of perfect neighbourhood preservation in the low-dimensional representation. The areas under the curves, providing global scores over all K, with more focus on the small neighbourhoods due to the log-scale, are reported in the legend, next to the methods' names.

A method like PCA (or CMDS) provides a linear projection of data, which typically reveals the global shape of the point cloud. The corresponding curve culminates indeed for rather large neighbourhood sizes. Small neighbourhoods are not well reproduced due to false neighbours, that is, intruders caused by the inevitable squashing entailed by linear projections. Showing a similar behaviour, nonmetric MDS performs even better, with the best rendition of large neighbourhoods. Another stress-based MDS method like Sammon's NLM, does not perform as well as NMDS, but is much simpler; the embedding is more cluttered and indistinctly circular. With some added complexity, CCA improves on small neighbourhoods; also some clusters become discernible. As a method closely related to spectral clustering through the use of the graph Laplacian matrix, Laplacian eigenmaps clearly separates the few main clusters, at the expense of collapsing them. Variants of SNE, like t-SNE, NeRV, JSE, and their multi-scale extensions, through the use of normalised similarities, can better cope with high-dimensional data than stress-based MDS, which tries to match distances directly, disregarding norm concentration effects [START_REF] François | The concentration of fractional distances[END_REF][START_REF] Lee | Shift-invariant similarities circumvent distance concentration in stochastic neighbor embedding and variants[END_REF]. Due to slightly discrepant similarity definitions (Student versus Gaussian), t-SNE tends to over-emphasize cluster separation, compared to SNE. With respect to neighbourhood preservation, SNE outperforms CCA, and it provides an embedding where clusters are less packed. NeRV behaves even better than SNE, but does not match the performance of t-SNE for very small neighbourhoods. Another SNE variant, JSE performs similarly as t-SNE, being slightly more faithful overall, at the cost of small neighbourhoods. Eventually, considering neighbourhoods over a large range of sizes allows multi-scale variants of SNE (Ms.SNE, Ms.NeRV, Ms.JSE) to perform the best, almost identically. Clusters are clearly visible, without magnifying their separation, leaving some room for intra-cluster details. The smallest neighbourhoods reach nearly 90%, whereas the largest ones remain above 50%.

These results show the improvements achieved over more than a century of research in the domain, both quantitatively and visually.

Alternative metrics and distances

Mahalanobis

Introduced by [C + 36], the Mahalanobis Distance (M D ), as for Euclidean distance, can be calculated in the original HD space and in the principal component (PC) space. M D is based on the correlation between variables by which different models can be identified and analyzed. This is a useful way to determine the similarity between a series of known and unknown data. It differs from the Euclidean distance in that it takes into account the variance and correlation of the data series. Thus, unlike the Euclidean distance where all the components of the vectors are treated independently and in the same way, the M D gives a smaller weight to the most dispersed components. In the case of signal analysis, and assuming that each component is a Gaussian random variable, this amounts to minimizing the influence of the most noisy components (those with the largest variance). The M D is often used for the detection of outliers in a dataset, or to determine the consistency of data provided by a sensor for example: this distance is calculated between the received data and those predicted by a model. In practice, the MD of a multivariate vector x = (x 1 , x 2 , x 3 , ..., x M ) T to a set of mean value vectors µ = (µ 1 , µ 2 , µ 3 , ..., µ M ) T and having a covariance matrix C X is defined by D M (x) = (x -µ) T C -1

x (x -µ). Finaly, it is important to note that the squared Mahalanobis distance is equal to the sum of squares of the scores of all non-zero standardised principal components (see [START_REF] Richard G Brereton | The mahalanobis distance and its relationship to principal component scores[END_REF] for more details).

Figure 12: Illustration of the concept of geodesic distance. The data points are distributed on a three-dimensional molted swiss roll manifold. A K-neighborhood graph is constructed on the data set, its edges being depicted as gray dotted lines linking data points. An edge can be weighted by the Euclidean distance between the two dots it binds. Some short-circuits across the manifold can be observed due to the fact that the data set is finite. The geodesic distance between the two squared data points is then computed as the length of the shortest path between them in the graph. This shortest path is depicted by solid black lines. On the other side, the dashed black line represents their Euclidean distance. One can observe that the geodesic distance more faithfully describes the relationship between data points as it tries to approximate the manifold on which they are distributed.

Geodesic distances

Used in the well known ISOMAP method [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], geodesic distance corresponds to a useful nonlinear metric for estimating the intrinsic geometry of a data manifold. Compared to traditional Euclidean distance (Figure 12), Geodesic distance consists in the sum of edge weights along the shortest path between two nodes in a K-neighborhood graph of the observed data (computed using Dijkstra's algorithm, for example). In the context of DR, this consists in estimating the shortest paths in the graph for all pairs of data points in high dimension and/or low dimension.

Commute-time distances

See Section 5.4.2.

Gower's similarity

The Gower's similarity is a similarity which comes from ecologica analysis [START_REF] John | A general coefficient of similarity and some of its properties[END_REF]. This similarity is designed for heterogeneous data sets. We define an heterogeneous data set as a data set wherein the objects are described with categorical and numerical variables. Moreover, the Gower's similarity can deal with missing data. The Gower's similarity is calculated as following. Let two objects x i and x j , described by M Table 1: An example of heterogeneous data set

V 1 V 2 V 3 V 4 V 5 V 6 x 1 1,5 A 2,3 x 2 7 10 W x 3 A 8 9 W
variables, which can be categorical or numerical. The Gower's similarity S ij between x i and x j is:

S ij = M k=1 w ijk S ijk M k=1 w ijk (24)
In this formula, S ijk is a "local" similarity, calculated between x i and x j for the variable k, such as:

• if k is a categorical variable, then

S ijk = { 0 if x ik = x jk 1 if x ik = x jk (25)
• if k is a numerical variable, then

S ijk = 1 - x ik -x jk max(x k ) -min(x k ) (26)
In this formula, w ijk is a weight associated to the variable k for the pair (x i , x j ), such as:

• if x ik or x jk are missing data, then w ijk = 0

• else w ijk = 1

The Gower's similarity has values between 0 and 1: 1 indicates a maximal similarity. The Gower's similarity can be easily transformed into a dissimilarity: D ij = 1 -S ij . The Gower's similarity is an interesting similarity to manage heterogeneous data sets with missing data. But this similarity should be used carefully, due to these properties.

In the table 1, we propose an example. In this example, S 1,2 cannot be calculated: x 1 and x 2 have no variable in common, because of missing data. But, S 1,3 and S 2,3 can be calculated. Moreover, S 1,3 = 1, the maximal possible value for a Gower's similarity.

In conclusion on missing data, the Gower's similarity can be used on data set with missing data, but the user of this similarity have to ensure than missing data do not compromise the representation of the data set.

Another point which schould be discussed, before using Gower's similarity, is the weighting of pair (x i , x j ). As mentionned below, these weights are originally used to manage missing data. But, we assume that these weights can be used to improve the calculation of the similarity. Two approaches can be considered. On the one hand, the weights can be chosen by the user of the similarity, depending on his expertise on the data set. On the other hand, the weights can be calculated, optimized in order to improve the representation of similarities between vectors in a data set.

Concerning this second approach, some works offer algorithms that set the weights in order to optimize a clustering process [vdH16] . But, for the moment, the optimization of variable weights in a reduction dimension perspective is actually not adressed.

Conclusions

Dimensionality reduction has been studied for over a century. Different paradigms and techniques have been investigated to create a plethora of methods. The oldest one is PCA, based on a variance-preserving linear projection. Classical MDS seeks a projection that preserves pairwise scalar products. Both work by eigenvalue decomposition. Other variants of MDS yield nonlinear mapping with ad hoc optimisation of so-called stress functions that quantify distance preservation (or (dis)similarity ranks in nonmetric MDS). Various neural networks, like self-organising maps and auto-encoders, can achieve DR as well, also with ad hoc optimisation techniques. Spectral optimisation, initially associated with linear projection only, has raised interest again in combination with the kernel trick, leading to methods like kernel PCA, Isomap, Laplacian eigenmaps, or locally linear embedding. Eventually, the current state of the art relies on stochastic neighbourhood embedding (SNE, t-SNE), where distances are wrapped in similarities with invariance properties that allow mitigating distance concentration, which affects stress-based MDS.

In parallel to DR method development, new tools to assess DR quality have been devised, trying to capture to essence of modern, nonlinear DR, which focuses more and more on the preservation of local neighbourhoods from the HD data space to the LD representation space.

Dimensionality reduction remains a branch of non-supervised machine learning, used mainly for exploratory data analysis and visualisation. Quality assessment and the choice of the metric to evaluate actual proximities between data are still open questions.

Figure 3 Figure 3 :

 33 Figure 3: Illustration of PCA seen as a multivariate regression process.

Figure 7 :

 7 Figure7: Illustration of auto-associative neural networks. The network is very similar to Figure3except that instead of directly joining two layers, the straight lines meet at some black dots, which represents nonlinear activation functions of the linear combinations at their inputs. Hence auto-associative neural networks are similar to PCA but model the transformations e θ and d θ as nonlinear functions.

Figure 9 :Figure 10 :

 910 Figure 9: Brendan Frey's faces. Some examples of facial expression images from the database, as well as the color code used further in the embedding scatter plots. As consecutive frames from a video clip, the thumbnail images should be read from top to bottom and then from left to right. The shades of colors in the bottom give an indication of time, per column.

Figure 11 :

 11 Figure 11: Quality assessment for the embeddings of Brendan Frey's faces in Fig. 10. Quality is evaluated through K-ary neighbourhood preservation. The neighbourhood size K varies along the horizontal axis.

He is currently a UCLouvain Professor and a Senior Research Associate with the Belgian fund of scientific research (Fonds de crossing of artificial intelligence, machine learning, information visualisation, image processing,