# Reconstructing Attractors with Autoencoders

F. Fainstein  Universidad de Buenos Aires, Facultad de Ciencias Exactas y
Naturales, Departamento de Física, Ciudad Universitaria, 1428 Buenos Aires,
Argentina.  CONICET - Universidad de Buenos Aires, Instituto de Física
Interdisciplinaria y Aplicada (INFINA), Ciudad Universitaria, 1428 Buenos
Aires, Argentina.  G. B. Mindlin  Universidad de Buenos Aires, Facultad de
Ciencias Exactas y Naturales, Departamento de Física, Ciudad Universitaria,
1428 Buenos Aires, Argentina.  CONICET - Universidad de Buenos Aires,
Instituto de Física Interdisciplinaria y Aplicada (INFINA), Ciudad
Universitaria, 1428 Buenos Aires, Argentina.  P. Groisman<EMAIL_ADDRESS>IMAS-CONICET and Departamento de Matemática, Facultad de Ciencias Exactas y
Naturales, Universidad de Buenos Aires, Buenos Aires, Argentina.

###### Abstract

We show how to train an autoencoder to reconstruct an attractor from recorded
footage, preserving the topology of the underlying phase space. This is
explicitly demonstrated for the classic finite-amplitude Lorenz atmospheric
convection problem.

Introduction.— In the natural sciences, many problems can be modeled in terms
of dynamic systems. For those problems, the state of the system at each time
corresponds to a point in a finite dimensional space, usually called the phase
space. The temporal evolution is represented by a trajectory contained on a
manifold $\mathcal{M}$. Elucidating the dynamics of the problem typically
involves characterizing the attractor (which represents the long-time
behavior) from the available data. This characterization can be done computing
topological invariants, which act as fingerprints of the system [1, 2].
However, a common challenge arises on how to reconstruct the attractor when we
have access to a subset of the variables or, as in a video recording, when we
do have access to a large set of measurements but the relevant variables are
given by an unknown nonlinear transformation.

Currently, numerous approaches leverage the breakthrough of artificial
intelligence to propose algorithms aimed at overcoming this challenge (see,
for instance, Refs. [3, 4, 5, 6, 7, 8, 9]). The question arises if the
reconstructed attractors preserve the topology of the original phase space of
the problem. In this Letter we address this issue using an autoencoder, a
prominent architecture in the realm of artificial neural networks [10]. We
focus on the case in which the data has been measured by an array of sensors
that provide a temporally ordered sequence of frames capturing the
spatiotemporal evolution of the system, and we discuss simple extensions for
other cases. Through a formal mathematical argument, we propose a training
loss function that guarantees that the reconstructed flow preserves the
topology of the underlying manifold. This is explicitly demonstrated in the
classic Lorenz model for finite-amplitude atmospheric convection [11].

The autoencoder.— An autoencoder is a type of artificial neural network
specifically designed for dimensional reduction. It comprises an encoder
($E_{w}$) and a decoder ($D_{w}$). Both functions are implemented through
layers of units (“neurons”), each one providing an output which is a nonlinear
function of a linear combination of its inputs. The inputs are weighted sums
of the outputs from the previous layer. The encoder maps the input data into a
lower-dimensional space, and correspondingly, the decoder reconstructs the
input from this encoded representation. The mid layer of the autoencoder,
often referred to as the latent space, serves as a compressed representation
of the data, where the most critical information about the input is condensed.
Autoencoders can be trained with segments of time traces, sets of arbitrarily
chosen features, or, as in our case, snapshots of our spatiotemporal data.

The training of a typical autoencoder proceeds as follows. Each frame is a set
of $d_{\mathbf{i}}$ numbers (the number of pixels of each frame) and the
weights in the network are chosen in a way that the sum of the differences
between outputs and inputs, added over the whole dataset, is as small as
possible. The interesting aspect of the method is that the number of units in
the middle layer $d_{\mathbf{l}}<d_{\mathbf{i}}$ is chosen as small as
possible. If the compression is pushed too much, there will be loss of
information and the decoder will not be able to assign a unique output to each
input frame. In other words, the successful recovery of a unique output for
each input requires a unique representation of each element of the input in
the latent space.

Carrying out this type of training, it was observed both by processing
synthetic movies and others from actual experiments, that the encoder induces
trajectories in the latent space whose topological organization is equivalent
to the ones of the original phase space of the problem being filmed [12, 13].
The reasons (if any) that guarantee this property to persist in general
situations remain veiled, and it is the purpose of this Letter to propose a
framework to solve this problem.

Formal problem statement and proposal.— Let us assume that the dynamics of our
phenomenon is described by a dynamical system

$\dot{y}(t)=\phi\left(y(t)\right),\qquad y\in\mathcal{M}.$ (1)

We do not have direct access to $y$ but our data is obtained through the
application of a function $\alpha\in
C^{2}(\mathcal{M}\times\mathcal{N},\mathbb{R}^{d_{\mathbf{i}}})$ such that for
every fixed $p$ in the manifold $\mathcal{N}$, $\alpha_{p}=\alpha(\cdot,p)$ is
an embedding. The variable $p$ represents all the additional data that is not
relevant for the problem but is saved when recording the movie. We will assume
the existence of a function $\beta$ such that for every $p$ and $y$,
$\beta(\alpha(y,p))=y$. This means that the video has been recorded in such a
way that all the relevant variables can be recovered from the movie frame by
frame, and $\beta$ is precisely the function that does so. For a given fixed
$p$, we have $\beta^{-1}=\alpha_{p}$. We do not have access neither to
$\alpha$ nor $\beta$. We just observe $x(t)=\alpha_{p}(y(t))$ (but not $y(t)$)
for a given background $p$ and trajectory $y(t)$.

The flow of a vector field $\phi$ on $\mathcal{M}$ consists on the
transformations $g_{\phi}\colon\mathcal{M}\to\mathcal{M}$ converting every
initial condition $y_{0}$ of Eq. 1 at time $0$ on the value $g_{\phi}y_{0}$ of
the solution at time $t$ [14]. The objective of our work is to build a low
dimensional representation of the flow from the movie (that we shall call
$g_{\varphi}$), ensuring that it is topologically equivalent to the (unknown)
flow ($g_{\phi}$). The proposal will be the representation obtained in the
latent space of a specifically trained autoencoder.

