Exploring Dimensionality Reductions with Forward and Backward Projections

Dimensionality reduction is a common method for analyzing and visualizing high-dimensional data across domains. Dimensionality-reduction algorithms involve complex optimizations and the reduced dimensions computed by these algorithms generally lack clear relation to the initial data dimensions. Therefore, interpreting and reasoning about dimensionality reductions can be difficult. In this work, we introduce two interaction techniques, \textit{forward projection} and \textit{backward projection}, for reasoning dynamically about scatter plots of dimensionally reduced data. We also contribute two related visualization techniques, \textit{prolines} and \textit{feasibility map} to facilitate and enrich the effective use of the proposed interactions, which we integrate in a new tool called \textit{Praxis}. To evaluate our techniques, we first analyze their time and accuracy performance across varying sample and dimension sizes. We then conduct a user study in which twelve data scientists use \textit{Praxis} so as to assess the usefulness of the techniques in performing exploratory data analysis tasks. Results suggest that our visual interactions are intuitive and effective for exploring dimensionality reductions and generating hypotheses about the underlying data.

INTRODUCTION
Dimensionality reduction (DR) is an effective technique for analyzing and visualizing high-dimensional datasets across domains, from sciences to engineering. Dimensionalityreduction algorithms such as principal component analysis (PCA) and multidimensional scaling (MDS) automatically reduce the number of dimensions in data while maximally preserving structures, typically quantified as similarities, correlations or distances between data points. This makes visualization of the data possible using conventional spatial techniques. For example, analysts generally use scatter plots to visualize the data after reducing the number of dimensions to two, encoding the reduced dimensions in a two-dimensional position.
DR Challenges: Most DR (also called manifold learning or distance embedding) algorithms are driven by complex numerical optimizations. Dimensions derived by these methods generally lack clear, easy-to-interpret mappings to the original data dimensions, forcing users to treat DR methods as black boxes. In particular, data analysts with limited experience in DR have difficulty in interpreting the meaning of the projection axes and the position of scatter plot nodes [39,7]. What do the axes mean? is probably users' most frequent question when looking at scatter plots in which points (nodes) correspond to dimensionally reduced data. Most scatter-plot visualizations of dimensionally reduced data are viewed as static images. One reason is that tools for computing and plotting these visualizations, such as Matlab and R, provide limited interactive exploration functionalities. Another reason is that few interaction and visualization techniques that go beyond brushing-and-linking or cluster-based coloring to allow dynamic reasoning with these visualizations.
Enriching User Experience with DRs: In this paper, we introduce two interactions, forward projection and backward projection (Figure 1), to help analysts explore and reason about scatter plot representations of dimensionally reduced data, facilitating a dynamic what-if analysis. We contribute two related visualization techniques, prolines and feasibility map, to facilitate the effective use of the proposed interactions. We also introduce Praxis, a new interactive DR exploration tool implementing our interaction and visualization techniques for data analysis.
Our techniques enable users to interactively explore: 1) the most important features that determine the vertical and horizontal axes of projections, 2) how changing feature values (dimensions) of a data point changes the point's projected location (two-dimensional representation) and 3) how changing the projected position of a data point changes the highdimensional values of that point. We demonstrate our techniques first using a PCA-based linear DR and then a nonlinear, deep autoencoder-based [21] DR.
We assess the computational effectiveness of our methods by analyzing their time and accuracy performances under varying sample and dimension sizes. We then conduct a user study in which twelve data scientists performed exploratory data analysis tasks using Praxis. The results suggest that our visual interactions are scalable and intuitive and can be effective for exploring dimensionality reductions and generating hypotheses about the underlying data.
We also observe that our techniques belong to a class of interactions that bidirectionally couple the data and its visual representation: dynamic visualization interactions [46]. We look at dynamic visualization interactions under the visual embedding model [13] and discuss the properties of effective interactions that the model suggests.

RELATED WORK
Our work is related to prior efforts in understanding and improving user experience with dimensionality reductions.

Direct Manipulation in DR
Direct manipulation has a long history in human-computer interaction [42,27,5] and visualization research (e.g. [40]). Direct manipulation techniques aim to improve user engagement by minimizing the perceived distance between the interaction source and the target object [22].
Developing direct manipulation interactions to guide DR formation and modify the underlying data is a focus of prior research [8,15,18,23,24,48]. For example, X/GGvis [8] supports changing the weights of dissimilarities input to the MDS stress function along with the coordinates of the embedded points to guide the projection process. Similarly, iPCA [23] enables users to interactively modify the weights of data dimensions in computing projections. Endert et al. [16] apply similar ideas to additional dimensionality-reduction methods while incorporating user feedback through spatial interactions in which users can express their intent by dragging points in the plane.
Our work is closely related to earlier approaches using direct manipulation to modify data in DR visualizations [23,45,38,10,31]. Like our forward projection and unconstrained backward projection techniques, iPCA enables interactive forward and backward projections for PCA-based DRs. However, iPCA recomputes full PCA for each forward and backward projection, and these can suffer from jitter and scalability issues, as noted in [23]. Using out-of-sample extrapolation, forward projection avoids re-running dimensionality-reduction algorithms. From the visualization point of view, this is not just a computational convenience, but also has perceptual and cognitive advantages, such as preserving the constancy of scatter-plot representations. For example, re-running (training) a dimensionality reduction algorithm with a new data sample added can significantly alter a two-dimensional scatter plot of the dimensionally reduced data, even though all the original inter-data point similarities may remain unchanged. In contrast to iPCA, we also enable users to interactively de-fine constraints on feature values and perform constrained backward projection.
We refer readers to a recent survey [37] for a detailed discussion of prior research on visual interaction with dimensionality reduction.

Visualization in DR Scatter Plots
Prior work introduces various visualizations in planar scatter plots of DRs, in order to improve the user experience by communicating projection errors [9,41,2,30], change in dimensionality projection positions [23], data properties and clustering results [41,12], and contributions of original data dimensions in reduced dimensions [17]. Low-dimensional projections are generally lossy representations of the original data relations: therefore, it is useful to convey both overall and per-point dimensionality-reduction errors to users when desired. Researchers visualized errors in DR scatter plots using Voronoi diagrams [2,30] and corrected (undistorted) the errors by adjusting the projection layout with respect to the examined point [9,41].
Biplot was introduced [17] to visualize the magnitude and sign of a data attribute's contribution to the first two or three principal components as line vectors in PCA. Prolines reduce to biplots when PCA is used for dimensionality reduction. Our proline construction algorithm is general and reflects the underlying out-of-sample extension method used. On the other hand, biplots are based on singular-value decomposition and always use PCA forward projection, regardless of the actual DR used. Additionally, prolines differ from biplots in being interactive visual objects beyond static vectors and are decorated to communicate distributional characteristics of the underlying data point.
Stahnke et al. [41] use a grayscale map to visualize how a single attribute value changes between data points in DR scatter plots. We use feasibility map, a grayscale map, to visualize the feasible regions in the constrained backward projection interaction.

Out-of-sample Extension and Back Projection for DR
We compute forward projections using out-of-sample extension (or extrapolation) [44]. Out-of-sample extension is the projection of a new data point into an existing DR (e.g. learned manifold model) using only the properties of the already computed DR. It is conceptually equivalent to testing a trained machine-learning model with data that was not part of the training set. For linear DR methods, out-of-sample extension is often performed by applying the learned linear transformation to the new data point. For autoencoders, the trained network defines the transformation from the high-dimensional to low-dimensional data representation [3].
Back or backward projection maps a low-dimensional data point back into the original high-dimensional data space. For linear DRs, back projection is typically done by applying the inverse of the learned linear DR mapping. For nonlinear DRs, earlier research proposed DR-specific backward-projection techniques. For example, iLAMP [14] introduces a backprojection method for LAMP [25] using local neighborhoods and demonstrates its viability over synthetic datasets [14]. Researchers also investigated general backward projection methods using radial basis functions [32,1], treating backward projection as an interpolation problem.
Autoencoders [21], neural-network-based DR models, are a promising approach to computing backward projections. An autoencoder model with multiple hidden layers can learn a nonlinear dimensionality reduction function (encoding) as well as the corresponding backward projection (decoding) as part of the DR process. Inverting DRs is, however, an illposed problem. In addition to augmenting what-if analysis, the ability to define constraints over a back projection can also ease the computational burden. Praxis also enables users to interactively set equality and boundary constraints over back projections through an intuitive interface.
We presented initial versions of forward projection, backward projection, and prolines earlier as part of Clustrophile, an exploratory visual clustering analysis tool [12]. We give here a focused discussion of our revised interaction and visualization techniques, introduce Praxis, a new visualization tool that implements our techniques for exploratory data analysis using DR, and provide a thorough computational and userperformance evaluation. The current work also introduces feasibility map, a new visualization technique to facilitate backward projection interactions.

