{'fraction_non_alphanumeric': 0.037051320083061404, 'fraction_numerical': 0.006437258973598339, 'mean_word_length': 4.749786798567286, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On genetic programming representations and fitness functions for interpretable dimensionality reduction

Dimensionality reduction (DR) is an important technique for data exploration and knowledge discovery. However, most of the main DR methods are either linear (e.g., PCA), do not provide an explicit mapping between the original data and its lower-dimensional representation (e.g., MDS, t-SNE, isomap), or produce mappings that cannot be easily interpreted (e.g., kernel PCA, neural-based autoencoder). Recently, genetic programming (GP) has been used to evolve interpretable DR mappings in the form of symbolic expressions. There exists a number of ways in which GP can be used to this end and no study exists that performs a comparison. In this paper, we fill this gap by comparing existing GP methods as well as devising new ones. We evaluate our methods on several benchmark datasets based on predictive accuracy and on how well the original features can be reconstructed using the lower-dimensional representation only. Finally, we qualitatively assess the resulting expressions and their complexity. We find that various GP methods can be competitive with state-of-the-art DR algorithms and that they have the potential to produce interpretable DR mappings.


INTRODUCTION
Dimensionality reduction (DR) is a key instrument for knowledge discovery [8,11]. Informally, we can define DR (see Section 2.1 for a formal definition) to be the area of study concerned with finding lower dimensional representations of the original data, such that some meaningful properties of the original dataset are preserved. It is based on the fact that most real-world datasets only have artificially high dimensionality due to, for instance, interdependent or noisy features. However, current state-of-the-art DR algorithms act as a black-box when going from high to low dimensions (e.g., t-SNE [21], isomap [25], LLE [25]), and do not provide any means to inspect the resulting lower dimensions in terms of the original variables, i.e., they do not provide a functional mapping that describes how the original dimensions are compressed into the lower-dimensional (or also called latent) representation. Other methods, such as kernel PCA [37] and neural autoencoders, do provide a functional mapping between the original data and its lower dimensional representation, however these mappings are hardly interpretable. Furthermore, in the absence of meaningful labels or classes, it becomes even more challenging to assess and inspect whether a lower-dimensional representation produced by DR is truly meaningful. Finally, DR is often used for data exploration, in which case the ability to draw insights from the data is key. For these reasons, producing a functional mapping that is human-interpretable can be very important, as that allows to assess whether the compression the mapping performs may be reasonable and, e.g., safe to use in high-stakes applications of machine learning [36].
Recently, the use of genetic programming [14] (GP) has been an increasingly popular way to produce interpretable models. For instance, GP has been used to learn programs encoded by simple symbolic expression representing physical laws [4,9]. Furthermore, recent works have investigated the use of such symbolic expressions produced by GP as interpretable models [7,17,31,33]. In this work, we leverage the interpretability potential of GP and apply it to non-linear DR.
GP has been shown to be a promising avenue to produce interpretable lower-dimensional representations in several recent works [17,18,20]. However, more representations and evaluation schemes can be imagined, but no work has been done to provide a comparison of these different options. We attempt to fill this gap by designing new GP-based methods for DR and establish a robust evaluation protocol across two objectives (predictive power of the newly constructed features and their ability to reconstruct the original input). The main contributions of this paper are threefold: (i) we investigate how GP can best be used for DR by considering standard DR loss functions and GP representations (see Section 3), (ii) additionally, we perform experiments with several DR loss functions that have previously never been applied to GP, such as weighted rank-correlation, using a neural-based autoencoder teacher, or directly using a GP-based autoencoder to reconstruct the input, and (iii), we show that our GP-based methods are able to learn non-linear and concise formulas of the lower-dimensional latent space. We benchmark our methods against three mainstream DR algorithms: principal components analysis [10] (PCA), locally linear embedding [25] (LLE), and isomap [29].

BACKGROUND 2.1 Dimensionality Reduction
The goal of DR is to find a lower-dimensional representation of the original data such that some properties of the original data, such as total variance (e.g., PCA), distances between points (e.g., isomap), ability to reconstruct the original features (e.g., autoencoders) are preserved. Formally, let us denote the original dataset by X ∈ R × where and represent the number of instances and features, respectively. Then, the goal is to find another representation X ∈ R × of X, such that ≪ , and such that the new representation retains some meaningful properties of the original data. Note that throughout this paper, we will refer to and as the th dimension (feature), and and˜as the th observation of the original and lower-dimensional representation, respectively.

Genetic Programming for Dimensionality Reduction
Recently, DR using GP has gained traction and has started being more thoroughly explored. In [17], the authors use GP to perform DR by using a fitness function that encourages the preservation of local neighborhoods. In [20], the authors build on their previous work [17] by taking a bi-objective approach where the first objective encourages to preserve local structure, while the second objective represents the number of dimensions in the latent representation. However, both in [17] and [20], there is no constraint imposed on the complexity of the evolved expressions, potentially rendering them non-interpretable. This issue is addressed in [18] where the authors use a bi-objective approach with the first objective being the cost function used in the t-SNE algorithm [21] and the second objective being a proxy for complexity. However, t-SNE is devised to strictly preserve local neighborhoods and does not result in a functional mapping, thus it is unclear whether the cost function of t-SNE is well-suited to be used for GP. In addition, GP has recently been used to perform feature extraction and transfer learning by constructing non-linear combinations of features to then be used in subsequent supervised learning tasks [16,31,38]. This type of application requires to specify a predictive task with ground truth labels. However, we may be interested in discovering intrinsic relationships between the variables or find an interpretable functional mapping for the compression of the original dataset without necessarily having an associated predictive task. One way to tackle the aforementioned problem was explored in [27], where the authors used GP to predict the value of each feature by using the remaining ones as input. However, this process does not belong to the realm of DR since the number of features remains the same.

METHODOLOGY
In GP, the main way to represent symbolic expressions is by defining it as a computational tree [24]. In this paper, we follow the literature and use trees to represent our GP programs. A detailed description of the tree-based representations we consider is given in the following subsection.

Computational Tree Representations
3.1.1 Single tree. The single-tree representation, as shown in Figure 1a, is one of the simplest way to represent a symbolic expression. In the single-tree representation, a single tree is responsible to output a unique single value. Formally, if we denote a single tree by , we have that : R → R, for = 1, . . . , . In words, each tree maps the original dimensions (the terminals at the leaves) to a single latent dimension (the root of the tree).
This representation is limited when it comes to modeling multioutput functions which is the case in DR, if > 1. Indeed, each latent dimension is evolved separately (see Figure 1a). Therefore, this representation cannot efficiently model non-separable fitness functions, where variables are dependent and should thus be jointly optimized, as is the case in DR.

Classic multi-tree.
Similarly to [17,18,20], we use a multitree representation as illustrated in Figure 1b. Each tree of the multi-tree takes the original data as input and outputs a single dimension of the lower-dimensional representation. During evolution, the trees of the multi-tree evolve simultaneously, different from single-tree representation, where each tree evolves separately during independent runs. This representation has two main advantages: (i) it only requires a single evolution run, (ii) it is able to learn non-separable fitness functions, where the variables are dependent, which is the case in cost functions that are typically used in DR.
In this paper, we will denote the multi-tree representation by represents the th tree of the multitree and is the number of dimensions in the lower-dimensional representation.
3.1.3 Autoencoder multi-tree. In addition to the classic multi-tree, we propose a new representation for DR (inspired from automatically defined functions [13]), composed of two classic multi-trees, with one multi-tree acting as encoder and one as decoder. Specifically, the encoder is a multi-tree composed of trees and is responsible to map the original input to its lower-dimensional representation. On the other hand, the decoder is a multi-tree composed of trees where is the number of original features. It is responsible to map the output of the encoder to reconstruct the original features. The encoder and the decoder are jointly evolved in a single run. Formally, let us denote the encoder and decoder of the autoencoder multi-tree by AMT ( ) and AMT ( ) respectively, then we have that AMT ( ) The purpose of the autoencoder multi-tree is to test whether a fully GP-based approach, i.e., one where both the encoder and decoder are evolved at the same time, is capable of performing DR effectively, as explained in the following subsection.

Fitness Functions
We investigate several fitness functions, denoted by F , to guide the evolutionary process of the multi-tree GP. As mentioned in Section 2.1, we denote the lower-dimensional representation of the original dataset X by X, where the th data point is given by = MT ( ), = 1, . . . , , for a given multi-tree individual MT . Furthermore, we denote the distance matrices of the original dataset and its Autoencoder multi-tree Figure 1: Three ways of mapping the original data to its lower-dimensional representation using computational trees. Note that the trees are read top-down, from leaves to output. (a) Classic: A GP individual is represented by a single tree. Each tree maps the input to a unique latent dimension, and is evolved separately. It thus requires independent runs. (b) Multi-tree: A GP individual is a collection of trees that are evolved jointly in a single evolution run. (c) Autoencoder multi-tree multitree: A GP individual is made of two parts: the encoder and the decoder. The encoder is a multi-tree composed of trees and is responsible to map the original input to its lower-dimensional representation. The decoder is a multi-tree composed of trees, where is the number of original features, and is responsible to map the output of the encoder to reconstruct the original features. The encoder and the decoder evolve jointly.
lower-dimensional representation by and , respectively, where , = , and , = , . For the distance preserving and rank preserving fitness functions, we investigate two distance metrics: (i) Euclidean distance, and, (ii) geodesic distance. The geodesic distance is particularly suited to represent the local low-dimensional geometry of the data manifold [29] that the Euclidean distance may not be able to capture.
Formally, for a fitness function denoted by F (X, ·), we seek MT * such that: (1) Note that for the GP-based autoencoder fitness function, we use an autoencoder multi-tree and the best individual would thus be denoted by AMT * .

Distance preserving.
This fitness function is based on preserving distances between pair of instances when mapping the original data to its lower-dimensional representation. To this end, we use Sammon mapping [26] which aims to minimize the following: The effect of , in the denominator is to put emphasis on points that are close to each-other in the original space. In other words, we give more importance in preserving local structure than global.

Rank preserving.
This fitness function is similar to preserving distances except we are only interested in the ranks of the distances. Essentially, this is a measure of the monotonicity of the relationship between the distances in the original and lower-dimensional representations. This has the advantage of being more flexible, less sensitive to outliers (i.e., large distances) and thus more easily learnable than the distance preserving fitness. In [20], the authors use Spearman's rank correlation coefficient [35] as the fitness function to perform dimensionality reduction using GP. In this paper, we use Kendall's which is preferred to Spearman's rank correlation coefficient for smaller samples [12]. Our aim is to maximize Kendall's : where is the number of concordant pairs, is the number of discordant pairs, and is the total number of pairs. Let us denote the th row of the distance matrices and by˜and respectively. Taking the average across instances, the final fitness function (to be minimized) is defined as: where , , and are computed using and . However, this approach fails to take into account that discordances between observations with high rank (a rank of 0 being the highest rank) are more important than those between items with low rank. This is because we want to emphasize preserving local structure over global structure, and so, if the distance between two points in the original space is large (i.e. the rank is low), we want it to have a lesser importance in the final fitness function. This serves similar purpose to the denominator in the Sammon mapping in Equation (2). Therefore, we use a weighted version of Kendall's proposed in [30], where the weighting is defined by ranking the original distances and assigning the highest rank of 0 to the shortest distance, and so on. Then, a weighting function assigns a weight depending on the rank . In this paper, we use a hyperbolic weighting function, and thus, we have that = 1 +1 .

Neural-based autoencoder teacher.
This fitness function is inspired by the student-teacher literature, where a (usually) simple model (e.g., linear regression, low-depth trees, symbolic regression) is trained to approximate the output of another (complex) model. This is a two-step process where we first train a base model and then train the surrogate model on the output (or hidden layer) produced by the base model. In our case, let us denote the latent layer of the neural-based autoencoder by L = { } =1 , where represents the th neuron and is the number of latent dimensions (e.g., the number of neurons in the latent layer). Our fitness function is then given by: where , = ( ) is the output of the th tree of the multi-tree MT = { } =1 . In words, the fitness function is the mean-squarederror between the output of the multi-tree and the latent layer of the autoencoder.

GP-based autoencoder.
This fitness function is similar to the neural-based autoencoder teacher fitness function in that it also makes use of an autoencoder. However, instead of evolving a classic multi-tree GP to match the latent layer of a neural autoencoder teacher, we let GP evolve a complete autoencoder, using the representation shown in Figure 1c. As mentioned, the GP-based autoencoder is composed of two multi-trees representing the encoder and the decoder of the autoencoder. In this case, the lower-dimensional representation of the original dataset X is given by the output of the encoding multi-tree. Formally, using the notation introduced in 3.1.3, we have that = AMT ( ) ( ) , = 1, . . . , , and the fitness function is given by: where , = , = ( ) ( ) is the output of the th tree of the GP's decoder AMT ( ) . In other words, the cost function is the reconstruction error between the output of the decoder and the original data, and the resulting lower-dimensional representation is given by the output of the encoder.

Baseline Methods
We compare the performance of the GP-based methods with three dimensionality reduction baselines: principal components analysis (PCA), locally linear embedding (LLE), and isomap. These methods are chosen because they learn a mapping on the training data that can be applied to unseen data, which is not possible with other methods such as t-SNE or multi-dimensional scaling (MDS). This enables a fair comparison since our GP-methods are trained to perform dimensionality reduction, and the evaluation is done on unseen data, as depicted in Figure 2. In addition, the baselines are heterogeneous as PCA seeks a global orthogonal linear transformation of the data, while LLE is a local and non-linear method, and isomap seeks to preserve the estimated intrinsic geometry of the data manifold using geodesic distance.

Performance Metrics
In this section, we discuss the choices made to evaluate the performance of the various fitness functions in guiding the evolution process of GP. This amounts to evaluating the quality of the resulting lower-dimensional representation of the original data. This is still an open problem in data science and will depend on the application and the questions that one wants to answer during the analysis. In addition, the fitness function being optimized during the dimensionality reduction process may be biased towards our evaluation metric or even be misaligned with it.
3.4.1 Predictive performance. In some situations, there exists a metric that truly represents our end goal, such as predictive performance on a supervised downstream task. If there exists such an objective (task) for a given dataset, then we do not have to worry whether the metric (e.g., classification accuracy or regression error) is biased towards any of our GP methods since it is a metric that we are objectively interested in. Using the predictive power of the low-dimensional representation on downstream supervised tasks to evaluate dimensionality reduction algorithms has been done in recent studies [17,18,20]. However, there are three main drawbacks with this approach. Firstly, we have to assume that the response (e.g., class structure for classification or target distribution for regression) is a significant factor within the manifold structure. Arguably, this assumption holds for the majority of curated datasets, and datasets for which supervised learning algorithms can perform well (i.e., one can learn labels from the data). Secondly, there may be some structure in the data that is not captured by the labels but that may still be of interest for knowledge discovery, and thus using predictive performance as a proxy would be insufficient. Finally, and perhaps the most critical drawback is that we may have datasets on which we want to perform dimensionality reduction but that have no labels.
In this work, we use balanced accuracy [3] as a performance metric. In particular, we follow [17,20] and use a random forest classifier as evaluation model for accuracy (denoted by in Section 4.4), trained on the low-dimensional representation to predict the labels and compute the predictive performance.

Reconstruction error.
Besides accuracy, we also consider the reconstruction error, as an additional metric, and, in this case, we use a neural-based decoder as the evaluation model. Note that, however, the reconstruction error is biased to favor the methods using the neural-based autoencoder teacher (F AE ) and the GP-based autoencoder (F GP ) fitness functions since they both minimize the reconstruction error during optimization.

EXPERIMENTAL SETUP 4.1 Datasets
We test our methods on 5 datasets from the UCI repository [6] for which a brief description is provided in Table 1. Note that for the purpose of our analysis, we only keep the continuous attributes, which is why the number of features is different from the original one for datasets that contain categorical variables. This is because we only use arithmetic operators and do not use any conditional or logical operators (e.g., if, and, or) in our GP representation (see Section 4.2). In addition, some of our fitness functions are based on the Euclidean distance, which would not be appropriate for a mix of continuous and categorical variables. Note that the fitness functions described in Section 3.2 make use of all the variables from the original data. This happens either at the time of computing the distance matrices or when computing the reconstruction error. In doing so, we may be using noisy variables as well as redundant variables (e.g., co-linear variables) which could be dominating other important variables during the optimization process. A simple and well established solution is to perform PCA before applying any sort of analysis that requires to compute distances between points, such as clustering analysis [2,5]. Thus, in this paper, we transform the original dataset using PCA (retaining 99% of the original variance) to compute the fitness functions in Section 3.2. This amounts to replacing X with X pca in F (X, ·), where X pca is the transformation of X after applying PCA and keeping a sufficient number of principal components to account for 99% of the original variance. Note that in order to keep the GP expressions interpretable in terms of the original input, the input to the multi-tree GP is still the original data X, as opposed to the PCA transformed data X pca .

Genetic Programming Parameters
We use standard GP parameter settings, as shown in Table 2. A basic tournament selection [24] process is used to chose the individuals that will act as parents to produce offspring via the variation operators. Each parent can be subject to up to three different genetic operators: (i) crossover, (ii) subtree mutation, and (iii) one-point mutation, with probabilities of , , and , respectively. These can be easily experimented with, using our publicly available code 1 . Note that the genetic operators, maximum tree size, and tree depth are applied to single trees within a multi-tree. This means that, for instance, each tree within a multi-tree has a probability of having a subtree mutation. For the crossover operator, where subtrees are exchanged between two trees, one has to decide whether crosspollination between trees at different index from the two multi-tree individuals is allowed. In [19], the authors randomly select a tree from each individual and perform standard crossover between them (random-index crossover), while in [1], the authors only allow trees at the same position in the multi-tree to crossover with each other (same-index crossover). Here, we use same-index crossover when mixing trees from different multi-trees. The reason we make this choice is that same-index crossover encourages trees at the same index to specialize [1] and converge towards learning the same latent dimension of the lower-dimensional representation.
Finally, the population is initialized using the popular ramped half-and-half method [14], and ephemeral constants are drawn from a N (0, 1). The terminal set (i.e., the set of all possible leaf nodes) is T = { , R}, where R represents the set of random constants and = { } =1 is the set of variables. The function set (i.e., the set of possible non-leaf nodes) is F = {−, +, ×}. This means that GP will evolve polynomials, with arbitrarily complex interactions only bounded by the depth of the tree. This results in simple expressions that can be used as a starting point to further interpretability studies on performing DR with GP. Furthermore, we conducted the same experiments using the extended function set F ′ = {−, +, ×, cos, log }, where log(·) = log(| · | + ), but obtained significantly worse results than when using F . Due to space limitations, the results obtained with F ′ are not presented here but can be found in the code repository.

Mini-batch Training
Note that during the evolutionary process, for each generation, we have to compute the distance matrix for all the multi-tree individuals in the population, which has an overall computational cost of O ( 2 × ), where is the population size used by GP, is the number of data points, and the number of latent variables of the lower-dimensional representation. For large values of , this can be prohibitively slow. Therefore, at each generation, we randomly sample a mini-batch of < data points from the data and compute the distance matrices on that mini-batch. In addition, mini-batch training is thought to improve generalization by changing the loss landscape at each batch, and thus escaping local optima [23].

Evaluation Protocol
A detailed account of the overall evaluation process, from the original data to the final predictions, is given in Figure 2. First, the original data is split into two subsets: (i) the training data for the DR model , and, (ii) a held-out set to apply on and obtain its lower-dimensional representation. In turn, this lower-dimensional data is split in two subsets: (i) the training data for the evaluation model , and, (ii) the hold-out set to apply the evaluation model to compute the performance metrics. The evaluation model (e.g., a neural decoder for reconstruction error and a random forest classifier for accuracy) is trained using a 10-fold cross validation and we report the results on the test set (see Section 5). Note that Figure 2: Evaluation protocol for our DR pipeline. First, we learn a DR model on a subset of the original data. This model is then applied on the held-out subset of the original data to obtain a lower-dimensional representation. We then learn an evaluation model on a subset of the lower-dimensional data and apply it to the held-out subset in order to compute the final performance metric.
the models used for the evaluation phase are described in Section 3.4.

Quantitative Analysis
In Table 3 (a), we display the average balanced accuracy and associated standard deviation computed across 30 independent runs for all the GP-based methods and the baselines. Similarly, in Table 3 (b), we display the results for the reconstruction error. In addition, we compute statistical significance using the Mann-Whitney U rank test [22], denoted by the asterisk notation.
We can see that the non-linear baseline methods LLE and isomap seem to yield the better performances, in particular when = 2. For = 3, F dist (Euclidean) performs best on 2 of the datasets but is not always statistically significant. Overall, it seems that for = 3, the difference between the methods is less pronounced than for = 2. This is intuitive since the more dimensions we keep the easier it is to capture relevant information from the original data, thus attenuating potential differences between methods. In fact, for more complex datasets with a higher number of features, such as Ionosphere and Libras Movement, we can still see a significant difference between the best performing method (isomap) and the other methods. Additionally, PCA (the only linear dimensionality reduction method in our experiments) does not perform well on these two datasets, suggesting that non-linearity allows for capturing relevant relationships between the original variables.
In Table 3 (b), we can see that F AE and F GP have the lowest reconstruction error out of the GP-based methods. This is expected because these methods minimize the reconstruction error during the optimization process. In addition, we can see that PCA and isomap yield good results across several datasets. Interestingly, while LLE performs well when balanced accuracy is used as metric, it is arguably the worst performing method to reconstruct the original data. This may be explained by the fact that co-linear features are equally weighted when computing the distance matrices and . This is detrimental to maximizing accuracy as there may be independent variables (e.g., non co-linear) that are important for the downstream classification task but which would be outweighted when computing the distance matrices. For this reason, we performed PCA before computing distances, as mentioned in Section 4.1.1. However, we only experimented with keeping 99% of the original variance and one would need to experiment further to identify the effect of performing PCA on the balanced accuracy and reconstruction error. Finally, while F AE and F GP do perform well, there exists drawbacks to these methods. For F AE we have to train a neural network autoencoder, which adds a layer of complexity to the overall method. That is, if our autoencoder teacher does not perform well in the first place, the GP-based student is bound to fail. On the other hand, while F GP does not require an autoencoder teacher model, it does not scale well with the number of input features since the GP-decoder requires a number of trees equal to the number of input features, as depicted in Figure 1c.
Overall, our results show that GP-based dimensionality reduction methods can be on a par with proven baseline methods. In addition, a major advantage of GP-based methods is that they have the potential to produce interpretable mappings in the form of symbolic expressions, which we discuss next.

Qualitative Analysis
In Table 4, we display the expressions that map the original data into a lower-dimensional representation, for the various GP-based methods. We can see that the resulting embeddings are small enough to have the potential of being interpretable. This is in contrast to the baseline methods where only PCA provide a global functional mapping from the original data to its lower-dimensional representation. However, PCA strictly returns a linear combination of all the variables, limiting its complexity and rendering it unintelligible Table 3: Score of (a) balanced accuracy and (b) reconstruction error (mean ± standard deviation of 30 runs). For a confidence level of 1 − , statistical significance is denoted as follows: * , * * , * * * , for = 0.1, 0.05, 0.01, respectively. Statistical testing is performed between the best method and the remaining 8 methods. If the mean performance is equal across several methods, the method with the lowest standard deviation is chosen as the best method. The best method and those that are not statistically significantly different are in bold. (

DISCUSSION
In this paper, we have seen that GP can be competitive with widely used DR methods, while producing small, non-linear, and potentially interpretable functional mappings of how the original data can be compressed into the latent dimensions. In particular, when the fitness used for GP is the same as the objective we want to optimize for, (e.g., reconstruction error), GP often outperforms baselines. A limitation of this paper is that we adopted classic GP evolutionary operators for variation and selection, yet state-of-the-art GP methods include more interesting mechanisms. For instance, one could use a GP version that is more expressive such as differentiable Cartesian GP [9] (DCGP). This is because in DCGP, one can optimize internal weights of the expression using gradient descent. Another example could be to adopt GP-GOMEA [32] as it is state-of-the-art to discover accurate yet small expressions [15].
Another limitation of this study is that we only used a fixed set of parameters and did not study how different GP representations and fitness functions inter-relate to different evolutionary budgets and settings. For example, we used a classic population size of 1000. However, the population size is a key parameter because it determines the supply of (high-order) building blocks available for the evolution. Some recent works in GP such as [28,32] suggest that better performance is only achievable when the population size is of tens of thousands or more. In light of this, future work on understanding the potential of GP for DR should include an analysis of the effect of important parameter settings. Moreover, future work should also include more baselines and loss functions in the comparison, including, e.g., the recenly introduced PaCMAP [34].
An interesting aspect emerging from our results is that DR methods achieving high accuracy often performed worse in terms of reconstruction error (and vice versa). This can happen because, when performing reconstruction, each reconstructed feature is equally important; whereas to answer a predictive task as per classification, only a subset of them may be.

CONCLUSION
In conclusion, we have considered different versions, i.e., representations and fitness functions of which two are novel (GP autoencoder representation and autoencoder teacher fitness), of genetic programming for dimensionality reduction. We have found that genetic programming is a competitive approach to obtain interpretable mappings for dimensionality reduction. We thus believe that genetic programming for dimensionality reduction is a promising research avenue.