Two flows $g_{\varphi}\colon\mathcal{M}_{1}\to\mathcal{M}_{1}$ and
$g_{\Phi}\colon\mathcal{M}_{2}\to\mathcal{M}_{2}$ are topologically equivalent
if there is a homeomorphism $h\colon\mathcal{M}_{1}\to\mathcal{M}_{2}$ mapping
orbits of $g_{\varphi}$ to orbits of $g_{\Phi}$, preserving orientation of the
orbits [14]. Note that from our assumptions the movie and the original
phenomenon are topologically equivalent so we can work directly with the
recorded data. This means that if we want to ensure that $h=E_{w}$ provides a
flow in the latent space topologically equivalent to $g_{\varphi}$, we need to
show that: (1) the encoder $E_{w}$ is an homeomorphism, (2) the image of the
trajectories in the input space are the solutions of a dynamical system (ruled
by a vector field $\Phi$).

We consider that our data is dense enough, i.e. either our dynamics support a
strange attractor so that we are sampling an orbit that is dense or we have
samples of several trajectories that fill up $\mathcal{M}$. We will assume
that we are in the first case and consider only one orbit, $x_{i}=x(t_{i})$,
$t_{i}\in[0,T]$, with the labels ordered in time (i.e.: $t_{i+1}>t_{i}$ for
every $i$). If we are in the second case we can proceed similarly but taking
several trajectories instead of only one. We will guide the training with the
loss function
$\mathcal{L}=\lambda_{1}\mathcal{L}_{1}+\lambda_{2}\mathcal{L}_{2}$, composed
by the usual mean squared error

$\mathcal{L}_{1}(w)=\sum_{x_{i}}\left|\left(D_{w}\circ
E_{w}\right)\left(x_{i}\right)-x_{i}\right|^{2},$ (2)

and a second term given by

$\mathcal{L}_{2}(w)=\sum_{x_{i}}|(D_{w}\circ E_{w})(x_{i+1})-(D_{w}\circ
E_{w})(x_{i})\\\ -(x_{i+1}-x_{i})|^{2}.$ (3)

Here $w$ represents the weights of the network, the sum is over all the
training set, and $\lambda_{1}$ and $\lambda_{2}$ are scaling factors. We will
assume that the autoencoder generalizes, that means that $\mathcal{L}$ is
small not just for the training set, but in all the input space
$\alpha_{p}(\mathcal{M})\in\mathbb{R}^{d_{\mathbf{i}}}$ (for a discussion on
this assumption see [15, 16, 17]). This forces $E_{w}$ and $D_{w}$ to be
continuous, and this means that $E_{w}$ is, asymptotically, an homeomorphism
(with inverse $D_{w}$).

What is not yet guaranteed is that the trajectories in the latent space
$z(t)=E_{w}\left(x(t)\right)$ are in fact solutions of a dynamical system. The
second term of our loss, $\mathcal{L}_{2}$, has been added precisely to
enforce this. To see this, assuming our autoencoder is expresive enough, let
us compute the vector field in the latent space by taking the directional
derivative:

$\Phi\left(z\right)=\frac{\partial
E_{w}}{\partial\varphi\left(x\right)}\left(x\right)\approx\frac{\partial
E_{w}}{\partial\varphi\left(D_{w}\left(z\right)\right)}\left(D_{w}\left(z\right)\right).$
(4)

The second term states that

$\frac{\partial\left(D_{w}\circ
E_{w}\right)}{\partial\varphi\left(x\right)}\left(x\right)=\frac{\partial
D_{w}}{\partial\Phi\left(z\right)}\left(z\right)\approx\varphi\left(x\right).$
(5)

This condition guarantees that nonzero velocities in the space of frames are
not mapped into zero velocities in the latent space, that is $\Phi(z)\neq 0$
if $\varphi(x)\neq 0$. Note that it also penalizes strong spatial variations
of $\Phi$ in the latent space, as one would expect for a smooth enough vector
field defining a dynamical system. Therefore, in the limit as the size of the
network goes to infinity, the phase flows $g_{\varphi}$ and $g_{\Phi}$ are
topologically equivalent: the homeomorphism $E_{w}$ maps orbits of
$\dot{x}=\varphi(x)$ to orbits of $\dot{z}=\Phi(z)$ homeomorphically and
preserving orientation of the orbits.

Numerical work.— To numerically test our method, we generate a synthetic movie
motivated by the classic model developed by Lorenz for an atmospheric
convection problem [11]. Considering a layer of fluid of uniform depth $H$ and
aspect ratio $a$, Lorenz proposes a modal decomposition for the stream
function ($\psi$) and the departure of temperature from the non-convective
state ($\theta$). He writes

$\psi=c_{1}X(t)\sin\left(\frac{\pi ax}{H}\right)\sin\left(\frac{\pi
z}{H}\right),$ (6) $\theta=c_{2}Y(t)\cos\left(\frac{\pi
ax}{H}\right)\sin\left(\frac{\pi z}{H}\right)-c_{3}Z(t)\sin\left(\frac{2\pi
z}{H}\right),$ (7)

where $X$, $Y$, $Z$ are functions of time alone, and $c_{1}$, $c_{2}$ and
$c_{3}$ are constants. Here $x$ and $z$ denote the spatial coordinates of the
problem. The dynamics of the modal amplitudes is determined by the famous
dynamical system

$\displaystyle\dot{X}$ $\displaystyle=\sigma\left(Y-X\right)$ (8)
$\displaystyle\dot{Y}$ $\displaystyle={rX}-Y-{XZ}$ $\displaystyle\dot{Z}$
$\displaystyle={XY}-{bZ}.$

By numerically integrating this system with $\sigma=10$, $b=8/3$ and $r=28$,
we generate temporal series with $40,000$ points (time step of $0.01$). We
discard the first $1000$ points to avoid the transient state. Setting $a=H=1$,
$c_{1}=c_{2}=2c_{3}=1/20$, and a spatial discretization of $0.025$, we use
these time traces and Eqs. 6-7 to generate a sequence of $40\times 40$ pixel
frames describing the spatiotemporal pattern for $\psi$ and $\theta$.