TIONS
We demonstrate our methods first on principal component analysis (PCA), one of the most frequently used linear dimensionality-reduction techniques; note that the discussion here applies as well to other linear dimensionality-reduction methods. PCA computes (learns) a linear orthogonal transformation of the empirically centered data into a new coordinate frame in which the axes represent maximal variability. The orthogonal axes of the new coordinate frame are called principal components.
To reduce the number of dimensions to two, for example, we project the centered data matrix, rows of which correspond to data samples and columns to features (dimensions), onto the first two principal components, e 0 and e 1 . Details of PCA Figure 2: Through forward projection, a user can quickly explore how much the difference in the StudentSkills index value explains the planar projection difference between Portugal (blue node) and Korea. along with its many formulations and interpretations can be found in standard textbooks on machine learning or data mining (e.g., [4,20]).

Forward Projection
Forward projection enables users to interactively change the feature values of a data point x and observe how these hypothesized changes in data modify the current projected location y ( Figure 2). This is useful because understanding the importance and sensitivity of features (dimensions) is a key goal in exploratory data analysis. In the case of PCA, we obtain the two-dimensional position change vector ∆y by projecting the data change vector x onto the principal components:

Prolines: Visualizing Forward Projections
It is desirable to see in advance what forward projection paths look like for each feature. Users can then start inspecting the dimensions that look interesting or important.
Prolines visualize forward projection paths using a linear range of possible values for each feature and data point (Figures 3). Let x i be the value of the ith feature for the data point x. We first compute the mean µ i , standard deviation σ i , minimum min i and maximum max i values for the feature in the dataset and devise a range I = [min i , max i ]. We then iterate over the range with step size cσ i , compute the forward projections as discussed above, and then connect them as a path. The constant c controls the step size with which we iterate over the range and is set to c = σ i /8 for the examples shown here.
In addition to providing an advance snapshot of forward projections, prolines can be used to provide summary information conveying the relationship between the feature distribution and the projection space. We display along each proline a small light-blue circle indicating the position that the data point would assume if it had a feature value corresponding to the mean of its distribution; similarly, we display two small arrows indicating a variation of one standard deviation (σ i ) from the mean (µ i ). The segment identified by the range with thick green and red line segments, providing a basic directional statistical context. µ i is the mean of the ith dimension in the dataset, the green segment represents forward projections for is highlighted and further divided into two segments. The green segment shows the positions that the data point would assume by increasing its feature value; the red one indicates a decreasing value. This enables users to infer the relationship between the feature space and the direction of change in the projection space.

Backward Projection
Backward projection as an interaction technique is a natural complement of forward projection. Consider the following scenario: a user looks at a projection and, seeing a cluster of points and a single point projected far from this group, asks what changes in the feature values of the outlier point would bring the outlier near the cluster. Now, the user can play with different dimensions using forward projection to move the current projection of the outlier point near the cluster. It would be more natural, however, to move the point directly and observe the change.
The formulation of backward projection is the same as that of forward projection: ∆y = ∆x E. In this case, however, ∆x is unknown and we need to solve the equation.
As formulated, the problem is underdetermined and, in general, there can be an infinite number of data points (feature values) that project to the same planar position. Therefore, we propose both unconstrained and constrained backward projections, for which users can define equality as well as inequality constraints.
In the case of unconstrained backward projection, we find ∆x by solving a regularized least-squares optimization problem: Figure 4: Forward projection with prolines. StudentSkills is revealed as key feature differentiating Portugal from Korea. Observe that a value of 515 for StudenSkills would be reasonable with respect to the feature distribution (µ i < 515 < µ i + σ i ), but not enough to make Portugal close to Korea in the projection plane. By visually comparing the lengths (variability) of different proline paths, the user can easily recognize which dimensions contribute most to determining the position of points in the dimensionally reduced space. For instance, a change in the feature value of WorkingLongHours (the shortest proline) would produce only a very small change in the projection.
Note that this is equivalent to setting ∆x = ∆y E T .
For constrained backward projection, we find ∆x by solving the following quadratic optimization problem: minimize ∆x ∆x E − ∆y 2 subject to C∆x = d lb ≤ ∆x ≤ ub C is the design matrix of equality constraints, d is the constant vector of equalities, and lb and ub are the vectors of lower and upper boundary constraints.

Guiding Backward Projections
Projection Marks: It is important to note that, since more than one data point in the multidimensional space can project to the same position, forward and backward projections may not always correspond. For this reason, we add to our prolines visualization a set of projection marks ( Figure 6) dynamically indicating the current value for each feature while the user performs backward projection. At the same time, dragging a data point highlights the green or red segment of each proline based on the increase or decrease of each feature, showing which dimensions correlate to each other. By combining forward projection paths to backward projection, the user can infer how fast each value is changing in relation to its feature distribution.
Feasibility Maps: Constrained backward projection enables users to semantically regulate the mapping into unprojected high-dimensional data space. For example we don't expect an Age field to be negative or bigger than 150, even though such a value can can constitute a more optimal solution in an unconstrained backward projection scenario.
We propose the feasibility map visualization as a way to quickly see the feasible space determined by a given set of constraints. Instead of manually checking if a position in the projection plane satisfies the desired range of values (considering both equality and inequality constraints), it is desirable to know in advance which regions of the plane correspond to admissible solutions. In this sense, feasibility map is a conceptual generalization of prolines to the constrained backward projection interaction.
To generate a feasibility map, we sample the projection plane on a regular grid and evaluate the feasibility at each grid point based on the constraints imposed by the user, obtaining a binary mask over the projection plane. We render this binary Figure 6: With the addition of projection marks, prolines can be used as guides while the user performs backward projection. The figure shows the same constrained backward projection example as Figure 5: while the user drags the projected data point, prolines indicate with a green or red color if their corresponding feature value is increasing or decreasing. Projection marks, represented as little blue circles, move along each proline as the user performs backward projection and indicate the current feature values. This is particularly useful for checking the position of each value with respect to its feature distribution. In the case of inequality constraints, values that do not satisfy a constraint generate a black proline. mask over the projection as an interpolated grayscale heatmap, where darker areas indicate infeasible planar regions (Figure 7). With accuracy determined by the number of backward projection samples, the user can see which areas a data point can assume in the projection plane without breaking the constraints. Note that, when dealing with linear dimensionality reduction techniques, the feasibility map originates boundaries that are orthogonal to the prolines of the constrained features. Nevertheless, generating the feasibility map by sampling the projection plane has the advantage of being independent of the dimensionality technique used.
In backward projection, if a data point is dragged to a position that does not satisfy a constraint, its color and the color of its corresponding projection marks turn to black. If the user drops the data point in an infeasible position, the point is automatically moved through animation back to the last feasible position to which it was dragged (Figure 8).

PRAXIS
To study the usage of our interaction and visualization techniques, we introduce Praxis, a new interactive tool integrating them. Through a data panel (Figure 9b), users can load a dataset in CSV format and visualize its PCA projection as a scatter plot (Figure 9a), using the first two principal components as axes of the projection plane.
Results of forward and backward projection, along with the two visualizations prolines and feasibility map, are displayed in the projection plane. The id (name) of a data point is shown on mouse hover, while clicking performs selection, showing its feature values in a dedicated sidebar panel (Figure 9c). In particular, the Selection Details panel is used for perform-  ing forward projection (clicking on a dimension makes its value modifiable) and for inspecting changes in feature values when backward projection is used. Three buttons enable (respectively for each feature) 1) reset it to its original value, 2) enable/disable inequality constraints and 3) lock its value to a specific number (equality constraint).
Double-clicking the row associated with a feature displays a histogram representing its distribution below the selected row, showing some basic statistics (Figure 9d). The current value of the feature is represented by a blue line and a cyan line indicates the distribution mean. Bins of the histogram are colored similarly to prolines: green for increasing values, red for decreasing values with respect to the original feature value. Dragging one of the two black handles lets the user set or unset lower and upper bounds for a feature distribution, thus defining a set of constraints for a specific data point.
Finally, selecting a data point in the projection plane displays two buttons that respectively enable 1) resetting its feature values (and position) to their original value and 2) showing a tooltip on top of its k currently nearest neighbors, in order to facilitate reasoning about the similarity with other data samples (especially when performing backward projection).
Praxis is a web application based on a client-server model. We implemented its interface using Javascript with the help of D3 [6] and ReactJS [35] libraries. A separate analytics server carries out the computations required by the projection computations. We implemented the server in Python, using the SciPy [26], NumPy [47], scikit-learn [34] and CVX-OPT [11] libraries. We solve the backward projection generated quadratic optimization problems using CVXOPT [11].