To choose the proper training data set for our autoencoder, we note that not
every movie is an embedding. For instance, the movie of a moving mass attached
to a spring it is not. However, if we place a flag attached to the mass, it
is, since the velocity could be recovered from each frame. In our case, we
need information of the stream function and the temperature, so we train an
autoencoder with a data set composed of inputs that contain both frames (see
Fig. 1(a)). The first 30000 frames are used for training, and the remaining
9000 frames for testing.

Figure 1: (a) The autoencoder. Input data is encoded to and decoded from the
latent space. (b) Evolution of both terms of our loss function, measured in
the train and test datasets, left axis. Mean values over $20$ trained
networks, with their standard error are shown. Evolution of the percentage of
autoencoders with correct topology in the latent space as a function of the
training epoch is shown in black dashed line, right axis. (c) Reconstructed
flow in the latent space of an autoencoder. Five Lorenz unstable orbits and
their linking numbers in phase space (d) and latent space (e). Dot color
indicates which orbit is above at the crossing point.

The network architecture consists of initial and final fully connected layers
of $40\times 40\times 2$ units. Intermediate fully connected layers containing
$64$, $32$, and $16$ units encode the input to a three-dimensional latent
space. Sine activation functions are used for all the layers except the latent
and the output layers, that have no activation functions. The training process
is conducted with a batch size of $600$ samples and $10^{3}$ epochs. The
network is implemented in Python using PyTorch library (version 2.1.1). Adam
optimization algorithm with learning rate $10^{-4}$, $\beta_{1}=0.9$ and
$\beta_{2}=0.999$ is used, and scaling factors for the loss function are set
to $\lambda_{1}=1$, $\lambda_{2}=50$ to ensure that both terms have similar
magnitude. We do not perform any systematic search in the hyperparameter space
to optimize the loss on the test set.

In Fig. 1(b) we show the average evolution, over $20$ trained networks, of
both terms of the cost function. In Fig. 1(c) the reconstruction obtained in
the latent space for one of the networks is shown. We remark that Fig 1(c) is
not the famous attractor given by Lorenz equations but the one reconstructed
by the autoencoder.

To make a topological description, we use approximations of the unstable
periodic orbits coexisting with the strange attractor. The way in which these
orbits are intertwined is a fingerprint of the system. The organization of the
orbits can be algebraically described by their relative rotation rates and
linking numbers [18]. For two orbits of period one, one of period two, and two
of period three, selected from the test set, we compute the relative linking
numbers [19]. In Figure 1 (d)-(e) we show the result in the original phase
space and in the latent space, respectively. The evolution of the percentage
of networks with correct topology during the training epochs is shown on the
right axis of Figure 1(b). All networks learnt the correct topology.

Conclusions.— In this study, we propose a method based on autoencoders to
reconstruct the dynamics of a system from recorded footage. The method is
applicable to any dataset which is an embedding, for instance, long-enough
segments of time series data [20]. Recent research suggested that such
processing leads to a flow exhibiting a topological structure equivalent to
the dynamical system governing the behavior captured by the movie [13]. Here,
we showcase a modification to the loss function commonly used to guide the
parameter tuning of an autoencoder, resulting in a trained network capable of
reconstructing a flow topologically equivalent to the original one. It is not
guaranteed to be an embedding: we are not tuning the encoder to have a
differential that is invertible at any point but to be injective only in the
direction of the flow. Yet, it turns out that this is enough to recover every
topological property of the system.

Without this modification in the loss functions the autoencoder might still be
able to produce a topologically equivalent dynamical system in the latent
space, in a significant number of fitting runs. This might occur thanks to
some kind of inducted bias/implicit regularization [15, 16, 17] but it is not
clear yet what kind of regularization might take place in particular runs. By
adding the second term we force the autoencoder to have the regularization we
need to ensure topological equivalence.

Recently, autoencoders were used to derive dynamical systems that could serve
as models for data [5]. The efforts provided algorithms for reconstructing
vector fields in the latent space. Our work builds confidence on those
efforts, by showing that properly trained, the flows in the latent space are
topologically equivalent to the dynamics responsible for the data.

Over recent years, machine learning techniques have provided tools capable of
solving extraordinarily complex problems, including the prediction of
dynamics. Despite its apparent lack of necessity to illuminate the underlying
mechanisms behind the learned dynamics, this work illustrates how the depth of
the dynamics is embedded in some representation of the network. This study
aligns with a series of observations made across various network
architectures, and it would not be surprising if similar representation
mechanisms could be formulated in those cases as well. This research
contributes to our understanding of the nuanced interplay between artificial
intelligence and dynamic systems, opening avenues for further exploration and
application across diverse architectural structures.

## References

  * Guckenheimer and Holmes [2013] J. Guckenheimer and P. Holmes,  _Nonlinear oscillations, dynamical systems, and bifurcations of vector fields_ , Vol. 42 (Springer Science & Business Media, 2013).
  * Solari  _et al._ [1996] H. G. Solari, M. A. Natiello, and G. B. Mindlin,  _Nonlinear dynamics: a two-way trip from physics to math_ (Taylor & Francis Group, 1996).
  * Lu  _et al._ [2018] Z. Lu, B. R. Hunt, and E. Ott, Attractor reconstruction by machine learning, Chaos 28 (2018).
  * Brunton and Kutz [2019] S. L. Brunton and J. N. Kutz,  _Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control_ (Cambridge University Press, Cambridge, England, 2019).
  * Champion  _et al._ [2019] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton, Data-driven discovery of coordinates and governing equations, Proc. Natl. Acad. Sci. USA 116, 22445 (2019).
  * Agostini [2020] L. Agostini, Exploration and prediction of fluid dynamical systems using auto-encoder technology, Physics of Fluids 32, https://doi.org/10.1063/5.0012906 (2020).
  * Liu and Tegmark [2021] Z. Liu and M. Tegmark, Machine learning conservation laws from trajectories, Physical Review Letters 126, 180604 (2021).
  * Floryan and Graham [2022] D. Floryan and M. D. Graham, Data-driven discovery of intrinsic dynamics, Nature Machine Intelligence 4, 1113 (2022).
  * Young and Graham [2023] C. D. Young and M. D. Graham, Deep learning delay coordinate dynamics for chaotic attractors from partial observable data, Physical Review E 107, 034215 (2023).
  * Hinton and Salakhutdinov [2006] G. E. Hinton and R. R. Salakhutdinov, Reducing the dimensionality of data with neural networks, Science 313, 504 (2006).
  * Lorenz [1963] E. N. Lorenz, Deterministic nonperiodic flow, Journal of atmospheric sciences 20, 130 (1963).
  * Uribarri and Mindlin [2020] G. Uribarri and G. B. Mindlin, The structure of reconstructed flows in latent spaces, Chaos 30, 093109, 8 (2020).
  * Fainstein  _et al._ [2023] F. Fainstein, J. Catoni, C. P. Elemans, and G. B. Mindlin, The reconstruction of flows from spatiotemporal data by autoencoders, Chaos, Solitons & Fractals 176, 114115 (2023).
  * Arnold [2012] V. I. Arnold,  _Geometrical methods in the theory of ordinary differential equations_ , Vol. 250 (Springer Science & Business Media, 2012).
  * Belkin  _et al._ [2019] M. Belkin, D. Hsu, S. Ma, and S. Mandal, Reconciling modern machine-learning practice and the classical bias-variance trade-off, Proc. Natl. Acad. Sci. USA 116, 15849 (2019).
  * Radhakrishnan  _et al._ [2020] A. Radhakrishnan, M. Belkin, and C. Uhler, Overparameterized neural networks implement associative memory, Proc. Natl. Acad. Sci. USA 117, 27163 (2020).
  * Belkin [2021] M. Belkin, Fit without fear: remarkable mathematical phenomena of deep learning through the prism of interpolation, Acta Numerica 30, 203 (2021).
  * Mindlin  _et al._ [1990] G. B. Mindlin, X.-J. Hou, H. G. Solari, R. Gilmore, and N. Tufillaro, Classification of strange attractors by integers, Physical Review Letters 64, 2350 (1990).
  * Birman and Williams [1983] J. S. Birman and R. F. Williams, Knotted periodic orbits in dynamical systems i: Lorenz’s equations, Topology 22, 47 (1983).
  * Takens [2006] F. Takens, Detecting strange attractors in turbulence, in  _Dynamical Systems and Turbulence, Warwick 1980_ (Springer, 2006) pp. 366–381.


# DI2: prior-free and multi-item discretization of biomedical data and its
applications

LALeonardo Alexandre  RSCRafael S. Costa  RHRui Henriques  IDMEC, Instituto
Superior Técnico, Universidade de Lisboa, Portugal  INESC-ID, Portugal
Instituto Superior Técnico, Universidade de Lisboa, Portugal  LAQV-REQUIMTE,
DQ, Faculty of Science and Technology, NOVA University of Lisbon, Campus de
Caparica, 2829-516 Caparica, Portugal

###### Abstract

Motivation: A considerable number of data mining approaches for biomedical
data analysis, including state-of-the-art associative models, require a form
of data discretization. Although diverse discretization approaches have been
proposed, they generally work under a strict set of statistical assumptions
which are arguably insufficient to handle the diversity and heterogeneity of
clinical and molecular variables within a given dataset. In addition, although
an increasing number of symbolic approaches in bioinformatics are able to
assign multiple items to values occurring near discretization boundaries for
superior robustness, there are no reference principles on how to perform
multi-item discretizations.
Results: In this study, an unsupervised discretization method, DI2, for
variables with arbitrarily skewed distributions is proposed. DI2 provides
robust guarantees of generalization by placing data corrections using the
Kolmogorov-Smirnov test before statistically fitting distribution candidates.
DI2 further supports multi-item assignments. Results gathered from biomedical
data show its relevance to improve classic discretization choices.
Software: available at https://github.com/JupitersMight/DI2

multi-item discretization,

prior-free discretization,

heterogeneous biomedical data,

###### keywords:

Research

##  1 Introduction

Approaches to discretization of continuous variables have long been discussed
alongside their pros and cons. Altman et al. [1] and Bennette et al. [2] both
discuss the relevance and impact of categorizing continuous variables and
reducing the cardinality of categorical variables. Liao et al. [3] compares
various categorization techniques in the context of classification tasks in
medical domains, without using domain knowledge of field experts. The
relevance of discretization meets both descriptive and predictive ends,
encompassing state-of-the-art approaches such as pattern-based biclustering
[4] and associative models such as XGBoost [5].

In this work we present DI2, a Python library that extends non-parametric
tests to find the best fitting distribution for a given variable and
discretize it accordingly. DI2 offers three major contributions: i)
corrections to the empirical distribution before statistical fitting to
guarantee a more robust approximation of candidate distributions; ii)
efficient statistical fitting of 100 state-of-the-art theoretical
distributions; and, finally, iii) assignment of multiple items according to
the proximity of values to the boundaries of discretization, a possibility
supported by numerous symbolic approaches [6, 4, 7]. The assignment of
multiple items [8, 9], generally referred as multi-item discretization,
conferes the possibility to avail the wealth of data structures and algorithms
from the text processing and bioinformatics communities without the risks of
the well-studied item-boundaries problem.

##  2 Background

Discretization methods have wide taxonomy[10] with a determinant division in:
1) supervised, where the method uses the class variable to bin the data, and,
2) unsupervised, where the method is independent of the class variable. DI2
places itself on the latter, it works independently on the class variable.
Other characteristics of DI2 are: 1) static, where discretization of the
variables takes place prior to an algorithm; 2) global, uses information about
the variable as a whole to make the partitions and can still be applied with a
scarce number of observations; 3) direct and splitting, splits the whole the
range of values into k intervals simultaneously; and 4) multivariate and
univariate, DI2 can use either the whole dataset to create the intervals and
discretize each variable or use each variable individually to create the
respective intervals.