EVALUATION
To evaluate our methods, we first conduct a user study with analysts and then perform a computational model analysis assessing accuracy and scalability.

User Study
We evaluate user experience with our techniques through a user study with twelve data scientists. We have two goals. First, to assess the effectiveness of our projection and visualization techniques in what-if analysis of dimensionally reduced data. Second, to understand how the use of the techniques differs for changing task type and complexity.
Participants: We recruited twelve participants with at least two years' experience in data science. Their areas of expertise included healthcare analytics, genomics and machine learning. Participants ranged from 25 to 55 years old, and all had at least a Master's degree in science or engineering. Ten participants regularly applied dimensionality reductions in their data analysis, using Matlab, R and Python. All participants were familiar with using PCA, while six had used MDS before. Four participants cited additional projection methods that they had previously used, including t-SNE and autoencoder.
Tasks and Data: Participants were asked to perform the following six high-level tasks using Praxis. We used a tabular dataset [41,33] containing eight socio-economic development indices for thirty-four countries belonging to the Organisation for Economic Co-operation and Development (OECD). The dataset was a CSV file with 34 rows and nine columns, where one column represented country names. Participants were free to use any combination of interactions and visualizations to complete a given task.  Procedure: The study took place in the experimenter's office; one user at a time used Praxis running on the experimenter's laptop. The study had three stages. In the first, participants were briefed about the experiment and filled out a pre-experiment questionnaire eliciting information about their experience in data science and use of dimensionality-reduction techniques. In the second stage, participants were introduced to the Praxis interface and to our techniques via a training dataset. Five minutes were dedicated to showing each user how to perform forward projection, backward projection and constraints formulation. Participants were then introduced to a new dataset and asked to complete the six tasks above. Task duration was manually timed and subject responses were collected through think-aloud protocol. Participants had at most two trials to perform each task; in the event of a failure, they moved on to the next task. For each task we also recorded whether the task was completed with or without the experimenter's help. In the last stage, participants were asked to performance techniques used Task C H I forward prolines backward feasibility  Table 1: Results of user study. Of a total of twelve participants, the table indicates for each task how many of them completed the task with no help (C), completed it with help (H) or did not complete the task (I). We also show the number of participants who used each of our proposed techniques to perform single tasks. Note that both prolines and feasibility map can give users visual information without requiring them to perform forward projection or backward projection (whereas, for instance, forward projection is intrinsically bound to prolines).
complete a post-experiment questionnaire to gather subjective feedback.
Results and Discussion: We adopt a task performance criterion similar to that employed in [41]. For each task, we count the number of participants who completed the task (C), completed the task with help (H) or were unable to complete the task (I). Similarly, we also report the frequency of the interaction (forward projection and backward projection) and visualization (prolines and feasibility map) techniques used by users to complete each task. We list these results in Table 1. Figure 10 shows the average log time spent on each task, inclusive of the cases in which the participant didn't complete the task. All the completed tasks were performed in less than 30 seconds. Task completion times and their standard errors reflect the intrinsic complexities of the tasks.
Overall, prolines proved to be a simple yet powerful visualization technique for exploring dimensionally reduced data as well as reasoning about the dimensionality reduction. Prolines were used 59 times by participants over the course of the six tasks performed. One participant mentioned he particularly liked "how prolines generate meaningful axes in a scatter plot where a clear mapping to data dimensions is unclear," whereas another described prolines as a "great way to understand dimensionality reductions, especially for people who used to treat them as a black box." The second most frequently used technique was backward projection, (used 19 times, followed by forward projection (15 times). Note that the use of forward projection always involves displaying prolines, which incorporates paths of forward projection locations for hypothetical feature values. Participants used forward projection when they wanted to interactively change the feature values and see precisely the projection change of the corresponding data point. In particular, one participant declared, "I feel backward projection is more natural to use and useful to see which features correlate to each other, but I would prefer forward projection for more precise control over feature values." Also, despite its lower incidence of use, feasibility map was employed by the participants when the task was sufficiently complex (T6).

Model Analysis
Since they are intended for use in interactive applications, the computational complexity of the proposed techniques needs to adhere to certain responsiveness requirements. At the same time, forward and backward projection methods need to be accurate enough at estimating changes in the dimensionally reduced space as well as in the underlying multidimensional data so as not to lead the user to false assumptions. We evaluate our proposed techniques for PCA in terms of time and accuracy over varying number of samples and dimensions of the input dataset and also over the amount of change introduced by the user (i.e. how much a feature value is modified in the case of forward projection, how much a data point is moved in the case of backward projection).
In our evaluation we iteratively perform forward projection and unconstrained backward projection on automatically gen-erated Gaussian random multivariate distributions, changing either the number of data samples or the number of data dimensions and leaving the other one fixed. We apply our techniques for each data point of the original distribution. The forward projection algorithm is then applied to each dimension with an amount of change in Time Performance: Figure 11 shows that the execution time for both forward projection (e,f) and backward projection (g,h) is on the order of microseconds and is not influenced by the number of samples nor by the amount of change; charts i and h show a linear dependence on the number of dimensions that does not, however, significantly affect the time performance. Even when dealing with larger datasets (e.g. > 500 samples, > 100 dimensions), both techniques are suitable for interactive data analysis tools. Figure 11 also shows the time performance of prolines (i,j) and feasibility map (k,l), respectively assuming the computation of each proline with an average resolution of 5 forward projection samples and the generation of the feasibility map with a resolution of 100 backward projection samples. In particular, we note that the time to compute prolines depends linearly on the number of dimensions.
Accuracy: To assess the accuracy of our techniques, we introduce a new similarity criterion for data-point neighborhoods in dimensionally reduced spaces. For each execution of the algorithms on a data point, we compute two sets of neighbors: 1) the n closest neighbors in the projection plane after performing forward projection or after moving the data point with backward projection, and 2) the n closest neighbors in the projection plane after performing the dimensionality reduction on the multidimensional data, after it has been modified through forward projection or by backward projection. Optimally, these two neighborhoods should contain the same elements, which should have the same relative distance from the data point on which the technique is performed. We define a neighborhood correlation index c n = c e × c o , where c e is the percentage of elements that appear in both neighborhoods, whereas c o is the percentage of elements whose distance from the data point considered remains in the same order. The index varies between 0 and 1, with 1 corresponding to very similar neighborhoods. Figure 11 shows that the accuracy of both forward projection (a,b) and backward projection (c,d) is mostly insensitive to the number of samples or dimensions. Instead, we notice a strong dependence on the amount of change introduced by the user. This shows that our proposed techniques are well suited for local changes in the data, and greater user modifications could possibly alter the properties of the dimensionality reduction.

DUCTIONS
We have so far demonstrated our interaction methods on PCA, a linear projection method. What about interacting with nonlinear dimensionality reductions? There are out-of-sample extrapolation methods for many nonlinear dimensionalityreduction techniques that make the extension of forward projection with prolines possible [3]. As for backward projection, its computation will be straightforward in certain cases (e.g. when an autoencoder [21] is used). In general, however, some form of constrained optimization specific to the dimensionality-reduction algorithm will be needed. Nonetheless, it is highly desirable to develop general methods that apply across dimensionality-reduction algorithms.
Below we discuss an application of our techniques to an autoencoder-based nonlinear dimensionality reduction.

Autoencoder-based dimensionality reduction
An autoencoder is an artificial neural network model that can learn a low-dimensional representation (or encoding) of data in an unsupervised fashion [36]. Autoencoders that use multiple hidden layers can discover nonlinear mappings between high-dimensional datasets and their low-dimensional representations. Unlike many other dimensionality-reduction methods, an autoencoder gives mappings in both directions between the data and low-dimensional spaces [21], making it a natural candidate for application of the interactions introduced here.
Figures 12 shows how backward projection can be applied in exploring a two-dimensional, autoencoder-based projection of handwritten digits from the MNIST database [29]. To this end, we first train an autoencoder with seven fully connected layers using the 60,000-digit images from the training set of the MNIST database.
Each image had 28 px × 28 px = 784 features. The seven layers of the autoencoder had sizes 784, 128, 32, 2, 32, 128 and 784. After training the network, we plot the encoded lowdimensional representations of 100 digits from the MNIST test set as circular nodes in the plane. We color each node by their associated digit. By performing backward projection on a data point, it is possible to observe the change in its feature values (pixel intensities) as the point is moved around the projection plane. Our technique shows, in this case, how one handwritten digit can gradually be transformed into another.
Using the same neural network model, we integrate support for autoencoder-based dimensionality reduction in our tool Praxis. In particular, Figure 13 shows prolines are not straight lines for non-linear dimensionality reduction methods. Figure 12: Backward projection with autoencoder. A user explores a dimensionality reduction of handwritten digits from the MNIST database [29] using backward projection. The twodimensional projection is obtained with a deep autoencoder. Projected data points are colored based on the digit represented.
By moving a node of the digit 7, back projection enables the user to see how its feature values (pixels) are updated. In particular, the user observes a smooth transition from 7 to 1, with 9 as an intermediate representation. f : X ! Y Figure 14: Visual embedding is a function that preserves structures in the data domain X within the embedded perceptual space Y (adapted from [13]). Figure 15: Bidirectionally coupling data and its visual representations. Visual embedding suggests that a change in the data should be reflected with a proportional change in its visual representation using f , the visual embedding function. Conversely, change in a visual representation should be reflected by a proportional change in the corresponding data using f −1 .