Some examples of unsupervised discretization methods are PD (Proportional
Discretization), FFD (Fixed Frequency Discretization)[11], equal-
width/frequency, k-means [12]. In this work, DI2 is compared with such classic
discretization methods. These are illustrated in Figures I, II, and III.

Figure I: Illustration of equal-frequency method with 9 points along an axis
and 3 categories. This method is based on the frequency of the items, where
each category has the same number of items, in order to set the intervals.
Figure II: Illustration of equal-width method with 9 points along an axis and
3 categories. This method is based on the range taken by the items, where each
category has the same width interval. Figure III: Illustration of K-means
method with 9 points along an axis and 3 categories. This method is based in
the k-means clustering, where each category is defined by a centroid.

###  2.1 DI2: normalization and feature scaling

While not mandatory, DI2 supports the following techniques: 1) min-max
scaling,

$X^{\prime}=\frac{X-X_{min}}{X_{max}-X_{min}},$ (1)

where $X$ is an ordered set of observed values; 2) z-score normalization,
usually applied to samples normally distributed [13],

$X^{\prime}=\frac{X-\mu}{\sigma},$ (2)

where $X$ is an ordered set of observed values; and 3) mean normalization,

$X^{\prime}=\frac{X-\mu}{X_{max}-X_{min}}.$ (3)

###  2.2 DI2: statistical hypotheses

In order to discretize the data into intervals, DI2 provides two statistical
hypothesis tests to be applied: 1) $\tilde{\chi}^{2}$ test [14], and 2)
Kolmogorov-Smirnov goodness-of-fit test [15].

In the aforementioned tests the observed distribution is matched with a
theoretical continuous
distribution111https://docs.scipy.org/doc/scipy/reference/stats.html provided
by the SciPy open-source library [16] where the parameters are estimated
through maximum likelihood estimation. The binning of the distributions for
the $\tilde{\chi}^{2}$ test is based on the number of categories the user
inputs and are built using the cumulative distribution function. The user can
either choose the $\tilde{\chi}^{2}$ or the Kolmogorov-Smirnov goodness-of-fit
as the primary fitting test. The theoretical continuous distribution with the
lowest test statistic is picked as the best fit for the observed distribution.

###  2.3 DI2: outlier correction

The Kolmogorov-Smirnov goodness-of-fit test can optionally be used to remove
up to 5% outlier points, from the observed distribution, according to the
theoretical continuous distribution under assessment. Kolmogorov-Smirnov
goodness-of-fit test returns a statistic (D statistic) that represents the
maximum distance between the observed that and the theoretical continuous
distribution we are testing.

$D=max\\{\max_{1\leq j\leq n}\\{\frac{j}{n}-F(X_{j})\\},\max_{1\leq j\leq
n}\\{F(X_{j})-\frac{(j-1)}{n}\\},$ (4)

where $n$ is the number of samples, $j$ is the index of a given sample, and
$F$ is the frequency of sample $X_{j}$. Using this statistic we can pin point
where the farthest point between the distributions is and remove it. After up
to 5% of the samples have been removed, the iteration with the best
Kolmogorov-Smirnov statistic is picked (from 0 outliers removed to up to 5%).
The data produced by outlier removal is then used to run the main statistical
hypothesis test picked ($\tilde{\chi}^{2}$ or Kolmogorov-Smirnov). This
correction guarantees the absence of penalizations caused by abrupt yet
spurious deviations driven by the selected histogram granularity and help
consolidate the choice of the theoretical continuous distribution.

###  2.4 DI2: multi-item discretization

After selecting the theoretical continuous distribution that best fits the
continuous variable, DI2 proceeds with the discretization. Given a desirable
number of categories (bins), multiple cut-off points are generated using the
inverse cumulative distribution function of the theoretical continuous
distribution. The cut-off points guarantee an approximately uniform
distribution of observation per category, although empirical-theoretical
distribution differences can underlie imbalances.

DI2 supports multi-item assignments by identifying border values for each
category, this is exemplified in Figure IV. To this end, the user can
optionally also define a percentage (between 0 and 50% with 20% default) to
affect the width of the borders. These borders take an intermediate value
which symbolize that it belongs to both upper and lower category. Width
extremes, 0% (50%) correspond to none (one) additional category assigned to
every observation.

Figure IV: Illustration example of discretization with 9 points along an axis
and 3 categories considering border values (values which belong to 2
categories).

##  3 Implementation

DI2 tool is fully implemented in Python and is provided as an open-source
method at GitHub with well-annotated APIs and notebook tutorials for a
practical illustration of its major functionalities. The algorithm workflow is
shown in Algorithm 1 and the Kolmogorov-Smirnov correction is shown in
Algorithm 2.

Input: dataset, number_of_bins

Optional input: statistical_test=”chi2”, multi_item_cutoff_margin=0.2,
kolmogorov_opt=True, normalizer=”min_max”, distributions=[…],
single_column_discretization=True

Output: The dataset discretized

y_normalized = [];

if  _single_column_discretization_ then

for  _column in dataset.columns_ do

y_normalized = normalization(dataset[column],normalizer);

main_operation(…);

end for

else

for  _column in dataset.columns_ do

y_normalized.append(normalization(dataset[column],normalizer));

end for

main_operation(…);

end if

Function  _main_operation(_…_)_:

best_dist, test_statistic, data_used;

for  _distribution in distributions_ do

temp_statistic;

if  _statistical_test == ”chi2”_ then

temp_data = kolmogorov_goodness_of_fit(y_normalized, distribution,
kolmogorov_opt);

temp_statistic = chi_squared_goodness_of_fit(temp_data, distribution,
number_of_bins);

else

temp_statistic, temp_data = kolmogorov_goodness_of_fit(y_normalized,
distribution, kolmogorov_opt);

end if

if  _temp_statistic $<$ test_statistic_ then

test_statistic = temp_statistic;

best_dist = distribution;

data_used = temp_data;

end if

end for

dataset[column] = discretize(best_dist, multi_item_cutoff_margin, data_used,
dataset[column], number_of_bins, y_normalized);

return

Algorithm 1 DI2 main algorithm

Input: observed_distribution, theoretical_distribution, outlier_removal_flag

Output: The statistic of Kolmogorov test and the corresponding data

N5 = size(data) $\times$ 0.05 if outlier_removal_flag else 1;

results = [];

i = 0;

while  _i $<$ N5_ do

Estimate_Parameters(theoretical_distribution);

D_plus = D_minus = [];

idx_max_d_plus = idx_max_d_minus = [];

calculate_d_minus(D_minus, idx_max_d_minus);

calculate_d_plus(D_plus, idx_max_d_plus);

if  _len(results) == 0_ then

results = [max(D_plus[idx_max_d_plus], D_minus[idx_max_d_minus]),
data.copy()];

else

ks = max(D_plus[idx_max_d_plus], D_minus[idx_max_d_minus]) if  _ks $<$
results[0]_ then

results = [ks, data.copy()];

end if

end if

if  _D_plus[idx_max_d_plus] $>$ D_minus[idx_max_d_minus]_ then

delete data[idx_max_d_plus];

else

delete data[idx_max_d_minus];

end if

++i;

end while

return results;

Algorithm 2 Kolmogorov outlier correction

##  4 Results and Discussion

To illustrate some of the DI2 properties, we considered two published
datasets: 1) the breast-tissue dataset [17], containing electrical impedance
measurements in samples of freshly excised tissue from the breast, and 2) the
yeast dataset [18], containing molecular statistics variables. Both of these
are available at the UCI Machine Learning repository [19].

This section first discusses results on breast-tissue dataset, DI2 is executed
with $\tilde{\chi}^{2}$ as the main statistical test, with and without
Kolmogorov outlier removal, with single column discretization, and 5
categories per variable outputted. We will then present and discuss the
fitting of the distribution with and without Kolmogorov outlier removal and
the interval borders created compared with equal-frequency and equal-width.

We then consider yeast dataset, DI2 is executed with $\tilde{\chi}^{2}$ as the
main statistical test, with and without Kolmogorov outlier removal, with
single column and whole dataset discretization, and 5 categories per variable
outputted, We will then present and discuss the distribution of the values by
category with and without border values and compare the different executions
of DI2 with themselves and with K-means discretization category distribution.

Finally, still considering the yeast dataset, we will present the execution of
multiple algorithms and the accuracy achieved with different DI2
discretization configurations and the other aforementioned discretization
techniques.

###  4.1 breast-tissue dataset

V.a. Q-Q plot of empirical distribution (blue dots) against the fitted
recipinvgauss distribution (red line).

V.b. Q-Q plot of empirical distribution (blue dots) against the fitted chi2
distribution (red line).

V.c. Empirical distribution (gray bins) and corresponding cut-off points using
equal-width, equal-frequency and D2I statistical fitting with and without
Kolmogorov-Smirnov correction. Red and yellow lines correspond to category and
border boundaries.