A MODEL FOR DYNAMIC VISUALIZATION INTERACTIONS
The interaction techniques introduced here belong to a class of interactions that tightly couples data with its visual representation so that when users interactively change one, they can observe a corresponding change in the other. For example, through forward projection, users observe how the visual representation (2D position) changes as they change the value of a dataset's attributes. Conversely, users can see how the data changes through backward projection as they change the visual representation. This class of interactions is essential for realizing dynamic visualizations (e.g. [46,28]) and we call them dynamic visualization interactions or dynamic interactions for short.
We now look at dynamic interactions under the visual embedding model [13]. The visual embedding model provides a functional view of data visualization and posits that a good visualization is a structure-or relation-preserving mapping from the data domain to the range (co-domain) of visual encoding variables ( Figure 14). Visual embedding immediately gives us criteria on which dynamic interactions should be considered effective ( Figure 15): 1) a change in data (e.g., induced by user through direct manipulation) should cause a proportional change in its visual representation and 2) a perceptual change in a visual encoding value (e.g., by dragging nodes in a scatter plot or changing the height of a bar in a bar chart) should be reflected by a change in data value that is proportional to the perceived change. However, to enable a dynamic interaction on a visualization, we need to have access to both the visual-ization function f and its inverse f −1 . The visual embedding model also suggests why implementing back mapping to the data space can be challenging.
We consider three basic forms of the visualization function f in Figure 16 through examples using a toy dataset in Figure 17. When the visualization function f is one-to-one, a dynamic interaction over f is straightforward as f −1 exists. When f is one-to-many (still invertible but not necessarily a proper function), f −1 exists and is determined by the target of interaction. Consider the example in Figure 16. We visualize how X values change for each NAME with a line chart. We also visualize the correlation of X and Y with a scatter plot. If a user moves a point up or down in the line chart, the corresponding change in X can be easily computed and the scatter plot can be updated in a brush-and-link fashion. Essentially, the one-to-many case can be seen as a collection of multiple one-to-one visualizations.
The most interesting case is when f is many-to-one and hence not invertible. A frequent source of such visualization functions is summary data aggregations, which are lossy. We can consider dimensionality reduction as a form of aggregation. A simple example of a many-to-one visualization is the bar chart ( Figure 16) that shows the mean X for the data points grouped by TYPE, A and B. Now, in a dynamic interaction scenario, how should we update the data values if a user changes the heights of the bars? Our backward projection solves a similar problem under a more complex visualization function, dimensionality reduction. In general, constructing a dynamic interaction over many-to-one visualization functions would require imposing a set of assumptions over data in the form of, e.g., constraints or models. This presents a challenging yet important future research direction. Figure 16: Visualization classes. Three basic forms of the visualization function f . Implementing a dynamic interaction is clearly challenging when f −1 does not exist.

VISUAL ANALYSIS IS LIKE DOING EXPERIMENTS
This paper introduces forward projection, backward projection, and the related visualizations, prolines and feasibility map, to improve user experience in exploratory data analysis using DR. We demonstrate these techniques on PCA-and autoencoder-based DRs. We also contribute a new tool, Praxis, implementing our techniques for DR-based data exploration.  Figure 17: Toy tabular dataset used for the three examples in Figure 16.
We evaluate our work through a computational performance analysis, along with a user study. Our visual interactions are scalable, intuitive to use and effective for DR-based exploratory data analysis. We situate our techniques in a class of visualization interactions at large that we discuss under the visual embedding model.
Data analysis is an iterative process in which analysts essentially run mental experiments on data, asking questions, (re)forming and testing hypotheses. Tukey and Wilk [43] were among the first in observing the similarities between data analysis and doing experiments, listing eleven similarities between the two. In particular, one of these similarities states that "Interaction, feedback, trial and error are all essential; convenience is dramatically helpful." In fact, data can be severely underutilized (e.g., dead [19,46]) without what-if analysis.
However, to perform data analysis as if we are running data experiments, dynamic visual interactions that bidirectionally couple data and its visual representations must be one of our tools. Our work here contributes to the nascent toolkit needed for performing visual analysis in a similar way to running experiments.


extfieldsofstudy:
Computer Science
Mathematics
Biology

s2fieldsofstudy:
Computer Science
Biology

text:
Algorithmically probable mutations reproduce aspects of evolution such as convergence rate, genetic memory, modularity, diversity explosions, and mass extinction

We show that if evolution is algorithmic in any form and can thus be considered a program in software space, the emergence of a natural algorithmic probability distribution has the potential to become an accelerating mechanism. We simulate the application of algorithmic mutations to binary matrices based on numerical approximations to algorithmic probability, comparing the evolutionary speed to the alternative hypothesis of uniformly distributed mutations for a series of matrices of varying complexity. When the algorithmic mutation produces unfit organisms---because mutations may lead to, for example, syntactically useless evolutionary programs---massive extinctions may occur. We show that modularity provides an evolutionary advantage also evolving a genetic memory. We demonstrate that such regular structures are preserved and carried on when they first occur and can also lead to an accelerated production of diversity and extinction, possibly explaining natural phenomena such as periods of accelerated growth of the number of species (e.g. the Cambrian explosion) and the occurrence of massive extinctions (e.g. the End Triassic) whose causes are a matter of considerable debate. The approach introduced here appears to be a better approximation to actual biological evolution than models based upon the application of mutation from uniform probability distributions, and because evolution by algorithmic probability converges faster to regular structures (both artificial and natural, as tested on a small biological network), it also approaches a formal version of open-ended evolution based on previous results. The results validate the motivations and results of Chaitin's Metabiology programme. We also show that the procedure has the potential to significantly accelerate solving optimization problems in the context of artificial evolutionary algorithms.

ALGORITHMIC EVOLUTION
C Entral to modern synthesis and general evolutionary theory is the understanding that evolution is gradual and is explained by small genetic changes in populations over time [1]. Genetic variation in populations can arise by chance through mutation, with these small changes leading to major evolutionary changes over time.
Of interest in connection to the possible links between the theory of biological evolution and the theory of infor-• S. Hernández-Orozco is a PhD student at Posgrado en Ciencia e Ingeniería de la Computación, Universidad Nacional Autónoma de México (UNAM), Mexico. • To whom correspondence should be addressed: hector.zenil@algorithmicnaturelab.org Contributions: HZ and SHO conceived the project. HZ and NAK provided guidance, data and proposed experiments. SHO, HZ and NAK analyzed the data. SHO and HZ wrote code. SHO and HZ wrote the paper. The results can be reproduced using the Online Algorithmic Complexity Calculator at http://www.complexitycalculator.com/ mation is the place and role of randomness in the process that provides the variety necessary to allow organisms to change and adapt over time.
It has been suggested [2], [3], [4], [5] that the deeply informational and computational nature of biological organisms makes them amenable to being studied or considered as computer programs following (algorithmic) random walks in software space, that is, the space of all possibleand valid-computer programs. Here we numerically test this hypothesis and explore the possible consequences of its validation vis-a-vis our understanding of the biological aspects of life and natural evolution by natural selection, as well as for applications to optimization problems in areas such as evolutionary programming.

Classical v. Algorithmic Probability
The classical probability of producing the first n digits of a mathematical constant like π in binary by chance, e.g., by typing on a typewriter, is exponentially unlikely as a function of the number of digits. However, because π is not random, in the sense that it has a short description that can generate an arbitrary number of digits of π with the same short formula, the algorithmic probability of π is much higher than its classical probability. This is because the classical probability of producing a short computer program encoding a short mathematical formula by chance is more likely than typing the digits of π themselves one by one. A π-generating formula can be written as a computer program in no more than N bits, Thus it has the classical probability of occurring by chance for the first n digits of π, where log(n) only accounts for the size of the number of digits n to be introduced into the computer program before it halts, compared to 1/2 n for the classical probability of producing n digits of π directly. Clearly, 1/2 N + log(n) and 1/2 n are divergent and for n > N + log(n) on 1/2 N + log(n) < 1/2 n .

Chaitin's Evolutionary Model
In the context of his work on Metabiology, Chaitin developed a theoretical computational model that evolves relative to its environment significantly faster than classical random mutation [2], [6], [7]. Starting with an empty binary string, the system approximates Chaitin's constant Ω, defined as Ω = p∈HP 2 −|p| where HP is the set of all halting programs [8], in an expected time of O(t 2 (log t) (1+O(1) )), which is significantly faster than the exponential time that the process is expected to take when randomly choosing mutations from a uniform distribution. This speed-up is obtained by drawing mutations according to the Universal Distribution [9], [10], which is explained in the following chapter.
In a previous result [11], it was shown that the previous system exhibits open-ended evolution (OEE, [12]) according to a strict definition of the concept, and that no decidable system with computable dynamics can achieve OEE. Here we will present a system that, by following the Universal Distribution, approaches this behaviour.

Algorithmic Probability and the Universal Distribution
At the core of our approach is the concept of Algorithmic Probability introduced by Solomonoff [13], Levin [14] and Chaitin [15]. Denoted by P (s), the algorithmic probability of a binary string s is formally defined as: where p is a random computer program in binary (whose bits were chosen at random) running on a so-called prefixfree (in order to constrain the number of valid programs as it would happen in physical systems) universal Turing machine U that outputs s and halts. Algorithmic probability connects the algorithmic likelihood of s to the intrinsic algorithmic s. The less algorithmically complex s (like π), the more frequently it will be produced on U by running a random computer program p. If K(s) is the descriptive algorithmic complexity of s (also known as Kolmogorov-Chaitin complexity [16], [17]), we have it that P (s) = 1 2 K(s)+O (1) . The distribution induced by P (s) over all strings is called the Universal Distribution or Levin's semi-measure, because the measure is semi-computable and can only be approximated from below and its sum does not add up to 1 to be a full measure.

Theoretical Considerations
In Chaitin's evolutionary model [2], [6], [7], a successful mutation is defined as a computable function µ, chosen according to the probabilities stated by the Universal Distribution [9], [10], that changes the current state of the system (as an input of the function) to a better approximation of the constant Ω [17]. In order to be able to simulate this system we would need to compute the Universal Distribution and the fitness function. However, both the Universal Distribution and the fitness function of the system require the solution of the Halting Problem [8], which is incomputable. Nevertheless, as with Ω itself, this solution can be approximated [18], [19]. In this project we are proposing a model that, to the best of our knowledge, is the first computable approximation to Chaitin's proposal.
For this first approximation we have made four important initial concessions: one with respect to the real computing time of the system, and three with respect to Chaitin's model: • We assume that building the probability distributions for each instance of the evolution takes no computational time, while in the real computation this is the single most resource-intensive step.

•
The goal of our system is to approximate objects of bounded information content: binary matrices of a set size. • We use BDM and Shannon's entropy as approximations for the algorithmic information complexity K. • We are not approximating the algorithmic probability of the mutation functions, but that of their outputs.
We justify the first concession in a similar fashion as Chaitin: if we assume that the interactions and mechanics of the natural world are computable, then the probability of a decidable event occurring is given by the Universal Distribution. The third one is a necessity, as the algorithmic probability of an object is incomputable (it requires a solution for HP too). In the next section we will show that Shannon's entropy is not as good as BDM for our purposes. Finally, note that given the Universal Distribution and a fixed input, the probability of a mutation is in inverse proportion to the descriptive complexity of its output, up to a constant error. In other words, it is highly probable that a mutation may reduce the information content of the input but improbable that it may increase the information content. Therefore, the last concession yields an adequate approximation, since a low information mutation can reduce the descriptive complexity of the input but not increase it in a meaningful way.

Our Expectations
It is important to note that, when compared to Chaitin's metabiology model [2], we changed the goal of our system therefore we must also the expectations we had for its behaviour.
Chaitin's evolution model [2] is faster than regular random models despite targeting a highly random object, thanks to the fact that positive mutations have low algorithmic information complexity and hence a (relatively) high probability of being stochastically chosen under the Universal Distribution. The universally low algorithmic complexity of these positive mutations relies on the fact that, when assuming an oracle for HP, we are also implying a constant algorithmic complexity for its evaluation function and target, since we can write a program that verifies if a change on a given approximation of Ω is a positive one without needing a codification of Ω itself.
In contrast, we expected our model to be sensitive with respect to the algorithmic complexity of the target matrix, obtaining high speed-up for structured target matrices that decreases as the algorithmic complexity of the target grows. However, this change of behaviour remains congruent with the main argument of metabiology [2] and our assertion that, contrary to regular random mutations, algorithmic probability driven evolution tends to produce structured novelty at a faster rate, which we hope to prove in the upcoming set of experiments.
In summary, we expect that when using an approximation to the Universal Distribution: • Convergence will be reached in significantly fewer total mutations than when using the uniform distribution for structured target matrices.

•
The stated difference will decrease in relation to the algorithmic complexity of the target matrix.
We also aimed to explore the effect of the number of allowed shifts (mutations) on the expected behaviour.

The Unsuitability of Shannon's Entropy
As shown in [19], [20], when compared to BDM we can think of Shannon's entropy alone as a less accurate approximation to the algorithmic complexity of an object (if its underlying probability distribution is not updated by a method equivalent to BDM, as it would not be by the typical uninformed observer). Therefore we expect the entropyinduced speed-up to be consistently outperformed by BDM when the target matrix inclines away from algorithmic randomness and has thus some structure. Furthermore, as random matrices are expected to have a balanced number of 0's and 1's, we anticipated the performance of single bit entropy to be nearly identical to the uniform distribution on unstructured (random) matrices. For block entropy [21], [22], that is, the entropy computed over submatrices rather than single bits, the probability of having repeated blocks is in inverse proportion to their size, while blocks of smaller sizes approximate single bit entropy, again yielding similar results to the uniform distribution. The results support our assumptions and claims.

Evolutionary Model
Broadly speaking, our evolutionary model is a tuple S, S, M 0 , f, t, α , where: • S is the state space, • M 0 , with M 0 ∈ S, is the initial state of the system, A successful evolution is the sequence M 0 , (M 1 , t 1 ) = S(M 0 , f, t), ..., ( , t n )) and t i is the total evolution time. We say that the evolution failed, or that we got an extinction, if instead we finish the process by (⊥, t n ), with t i being the extinction time. The evolution is undetermined otherwise. Finally, we will call each element (M i , t i ) an instance of the evolution.

Experimental Setup: A Max One Problem Instance
For this experiment, our phase state is the set of all binary matrices of sizes n × n, our fitness function is defined as the Hamming distance f (M ) = H(M t , M ), where M t is the target matrix, and our convergence parameter is α = 0. In other words, the evolution converges when we produce the target matrix, guided only by the Hamming distance to it, which is defined as the number of different bits between the input matrix and the target matrix.
The stated setup was chosen since it allows us to easily define and control the descriptive complexity of the fitness function by controlling the target matrix and, therefore also control the complexity of the evolutionary system itself. Is important to note that our setup can be seen as a generalization of the Max One problem [23], where the initial state is a binary "target gene" and the target matrix is the "target gene"; when we obtain a Hamming distance of 0 we have obtained the gene equality.