Figure V: Distribution matching of DA variable from breast-tissue againt two
statistical distributions (recipinvgauss in V.a and chi2 in V.b, as well as
the corresponding discretization boundaries and border values (V.c).

The breast-tissue dataset contains 106 data instances and 9 continuous
variables (I0, PA500, HFS, DA, AREA, A/DA, MAX IP, DR, P). The gathered
results show the decisions placed by DI2 in the absence and presence of
Kolmogorov-Smirnov optimization.

Table 1 shows the best fitting distribution for each continuous variable of
the dataset. Variables ”I0”, ”PA500”, ”A/DA”, ”DR”, and ”P” remained unchanged
with a removal of up to 5% of outlier points. Variables ”HFS” and ”Area”
produced better results in the $\tilde{\chi}^{2}$ test with the removal of
outliers solidifying the distribution choice. Finally, the fitting choice
changed for variables ”DA” and ”Max IP” under the $\tilde{\chi}^{2}$ test,
revealing a more solid choice from the analysis of the residuals.

Table 1: Best fitting distributions for each continuous variable, without and with Kolmogorov-Smirnov correction. Both $\tilde{\chi}^{2}$ (primary) and KS statistics are shown. Variables |  Without opt. |  $\tilde{\chi}^{2}$ statistic |  p-value <0.05 |  D statistic |  With opt. |  $\tilde{\chi}^{2}$ statistic |  p-value <0.05 |  D statistic
---|---|---|---|---|---|---|---|---
I0 |  alpha |  8.8 |  True |  0.12 |  alpha |  8.8 |  True |  0.11
PA500 |  exponnorm |  2.98 |  True |  0.07 |  exponnorm |  2.98 |  True |  0.07
HFS |  foldcauchy |  2.25 |  True |  0.07 |  foldcauchy |  1.57 |  True |  0.07
DA |  recipinvgauss |  1.6 |  True |  0.06 |  chi2 |  1.01 |  True |  0.06
Area |  frechet_r |  0.5 |  True |  0.07 |  frechet_r |  0.25 |  True |  0.05
A/DA |  mielke |  1.17 |  True |  0.06 |  mielke |  1.17 |  True |  0.05
Max IP |  johnsonsu |  4.72 |  True |  0.05 |  alpha |  1.09 |  True |  0.07
DR |  johnsonsb |  1.2 |  True |  0.05 |  johnsonsb |  1.2 |  True |  0.05
P |  genextreme |  5.13 |  True |  0.09 |  genextreme |  5.13 |  True |  0.09

Considering variable ”DA”, Figures V.a and V.b show its Q-Q (quantile-
quantile) plot, offering a view on the adequacy of the statistical fitting. In
this context, we depict histograms for the observed data with 100 bins (blue
dots) and the best theoretical distribution picked without and with
Kolmogorov-Smirnov correction (red line). A moderate improvement from Figure
V.a to V.b can be detected, with the observed quantiles (blue dots) being
closer to the theoretical continuous quantiles (red line).

After the fitting stage, cut-off points are calculated to produce the final
categories. Figure V.c compares different discretization options: equal-
frequency, equal-width, and the two best fitting theoretical continuous
distributions (without and with Kolmogorov-Smirnov optimization). Cut-off
points are marked as red lines, and the border cut-off points in yellow. This
analysis shows how critical discretization can be for determining the
inclusion or exclusion of high density bins. The ability of DI2 to assign
multiple items using borders can thus be explored by symbolic approaches to
mitigate vulnerabilities inherent to the discretization process [20, 21].

###  4.2 yeast dataset

The yeast dataset contains 1484 data instances and 10 variables, including the
sample identification, class, and 8 molecular statistics variables (mcg, gvh,
alm, mit, erl, pox, vac, nuc). In the breast-tissue dataset we compared DI2
with two other unsupervised discretization, equal-frequency and equal-width.
With this dataset we will compare DI2 with k-means.

Table 2 displays the results of the statistical tests produced by DI2 when
applied to each variable independently and, the last row of the table, when
applied to the whole dataset together. Let’s use variable ”mit” as an example
for this dataset. Figure VI.a displays the distribution of values in the
variable ”mit” before outlier removal (brown and blue area of histogram) and
after outlier removal (brown area of histogram).

Figures VI.b and VI.c compare the distribution of the categories of both
discretization techniques (DI2 with outlier removal along different
discretizations and k-means discretization), and also assess the impact of
outlier removal had in categorizing in different executions of DI2 and
k-means. Figure VII presents the border values under different DI2
discretization settings.

VI.a. Variable ”mit” distribution before Kolmogorov-Smirnov outlier removal
(brown and blue area) after Kolmogorov-Smirnov outlier removal (brown area) .

VI.b. Variable ”mit” categories distribution after DI2 single column
discretization, Kolmogorov-Smirnov outlier removal and no border values (blue
columns) and after k-means discretization (pink columns).

VI.c. Variable ”mit” categories distribution after DI2 discretization with
different settings without border values. Single column discretization without
Kolmogorov-Smirnov outlier removal (dark blue columns), whole dataset
discretization with Kolmogorov-Smirnov outlier removal (light purple columns),
whole discretization without Kolmogorov-Smirnov outlier removal (dark purple
columns).

Figure VI: Variable ”mit” distribution, IV.a., categories distribution after
DI2, IV.b. and IV.c., and after k-means discretization, IV.b.

VII.a. Categories between 0 and 2.

VII.b. Categories between 2 and 4.

Figure VII: Variable ”mit” categories distribution after DI2 discretization with different settings with border values. Single column discretization with Kolmogorov-Smirnov outlier removal (light blue columns), single column discretization without Kolmogorov-Smirnov outlier removal (dark blue columns), whole dataset discretization with Kolmogorov-Smirnov outlier removal (light purple columns), whole discretization without Kolmogorov-Smirnov outlier removal (dark purple columns). Table 2: Best fitting distributions for each continuous variable, without and with Kolmogorov- Smirnov outlier removal. Variables |  Without opt. |  $\tilde{\chi}^{2}$ statistic |  p-value <0.05 |  D statistic |  With opt. |  $\tilde{\chi}^{2}$ statistic |  p-value <0.05 |  D statistic
---|---|---|---|---|---|---|---|---
mcg |  foldcauchy |  3.72 |  True |  0.08 |  exponnorm |  3.18 |  True |  0.02
gvh |  genlogistic |  3.57 |  True |  0.03 |  genlogistic |  2.02 |  True |  0.02
alm |  genlogistic |  17.00 |  False |  0.05 |  genlogistic |  12.08 |  False |  0.03
mit |  exponnorm |  19.23 |  False |  0.05 |  exponnorm |  6.11 |  True |  0.03
erl |  chi2 |  4.45 $\times 10^{-14}$ |  True |  0.99 |  chi2 |  4.23 $\times 10^{-14}$ |  True |  0.99
pox |  chi2 |  4.45 $\times 10^{-14}$ |  True |  0.99 |  gengamma |  4.23 $\times 10^{-14}$ |  True |  0.99
vac |  laplace |  20.99 |  False |  0.08 |  pearson3 |  14.18 |  False |  1.00
nuc |  exponnorm |  1116.63 |  False |  0.26 |  mielke |  795.28 |  False |  0.26
all variables |  genhalflogistic |  45.69 |  False |  0.25 |  genhalflogistic |  10.25 |  False |  0.21

The performed analysis for the yeast dataset shows how critical the category
border, previously discussed in more detail with the breast-tissue dataset,
can be, determining the inclusion or exclusion of values. The ability of DI2
to assign multiple items using borders can be explored by symbolic approaches
to mitigate vulnerabilities inherent to the discretization process as
discussed in the following subsection.

###  4.3 Predictive performance

To test the impact DI2 has when discretizing data we considered the yeast
dataset, 5 categories per variable, and six supervised classification methods:
Naive Bayes [22], Random Forest [23], SMO [24](Sequential Minimal
Optimization), C4.5 [25], MLRM [26](Multinomial Logistic Regression Model),
and FleBiC [27].

In Figure VIII the results for the aforementioned models are presented, with
the exception of FleBiC which will be discussed later on. Each bar represents
the average accuracy achieved with each discretization method, and the small
bracket the standard deviation of the accuracy. In each model, DI2, with
configurations of single column discretization and outlier removal, matched
with the highest accuracy achieved by other discretization methods, and for
the C4.5 model, DI2, with configurations of combined column discretization,
achieved the highest accuracy compared with other discretization methods.

Figure VIII: Average accuracy per classifier and discretization method
available without border values. From left to right in each group of bars:
K-means, Quartile, Uniform, DI2 (single, kol. correction), DI2 (single), DI2
(whole, kol. correction), DI2 (whole).

In order to fully test out the potential of DI2, we now consider border
values. From the aforementioned supervised classification methods only FleBiC
is able to test this feature. FleBiC can be executed normally or with Random
Forests (weight of decision 50% FleBiC, 50% Random Forests), which we will
designate as FleBiC Hybrid. In Figure IX the results of executing FleBiC and
FleBiC Hybrid are presented. In terms of average accuracy, Figure IX.a., both
FleBiC and FleBiC Hybrid predict with a higher average accuracy when using DI2
method than when using other discretization methods. Within the different
settings of DI2, FleBiC Hybrid presents a higher accuracy when the predictive
model considers border values. Finally in terms of sensitivity, Figure IXb.,
we can see how considering border values affects the prediction of class NUC,
making it possible to break through a ceiling on the prediction of class NUC
when other discretization methods couldn’t.

IX.a. Average accuracy of each FleBiC with different discretization methods.

IX.b. Average sensitivity when predicting class NUC of each FleBiC version.

Figure IX: Accuracy when executing different FleBiC versions, and Sensitivity
of when predicting class NUC, with multiple discretization methods. From left
to right the bars are: K-means, Quartile, Uniform, DI2 (single, kol.
correction), DI2 (single), DI2 (border values, single, kol. correction), DI2
(border values, single).

##  5 Conclusion

This work proposed a new unsupervised method for data discretization, DI2,
that takes into account data distribution, outliers within the data, and
border values.

Results show that DI2 is a viable and coherent form of discretizing data. When
compared with other well-known and frequently used unsupervised methods, DI2
brings out the same average accuracy in well-known and supervised classfiers.
If we consider border values, which other unsupervised methods do not
consider, and use a classifier that handles border values, FleBiC, then DI2 is
able to improve the average accuracy previously achieved without border
values.

## Acknowlegments

This work was supported by Fundação para a Ciência e a Tecnologia (FCT),
through IDMEC, under LAETA project (UIDB/50021/2020), IPOscore with reference
(DSAIPA/DS/0042/2018), and ILU (DSAIPA/DS/0111/2018). This work was further
supported by the Associate Laboratory for Green Chemistry (LAQV), financed by
national funds from FCT/MCTES (UIDB/50006/2020 and UIDP/50006/2020), INESC-ID
plurianual (UIDB/50021/2020) and the contract CEECIND/01399/2017 to R. S.
Costa.

## References

  * [1] Altman, D.G.: Categorizing continuous variables. Wiley StatsRef: Statistics Reference Online (2014)
  * [2] Bennette, C., Vickers, A.: Against quantiles: categorization of continuous variables in epidemiologic research, and its discontents. BMC medical research methodology 12(1), 21 (2012)
  * [3] Liao, S.-C., Lee, I.-N.: Appropriate medical data categorization for data mining classification techniques. Medical informatics and the Internet in medicine 27(1), 59–67 (2002)
  * [4] Henriques, R., Madeira, S.C.: Bicpam: Pattern-based biclustering for biomedical data analysis. AMB 9(1), 27 (2014)
  * [5] Chen, T., Guestrin, C.: Xgboost: A scalable tree boosting system. In: Proceedings of the 22nd Acm Sigkdd International Conference on Knowledge Discovery and Data Mining, pp. 785–794 (2016)
  * [6] Okada, Y., Okubo, K., Horton, P., Fujibuchi, W.: Exhaustive search method of gene expression modules and its application to human tissue data. IAENG international journal of computer science 34(1), 119126 (2007)
  * [7] Zhang, L., Shah, S.K., Kakadiaris, I.A.: Hierarchical multi-label classification using fully associative ensemble learning. Pattern Recognition 70, 89–103 (2017)
  * [8] Wang, T.: Multi-value rule sets. arXiv preprint arXiv:1710.05257 (2017)
  * [9] Wang, T.: Multi-value rule sets for interpretable classification with feature-efficient representations. In: Proceedings of the 32nd International Conference on Neural Information Processing Systems, pp. 10858–10868 (2018)
  * [10] Garcia, S., Luengo, J., Sáez, J.A., Lopez, V., Herrera, F.: A survey of discretization techniques: Taxonomy and empirical analysis in supervised learning. IEEE Transactions on Knowledge and Data Engineering 25(4), 734–750 (2012)
  * [11] Yang, Y., Webb, G.I.: Discretization for naive-bayes learning: managing discretization bias and variance. Machine learning 74(1), 39–74 (2009)
  * [12] Tou, J.T., Gonzalez, R.C.: Pattern recognition principles (1974)
  * [13] Dodge, Y., Commenges, D.: The Oxford Dictionary of Statistical Terms. Oxford University Press on Demand, ??? (2006)
  * [14] Lowry, R.: Concepts and applications of inferential statistics (2014)
  * [15] Gonzalez, T., Sahni, S., Franta, W.R.: An efficient algorithm for the kolmogorov-smirnov and lilliefors tests. ACM TOMS 3(1), 60–64 (1977)
  * [16] Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al.: Scipy 1.0: fundamental algorithms for scientific computing in python. Nature methods 17(3), 261–272 (2020)
  * [17] Jossinet, J.: Variability of impedivity in normal and pathological breast tissue. Medical and biological engineering and computing 34(5), 346–350 (1996)
  * [18] Horton, P., Nakai, K.: A probabilistic classification system for predicting the cellular localization sites of proteins. In: Ismb, vol. 4, pp. 109–115 (1996)
  * [19] Asuncion, A., Newman, D.: UCI machine learning repository (2007)
  * [20] Ushakov, N., Ushakov, V.: Recovering information lost due to discretization. In: XXXIV. International Seminar on Stability Problems for Stochastic Models, p. 102
  * [21] Chmielewski, M.R., Grzymala-Busse, J.W.: Global discretization of continuous attributes as preprocessing for machine learning. In: Third International Workshop on Rough Sets and Soft Computing, pp. 294–301 (1994)
  * [22] John, G.H., Langley, P.: Estimating continuous distributions in bayesian classifiers. In: Eleventh Conference on Uncertainty in Artificial Intelligence, pp. 338–345. Morgan Kaufmann, San Mateo (1995)
  * [23] Breiman, L.: Random forests. Machine Learning 45(1), 5–32 (2001)
  * [24] Platt, J.: Sequential minimal optimization: A fast algorithm for training support vector machines (1998)
  * [25] Quinlan, J.R.: C4. 5: Programs for Machine Learning. Elsevier, ??? (2014)
  * [26] le Cessie, S., van Houwelingen, J.C.: Ridge estimators in logistic regression. Applied Statistics 41(1), 191–201 (1992)
  * [27] Henriques, R., Madeira, S.C.: Flebic: Learning classifiers from high-dimensional biomedical data using discriminative biclusters with non-constant patterns. Pattern Recognition, 107900 (2021). doi:10.1016/j.patcog.2021.107900