Evolution Dynamics
The main goal of this project is to contrast the speed of the evolution when choosing between two approaches to determining the probability of mutations: • When the probability of a given set of mutations has a uniform distribution. That is, all possible mutations have the same probability of occurrence, even if under certain constraints.
• When the probability of a given mutation occurring is given by an approximation to the Universal Distribution (UD) [10], [14]. As the UD is non-computable, we will approximate it by approximating the algorithmic complexity K ( [16], [17]) by means of the Block Decomposition Method (BDM, with no overlapping) [19] based on the Coding Theorem Method (CTM) [24], [25], [26] (see appendix A). • We will also investigate the results by running the same experiments using Shannon Entropy instead of BDM to approximate K.
Each evolution instance was computed by iterating over the same dynamic. We start by defining the set of possible mutations as those that are within a fixed n number of bits from the input matrix. In other words, for a given input matrix M , the set of possible mutations in a single instance is defined as the set Then, for each matrix in M, we compute the probability P (M ) is defined as: |M| in the case of the Uniform Distribution. For implementation purposes, we used a minor variation to the entropy probability distribution to be used and compared to BDM. The probability distributions for the set of possible mutations using entropy were built using two heuristics: Let M be a possible mutation of M , then the probability of obtaining M as a mutation is defined as either, The first definition assigns a linearly higher probability to mutations with lower entropy. The second definition is consistent with our use of BDM in the rest of the experiments. The constant is an arbitrary small value that was included to avoid undefined (infinite) probabilities. For the presented experiments was set at 1 −10 .
Once the probability distribution is computed, we set the number of steps as 0 and then, using a (pseudo)random number generator (RNG), we proceed to stochastically draw a matrix from the sated probability distributions and evaluate its fitness with the function f , adding 1 to the number of steps. If the resultant matrix does not show an improvement in fitness, we draw another matrix and add another 1 to the number of steps, not stopping the process until we obtain a matrix with better fitness or reach the extinction threshold. We can either replace the drawn matrix or leave it out of the pool for the next iterations. A visualisation of the stated work flow for a 2 × 2 matrix is shown in Figure 1.
To get a complete evolution sequence, we iterate the stated process until either convergence or extinction is reached. As stated before, we can choose to not replace an evaluated matrix from the set of possible mutations in each instance, but we chose to not keep track of evaluated matrices after an instance was complete. This was done in order to keep open the possibility of dynamic fitness functions in future experiments.
In this case, the evolution time is defined as the sum of the number of steps (or draws) it took the initial matrix to reach equality with the target matrix. When computing the evolution dynamics by one of the different probability distribution schemes we will denote it by uniform strategy, BDM strategy or h strategy, respectively. That is, the uniform distribution, the distribution for the algorithmic probability estimation by BDM, and the distribution by Shannon entropy. Fig. 1. An evolution instance. The instances are repeated by updating the input matrix until convergence or extinction is reached.

The Speed-Up Quotient
We will measure how fast (or slow) a strategy is compared to the uniform by the speed-up quotient, which we will define as: Definition 1. The speed-up quotient, or simply speed-up, between the uniform strategy and a given strategy f is defined as where S u is the average number of steps it takes a sample (a set of initial state matrices) to reach convergence under the uniform strategy and S f is the average number of steps it takes under the f strategy.

Cases of Negative Speed-up
In order to better explain the choices we have made to our experimental setup, first we will present a series of cases where we obtained no speed-up or slow-down. Although these cases were expected, they shed important light on the behaviour of the system.

Entropy vs. Uniform on Random Matrices
For the following experiments, we generated 200 random matrices separated into two sets: initial matrices and target matrices. After pairing them based on their generation order we evolved them using 10 strategies: the uniform distribution, block Shannon's entropy for blocks of size 4×4, denoted below by h b , entropy for single bits denoted by h, and their variants where we divide by h 2 and h 2 b respectively. The strategies were repeated for 1-and 2-bit shifts (mutations).
The results obtained are summarized in the table 1, which lays out the strategy used for each experiment, the number of shifts/mutations allowed, the average number of steps it took to reach convergence, as well as the standard error of the sample mean. As we can see, the differences in the number of steps required to reach convergence are not significant, validating our assertion that, for random matrices, entropy evolution is not much different than the uniform evolution.
Because the algorithmic complexity of a network makes sense only in its unlabelled version in general, and in most of the cases. In [19], [26], [27] we showed, both theoretically and numerically, that approximations of algorithmic complexity of adjacency matrices of labelled graphs are a good approximation (up to a logarithmic term or the numerical precision of the algorithm) of the algorithmic complexity of the unlabelled graphs. This means that we can consider any adjacency matrix of a network a good representation of the network disregarding graph isomorphisms.

Entropy vs. Uniform on a Highly Structured Matrix
For this set of experiments, we took the same set of 100 8 × 8 initial matrices and evolved them into a highly structured matrix, which is the adjacency matrix of the star with 8 nodes. For this matrix, we expected entropy to be unable to capture its structure, and the results obtained accorded with our expectations. The results are shown in table 2.
As we can see from the results, entropy was unable to show a statistically significant speed-up compared to the uniform distribution. Over the next sections we show that we have obtained a statistically significant speed-up by using the BDM approximation to algorithmic probability distributions, which is expected because BDM manages to better capture the algorithmic structures of a matrix rather than just the distribution of the bits which entropy measures. Based on the previous experiments, we conclude that entropy is thus not a good approximation for K, and we will omit its use in the rest of the article.

Randomly Generated Graphs
For this set of experiments, we generated 200 random 8 × 8 matrices and 600 16 × 16 matrices, both sets separated into initial and target matrices. We then proceeded to evolve the initial matrix into the corresponding target by the following strategies: uniform and BDM within 2-bit and 3-bit shifts (mutations) for the 8 × 8 matrices and only 2-bit shifts for the 16 × 16 matrices due to computing time. The results obtained are shown in the Figure 2. In all cases, we do not replace drawn matrices and the extinction threshold was set at 2500.  From the results we can see two important behaviours for the 8 × 8 matrices. The matrices generated are of high BDM complexity and evolving the system using the uniform strategy tends to be faster than using BDM for these highly random matrices. Secondly, although increasing the number of possible shifts by 1 seems, at a first glance, a small change in our setup, it has a big impact on our results: the number of extinctions has gone from 0 for both methods to 92 for the uniform strategy and 100 for BDM. This means that most evolutions will rise above our threshold of 2500 drafts for a single successful evolutionary step, leading to an extinction. As for the 16 × 16 matrices, we can see a formation of two easily separable clusters that coincide perfectly with the Uniform and BDM distributions respectively.

The Causes of Extinction
For the uniform distribution, the reason is simple: the number of 3-bit shifts on 8 × 8 matrices gives a space of possible mutations of 8×8 3 = 41664 matrices, which is much larger than the number of possible mutations present within 2shifts and 1-shift (mutation), which are 8×8 2 = 2016 and 8 × 8 = 64 respectively. Therefore, as we get close to convergence, the probability of getting the right evolution, if the needed number of shifts is two or one, is about 0.04%, and removing repeated matrices does not help in a significant way to avoid extinction,since 41 664 is much larger than 2500.
Given the values discussed, we have chosen to set the extinction threshold at 2500 and the number of shifts at 2 for 8 × 8 matrices, as allowing just 64 possible mutations for each stage is a number too small for showing a significant difference in the evolutionary time between the uniform and BDM strategies, while requiring evolutionary steps of 41 664 for an evolutionary stage is too computationally costly. The threshold of 2500 is close to the number of possible mutations and has been shown to consume a high amount of computational resources. For 16 × 16 matrices, we performed 1-bit shifts only, and occasionally 2-bit shifts when computationally possible.

The BDM Strategy, Extinctions and Persistent Structures
The interesting case is the BDM strategy. As we can see clearly in Figure 3 for the 8×8 3-bit case, the overall number of steps needed to reach each extinction is often significantly higher than 2500 under the BDM strategy. This behaviour cannot be explained by the analysis done for the uniform distribution, which predicts the sharp drop observed in the blue curve. After analyzing the set of matrices drawn during failed mutations (all the matrices drawn during a single failed evolutionary stage), we found that most of these matrices have in common highly regular structures. We will call these structures persistent structures. Formally, regular structures can be defined as follows: The previous inequality is a consequence of the fact that the possible mutations are finite and only a small number of them, if any, can have a smaller algorithmic complexity than the mutations that contain Γ; otherwise we contradict the existence of C. In other words, as Γ has relatively low complexity, the structures that contain Γ tend to also have low algorithmic complexity, and hence a higher probability of being chosen.
Finally, as shown in the section 3.2, we can expect the number of mutations with persistent structures to increase in factorial order with the number of possible mutations and in polynomial order with respect to the size of the matrices that compose the state space.

Proposition 3.
As a direct consequence of the last statement, we have it that, for systems evolving as described in the section 2.6 under the Universal Distribution: • Once a structure with low descriptive complexity is developed, it is exponentially hard to get rid of it.

•
The probability of finding a mutation without the structure decreases in factorial order with respect to the set of possible mutations.
• Evolving from and to unrelated regular structures is also hard.
Given the fourth point, we will always choose random initial matrices from now on, as the probability of drawing a mutation other than an empty matrix (of zeroes), when one is present in the set of possible mutations, is extremely low (below 9 × 10 −6 for 8 × 8 matrices with 2 shifts).

Positive Speed-Up Instances
In the previous section, we established that the BDM strategy yields a negative speed-up when targeting randomly generated matrices, which are expected to be of high algorithmic information content or unstructured. However, as stated in section 2.3, that behaviour is within our expectations. In the next section we will show instances of positive speed-up, including cases where previously entropy failed to show significant speed-up or was outperformed by BDM.

Manually Built Structured Matrices
For the following set of experiments we manually built three 8 × 8 matrices that encode the adjacency matrices of three undirected non-random graphs with 8 nodes that are intuitively structured: the complete graph, the star graph and a grid. The matrices used are shown in Figure 4. After evolving the same set of 100 randomly generated matrices for the three stated matrices, we can report that we found varying degrees of positive speed-up, that correspond to their respective descriptive complexities as approximated by their BDM values. The complete graph, along with the empty graph, is the graph that has the lowest approximated descriptive complexity with a BDM value of just 24.01. As expected, we get the best speed-up quotient in this case. After the complete graph, the star intuitively seems to be one of the less complex graphs we can draw. However, its BDM value (105.434) is significantly higher than the grid (83.503). Accordingly, the speed-up obtained is lower. The results are shown in the Figure 5.  As we can see from the Figure 5, a positive speed-up quotient was consistently found within 2-bit shifts without replacements. We have one instance of negative speed-up with one shift with replacements for the grid, and negative speed-up for all but the complete graph with two shifts.
However, it is important to say that almost all the instances of negative speed-up are not statistically significant, as we have a very high extinction rate of over 90%, and the difference between the averages is lower than two standard errors of the mean. The one exception is the grid at 1-bit shift, which had 45 extinctions for the BDM strategy. The complete tables are presented in the appendix C.

Mutation Memory
The cause of the extinctions found in the grid are what we will call maladaptive persistent structures (definition 2), as they occur at a significantly higher rate under the BDM distribution. Also, as the results suggest, a strategy to avoid this problem is adding memory to the evolution. In our case, we will not replace matrices already drawn from the set of possible mutations.
We do not believe this change to be contradictory to the stated goals, since another way to see this behaviour is that the Universal Distribution dooms (with very high probability) populations with certain mutations to extinction, and evolution must find strategies to eliminate these mutations fast from the population. This argument also implies that extinction is faster under the Universal Distribution than regular random evolution when a persistent maladaptive mutation is present, which can be seen as a form of mutation memory. This requirement has the potential to explain evolutionary phenomena such as the Cambrian explosion, as well as mass extinctions: once a positively structured mutation is developed, further algorithmic mutations will keep it (with a high probability), and the same applies to negatively structured mutations. This can also explain the recurring structures found in the natural world. Degradation of a structure is still possible, but will be relatively slow. In other words, evolution will remember positive and negative mutations (up to a point) when they are structured.
From now on, we will assume that our system has memory and that mutations are not replaced when drawn from the distribution.

The Speed-Up Distribution
Having explored various cases, and found several conditions where negative and positive speed-up are present, the aim of the following experiment was to offer a broader view of the distribution of speed-up instances as functions of their algorithmic complexity.
For the 8 × 8 case, we generated 28 matrices by starting with the undirected complete graph with 8 nodes, represented by its adjacency matrix, and then we removed one edge at a time until the empty graph (the diagonal matrix) was left, obtaining our target matrix set. It is important to note that the resultant matrices are always symmetrical. The process was repeated for the 16 × 16 matrices, obtaining a total of 120 16 × 16 target matrices.
For each target matrix in the first target matrix set, we generated 50 random initial matrices and evolved the population until convergence was reached using the two stated strategies: uniform and BDM, both without replacements. We saved the number of steps it took for each of the 2800 evolutions to reach convergence and computed the average speed-up quotient for each target matrix. The stated process was repeated for the second target matrix set, but by generating 20 random matrices for each of the 120 target matrices to conserve computational resources. The experiment was repeated for shifts of 1 and 2 bits and the extinction thresholds used were 2500 for 8 × 8 and 10 000 for 16 × 16 matrices.
As we can see from the results in Figure 6, the average number of steps required to reach convergence is significantly lower when using the BDM distribution for matrices with low algorithmic complexity, but the difference drops along with the complexity of the matrices but never crosses the Sequence of Generated Target Matrices extinction threshold. This suggests that symmetry over the diagonal is enough to guarantee a degree of structure that can be captured by BDM. It is important to report that we found no extinction case for the 8 × 8 matrices, 13 in the 16 × 16 matrices with 1-bit shifts, all for the BDM distribution, and 1794 with 2-bit shifts, mostly for the uniform distribution.
This last experiment was computationally very expensive. Computing the data required for the 16 × 16, 2-bit shifts sequence took 12 days, 6 hours and 22 minutes on a single core of an i5-4570 PC with 8GB of RAM. Repeating this experiment for 3-bit shifts is infeasible with our current setup, as it would take us roughly two months shy of 3 years. Now, by combining the data obtained for the previous sequence and the random matrices used in section 3.1.3, we can approximate the positive speed-up distribution. Given the nature of the data, this approximation (Figure 7) is given as two curves, each representing the expected evolution time from a random initial matrix as a function of the algorithmic information complexity of the target matrix for both strategies, uniform and BDM respectively. The positive speed-up instances are those where the the BDM curve is below the uniform curve.
The first result we get from Figure 7 is a confirmation of an expected one: unlike the uniform strategy, the BDM strategy is highly sensitive to the algorithmic information content of the target matrix. In other words, it makes no difference for a uniform probability mutation space whether the solution is structured or not, while an algorithmic probability driven mutation will naturally converge faster to structured solutions.
The results obtained expand upon the theoretical development presented in section 3. mutations grows, so do the instances of persistent structures and the slow-down itself. This behaviour is evident given that, when we increase the dimension of the matrices, we obtain a wider gap within the intersection point of the two curves and the expected BDM value, which corresponds to the expected algorithmic complexity of randomly generated matrices. However, we also increase the number of structured matrices, ultimately producing a richer and more interesting evolution space.

Chasing A Biological Network
We now set as target the adjacency matrix of a biological network corresponding to the topology of an ERBB signalling network [28]. The network is involved in responses ranging from cell division, death, motility, and adhesion and when dysregulated it has been found to be strongly related to cancer [29], [30]. As one of our main hypotheses is that algorithmic probability is a better model for explaining biological diversity, it is important to explore whethernaturally occurring structures are more likely to be produced under the BDM strategy than the uniform strategy, which is equivalent to showing them evolving faster.
The binary target matrix is shown in Figure 8 and it has a BDM of 349.91 bits. For the first experiment, we generated 50 random matrices that were evolved using 1bit shift mutations for the Uniform and BDM distributions, without repetitions. The BDM of the matrix is at the right of the intersection point inferred by the cubic models shown in Figure 7. Therefore we predict a slow-down. The results obtained are shown in the table 3.
As the results show, we obtained a slow-down of 0.71, without extinctions. However, as mentioned above, the BDM of the target matrix is relatively high, so this result  is consistent with our previous experiments. However, the strategy can be significantly improved.

Evolutionary Networks
An evolutionary network N is a tensor of dimension 4 of nodes M i which are networks themselves with edges drawn if M k evolves into M k and weight corresponding to the number of times that a network M k has evolved into M k . Fig. 9 shows a subnetwork of the full network for each evolutionary strategy from 50 (pseudo-)randomly generated networks with the biological ERBB signalling network as target.
Mutations and overexpression of ERB receptors (ERBB2 and ERBB3 in this network) have been strongly associated to more than 10 types of tissue-specific cancers and they can be seen at the highest level egulating most of the acyclic network.
We call forward mutations, mutations that led to the target network, and backward mutations, mutations that get away from the target network through the same evolutionary paths induced by forward mutations.
The forward mutations in the neighbourhood of the target (the evolved ERBB network) for each strategy, are as follow.
One of the mutations by BDM involves the breaking of the only network cycle of size 6: EGFR→ERBB3, ERBB3→IGFR, IGFR→ESR, ESR→MYC, MYC→EGFR by deletion of the interaction MYC→EGFR, with probability 0.05 among the possible mutations in the BDM immediate neighbourhood of the target. In the cycle is involved ERBB3 which has been found to be related to many types of cancer when overexpressed [29] Fig. 9. Evolutionary convergence in evolutionary subnetworks closest to the target matrix with edges shown only if they were used more than once. The darker the edge colour the more times that an evolutionary path was used to reach the target, the highest number is 7 for the BDMbased mutation (top) and the lower 2 (e.g. all those in the uniform random distribution). BDM-based is the only disconnected graph meaning that it produced a case of convergent evolution even before reaching the target.
For the local BDM strategy, the following were the top 5 forward mutations: EGFR→ERBB2, 0.32; EGFR→ ERBB3, 0.107; IGFR→ CCNE, 0.0714; ERBB3→ ERBB2, 0.0714; EGFR→ ESR, 0.0714; with ERBB2 and ERBB3 heavily involved in 3 of the top 5 possible mutations with added probability 0.49 and thus more likely than any other pair and interaction of proteins in the network.
Under a hypothesis that mutations can be reversals to states in past evolutionary pathways, then mutations to such interactions may be the most likely backward mutations to occur.

The Case for Localized Mutations and Modularity
As previously mentioned in the proposition 3, the main causes of slow-down under the BDM distribution are maladaptive persistent structures. These structures will negatively impact the evolution speed in factorial order relative to the size of the state space. One direct way to reduce the size set of possible mutations is to reduce the size of the matrices we are evolving. However, doing so will reduce the number of interesting objects we can evolve towards too. Another way to accomplish the objective while using the same heuristic is to rely on localized (or modular) mutations. That is, we force the mutation to take place on a submatrix of the input matrix.
The way we implement the stated change is by adding a single step in our evolution dynamics: at each iteration, we will randomly draw, with uniform probability, one submatrix of size 4 × 4 out of the set of adjacent submatrices that compose the input matrix, with no overlap, and force the mutation to be there by computing the probability distribution over all the matrices that contain the bit-shift only at the chosen place. We will call this method the local BDM method.
It is important to note that, within 1-bit shifts (point mutations), the space of total possible mutations remains the same when compared to the uniform and BDM strategies. Furthermore, the behaviour of the uniform strategy would remain unchanged if the extra step is applied using the Uniform distribution.
We repeated the experiment shown in the table 3 with the addition of the local BDM strategy and the same 50 random initial matrices. Its results are shown in the table 4. As we can see from the results obtained, local BDM obtains a statistically significant speed-up of 1.25 when compared to the uniform. One potential explanation of why we failed to obtain speed-up for the network with the BDM strategy is that, as an approximation to K, the model depends on finding global algorithmic structures, while the sample is based on a substructure which might not have enough information about the underlying structures we hypothesize govern the natural world and make possible scientific models and predictions.
However, biology evolves modular systems [31], such as genes and cells, that in turn build building blocks such as proteins and tissues. Therefore, local algorithmic mutation is a better model. This is a good place to recall that local BDM was devised as a natural solution to the problem presented by maladaptive persistent structures in global algorithmic mutation. Which also means that this type of modularity can be evolved by itself given that it provides an evolutionary advantage, as our results demonstrate.
We will further explore the relationship between BDM and local BDM within the context of global structures in the next section. Our current setup is not optimal for further experimentation in biological and local structured matrices, as the computational resources required to build the probability distribution for each instance grows in quadratic order relative to matrix size, though these computational resources are not needed in the real world (c.f. Conclusions).

Chasing Synthetic Evolving Networks
The aim of the next set of experiments was to follow, or chase, the evolution of a moving target using our evolutionary strategies. In this case, we chased 4 different dynamical networks: the ZK graphs [20], K-ary trees, an evolving Nstar graph and a star-to-path graph dynamic transition artificially created for this project (see Appendix for code). These dynamical networks are families of directed labelled graphs that evolve over time using a deterministic algorithm, some of which display interesting graph-theoretic and entropy-fooling properties [20]. As the evolution dynamics of these graphs are fully deterministic, we expected BDM to be significantly faster than the other two evolutionary strategies, uniform probability and local BDM.
We chased these dynamics in the following way: Let S 0 , S 1 , . . ., S n , . . . be the stages of the system we are chasing. Then the initial state S 0 was represented by a random matrix and, for each evolution S i → S i+1 , the input was defined as the adjacency matrix corresponding to S i , while the target was set as the adjacency matrix for S i+1 . In order to normalize the matrix size, we defined the networks as always containing the same number of nodes (16 for 16 × 16 matrices). We followed each dynamic until the corresponding stage could not be defined in 16 nodes.
The results which were obtained, starting from 100 random graphs and 100 different evolution paths at each stage, are shown in Figure 10. It is important to note that, since the graphs were directed, the matrices used were nonsymmetrical. From the results we can see that local BDM consistently outperformed the uniform probability evolution, but the BDM strategy was the faster by a significant margin. The results are as expected and confirm our hypothesis: uniform evolution cannot detect any underlying algorithmic cause of evolution, while BDM can, inducing a faster overall evolution. Local BDM can only detect local regularities, which is good enough to outrun uniform BDM in these cases. However, as the algorithmic regularities are global, local BDM is slower than (global) BDM.

DISCUSSION AND CONCLUSIONS
The mathematical mechanisms of biological information, from Mendelian inheritance to Darwin's evolution and the discovery of the digital nature of the genetic code together with the mechanistic nature of translation, transcription and other inter cellular processes, not only suggest a strong computational basis from algorithmic principles underlying biological processes, but by taking it to the ultimate level, they suggest that evolution by natural selection may not be (very) different, and can thus be regarded and studied as evolving programs in software space in the terms suggested by Chaitin [6], [7] and by Wolfram [3].
Our findings demonstrate that key aspects of life may be better explained with this algorithmic account and that the empirical observation of the rate of biological evolutionary convergence emerges naturally, unlike the alternative-and generally assumed-hypothesis of uniformly distributed random mutations. Here we have answered some questions by abandoning this uniform distribution assumption, questions ranging from the apparition of sudden major stages of evolution, the emergence of 'subroutines' in the form of modular persistent structures and the need of an evolving memory carrying information organized in such modules that drive the speed-up in the process of evolution by selection.
The algorithmic distribution emerges naturally from the interaction of deterministic systems [10], [24], [25], [32]. In other words, we are simulating the conditions of an algorithmic/procedural world and there is no reason to believe that it requires greater real-world (thus highly parallel) computation than is required by the assumption of the uniform distribution. The Universal Distribution thus emerges naturally and is the null hypothesis in algorithmic space in which noise is only apparent in nature and a consequence of other interacting open deterministic systems in, e.g., the organism's environment from which it harnesses algorithmic structure to its own advantage [4].
The interplay of the evolvability of organisms from the persistence of such structures also explains two opposed phenomena: recurrent explosions of diversity and mass extinctions, phenomena which have occurred during the history of life on earth that have not been satisfactorily explained under the uniform mutation assumption. However, taking the informational and computational aspects of life based on modern synthesis to the ultimate and natural consequences, the present approach based on weak assumptions of deterministic dynamic systems offers a novel framework of algorithmic evolution within which to study both biological and artificial evolution.

APPENDIX A APPROXIMATIONS TO ALGORITHMIC COMPLEXITY
The algorithmic complexity of a string K(s) (also known as Kolmogorov-Chaitin complexity [16], [17]) is defined as the length of the smallest program that produces s as an output and halts. This measure of complexity is invariant-up to a constant value-with respect to the choice of reference universal Turing machine. Finding the exact value of K(s) for any s is a lower semi-computable problem. This means that there is no general effective method to find K(s) for any given string, but upper bounds can be estimated.
The CTM relies upon approximating the algorithmic probability of an object by running every possible a large set of small Turing machines, generating an empirical probability distribution for the produced strings by counting the number of small Turing machines that produce each string and halt. The algorithm can only be decided for a small number of Turing machines and for those that can be decide it runs in exponential time, therefore only approximations of K(s) for small strings are feasible. However, this computation only needs to be done once to populate a lookup table that allows its application in linear (constant in exchange of memory) time.
BDM is an extension of CTM defined as where each s i corresponds to a substring of s for which its CTM value is known and n i is the number of times the string s i appears in s. A thorough discussion of BDM is found in [19].

APPENDIX B RECURSIVELY GENERATED GRAPHS
To test the speed of algorithmic evolution on recursive dynamic networks we generated 3 other low algorithmic graphs different from the ZK graph as defined in [20] that is also of low algorithmic complexity. We needed graphs that evolved over time in a low algorithmic complexity fashion from and to low algorithmic complexity graphs. The 3 graphs were canonically labelled using the positive natural numbers up to n by maximizing the number of nodes with consecutive numbers, then a rule was applied from lowest to highest number until the transformation was complete. The Wolfram Language code used to generate these recursively (hence of low algorithmic complexity/randomness) evolving graphs are the following.For the ZK graph [20] recursively generated by: where n is the size of the star graph in the 2 previous evolving graphs.
The K-ary tree evolving graph was generated with the Wolfram Language built-in function KaryTree[n], where n is the size of the K-ary tree.

PLE GRAPHS
Details of the theoretical and numerical application of BDM to matrices and graphs are provided in in [26], [34].
The tables 5, 6 and 7 contain full statistical information for the speed-up obtained for the simple graphs 'complete', 'star' and 'grid'.