# Big Data Spectra Analysis Using Analytical Programming and Random Decision Forests

Petr Šaloun, Peter Drábik, Ivan Zelinka, Jaroslav Bucko

► To cite this version:

Petr Šaloun, Peter Drábik, Ivan Zelinka, Jaroslav Bucko. Big Data Spectra Analysis Using Analytical Programming and Random Decision Forests. 13th IFIP International Conference on Computer Information Systems and Industrial Management (CISIM), Nov 2014, Ho Chi Minh City, Vietnam. pp.266-277, 10.1007/978-3-662-45237-0_26 . hal-01405596

HAL Id: hal-01405596

https://inria.hal.science/hal-01405596

Submitted on 30 Nov 2016

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Big Data Spectra Analysis using Analytical Programming and Random Decision Forests

Petr Šaloun<fnref n="1" />, Peter Drábik<fnref n="1" />, Ivan Zelinka<fnref n="1" />, Jaroslav Bucko<fnref n="2" />

{petr.saloun, peter.drabik.st, ivan.zelinka}@vsb.cz, bucko.jaroslav@gmail.com

**Abstract.** Spectra analysis on large datasets is in focus of this paper. First of all we discuss a method useful for spectra analysis – analytical programming and its implementation. Our goal is to create mathematical formulas of emission lines from spectra, which are characteristic for Be stars. One issue in performing this task is symbolic regression, which represents the process in our application, when measured data fits the best represented mathematical formula. In past this was only a human domain; nowadays, there are computer methods, which allow us to do it more or less effectively. A novel method in symbolic regression, compared to genetic programming and grammar evolution, is analytic programming. The aim of this work is to verify the efficiency of the parallel approach of this algorithm, using CUDA architecture. Next we will discuss parallel implementation of random decision forest (RDF) to classify huge amounts of various spectra. The mathematical formulas obtained via AP will be used to reduce attributes of explored spectra. Our goal is to propose scalable algorithm for classification of such data, which will preferably need only one pass over data, while maintaining acceptable accuracy. Later we will try to create module compatible with VO and DA data Mining and Exploration project.

**Keywords:** analytical programming, spectra analysis, CUDA, evolutionary algorithm, differential evolution, parallel implementation, symbolic regression, random decision forest, data mining, virtual observatory

## 1 Introduction

Nowadays, astronomy is one of the scientific disciplines which produce enormous amounts of data per day, which needs to be processed automatically. With advancement of surveillance methods the data grows in amount and complexity. In our work we focused on the spectra analysis of Be star-candidates. Be stars are hot B-type stars (effective temperature 10 000 K to 30 000 K) with luminosity class III to V (i.e. not supergiant stars) whose spectrum has shown at least once an emission line – usually hydrogen in the Balmer line, see Figure 1. Sometimes, other emission

<footnote n="1">¹ VŠB-Technical University of Ostrava FEE&I, 17. listopadu 15/2172 CZ-708 33 Ostrava-Poruba, Czech Republic</footnote>

<footnote n="2">² Slovak Technical University Bratislava FIIT, Ilkovičova 2, SK-84216 Bratislava 4 Slovakia</footnote>
lines are visible, for example neutral helium. Even when the spectrum
goes back to normal, the star remains in the Be star class [9]. Some of
them are among the brightest stars in the sky [2]. With the high com-
plexity, the data attributes need to be reduced to less dimensions. The
spectra attributes will be reduced from value for every observed wave-
length in spectra to lines in certain wavelength ranges, which represent
features of object the spectra belongs to.

Non-image data, such as spectra, are also mainly distributed in astron-
omy in Flexible Image Transport System (fits) format, standardized in
1981 [1]. This format can have the *.fts, *.fits, *.fit extensions and
a major feature it has is that images metadata are store in a human
readable ASCII header. In our work we processed fits format data files
with fv FITS Editor, available from the standard Ubuntu repository.

Astronomical data are obtainable from various virtual observatories (VO).
A set of standards to access and manipulate such data was created by
International Virtual Observatory Alliance. They describe API to com-
municate with VO web services. The services provide not only FITS files,
but VOTable formatted data too.

<figure><img src="image_2.png" /><figcaption>Fig. 1. An example of spectra with a characteristic Be star emission.</figcaption></figure>

## 2 Data mining with random decision forests, an introduction

Random decision forests are a data mining structures proposed in 2001 by Leo Breiman [14] for supervised classification. Random decision forests are ensembles of classification and regression trees. They main difference from other decision trees and ensembles is a way of splitting nodes in trees to make decisions.
## 2.1 Training

During training the bagging is used. Each tree in ensemble is grown using subset consisting of cases, drawn random from training set. From this part of training data, a random subset of features is selected. Node split is selected randomly from best splits. The trees are not pruned, unlike in other classification trees. The trees grown during training can be tested for error using out-of-bag estimates.

## 2.2 Prediction

During prediction the sample is classified by every tree in the forest. After that, the final class of sample is class with the most "votes" from trees as shown in Figure 2, or can be selected from mean of probabilities of classes from each tree.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Classification with bagging</figcaption></figure>

## 2.3 Advantages of random forests

Random selection of features brings following advantages to RDF:
- robustness to noise and outliers,

- faster than bagging and boosting,

- simply and easily parallelized,

- no over-fitting.

RDF is considered one of the best learning algorithms, by its performance [16, 17].

## 2.4 Scalable implementation

Parallel implementation of RDF seems to be simple, by growing multiple trees at a time during training stage, and running prediction for multiple trees in forest at prediction stage. An optimal ratio between trees in ensemble and number of computing units should be found to minimize idle time waiting for new spectra samples to be loaded into memory. Since the working data is very large, a single pass over data is very time consuming. Everything should be done in that single pass, or at least in as few passes as possible. This makes decisions for best node split difficult. For this a histogram of values can be used like proposed by Ben-Haim for streaming parallel decision trees in 2010 [18] and used by Li [19]. In spectra analysis it is not necessary to worry about missing values, because missing value in astronomical spectra indicates, presence of certain chemical elements in object, what can be treated as another attribute. When not classifying spectra but other data, the missing values can be handled by substitution to median of values from histogram too. There are currently some scalable solutions for random forests:

**Apache Mahout:** is a data mining library based on Hadoop written in Java. With random forests allows the use of multiple classifiers.

**Scalable random forest by Li B.:** is a random forest implementation with MapReduce. It grows trees breadth-first. Uses histograms to calculate best split [19].

**H2O:** is a scalable machine learning platform, which uses L. Breimans implementation of random forests. It is capable to run on Apache Hadoop and Amazon EC2 (beta). Written in pure java.

## 2.5 Data complexity reduction

To improve speed and accuracy, data complexity reduction and relevant features of spectra selection is needed. The data reduction can be made with the help of analytical programming by identifying interesting features and how they are represented in spectra and describing them by mathematical functions. It will not be necessary to break every observed spectra to set of features, instead, in every tree node, the spectra can be compared against one of the created mathematical functions and based on the result, the decision of which child node descend to can be made.

# 3 Virtual Observatory

The Virtual Observatory is an international astronomical community-based initiative. It aims to allow global electronic access to the available
astronomical data archives of space and ground-based observatories and
other sky survey databases. To obtain data from VO we must follow
IVOA standards for data access layer. DAL services are provided as
HTTP REST web services. For spectra analysis are interesting following
services:

**Simple Spectral Access Protocol** defines a uniform interface to remotely discover and access one-dimensional spectra and aggregations of 1D spectra. Discovered datasets are returned in VOTable format,

**Simple Line Access Protocol** defines interface to spectral lines search
from Spectral Line Data Collections. Spectral lines are returned in
VOTable format.

**Simple Image Access Protocol** allows retrieving on-the-fly created
images of sky given the position and size of the desired output image.
After client has described image, service returns list of possible
images in VOTable format.Client then chooses image to retrieve.
Retrieved image is stored in FITS format.

**Simple Cone Access Protocol** allows to retrieve records from a catalogue of astronomical sources. The query describes sky position and an angular distance, defining a cone on the sky. The response returns a list of astronomical sources from the catalog whose positions lie within the cone, formatted as a VOTable.

We will use these protocols to obtain spectra and metadata of spectra
to process.

4 DData Mining and Exploration

DAME<fnref n="3" /> is an Italian project to provide Astrophysics community with
easy to use data mining suit. This suit is called DAMEWARE. Its focus
is on processing of massive data sets with machine learning methods. It
is based on S.Co.PE.<fnref n="4" />, a general-purpose supercomputing infrastructure
of the University of Naples Federico II. DAMEWARE model library is
extendable via plugins. The only requirement is to have executable of
model and to supply parameters for model or configuration files without
need to know underlying DAMEWARE architecture. We would like to
create a model plugin from our implementation.

5 Overview of analytical programming

Analytical programming, proposed in 2005 [4], was inspired by Hilbert
spaces and genetic programming. The principles and general philosophy
behind analytical programming (AP) stem from these two methods. Into
AP an idea about the evolutionary creation of symbolic solutions is taken
from GP while from Hilbert spaces the idea of functional spaces and the
building of the resulting function by means of the search process are
adopted in AP. The core of AP is based on a set of functions, opera-
tors and so-called terminals, which are usually constants or independent

<footnote n="3">³ http://dame.dsf.unina.it</footnote>

<footnote n="4">⁴ http://www scope.unina.it</footnote>
variables as well as in GP and GE. The main aim of AP is to synthesize
a suitable program which would fit the measured data as well as pos-
sible (with the given precision). For this reason, a discrete set handling
(DSH) idea [7, 8] was adopted in AP [5]. Discrete set handling creates
an interface between the Evolutionary Algorithm (EA) and the problem.
Therefore, we can use in AP almost any evolutionary algorithm. The
individual is represented as a nonnumeric value and a numerical value is
added to the evolutionary process as an integer index. This index repre-
sents an individual from the General Function Set (GFS).

5.1 Versions of AP

There are three versions of AP. AP<sub>basic</sub> is a basic version of AP and uses constants from the terminal set. AP<sub>meta</sub> – in this version there are constants not defined in the terminal set; in the terminal set there is only one general constant K and every constant K is estimated by a different or the same evolutionary algorithm. We can mention one disadvantage of this version – because of running evolution under evolution; it could be very slow for a large number of steps. There are a big number of evaluations of the fitness function. The last version is AP<sub>nf</sub> – constants are estimated by non-linear fitting algorithm.

5.2 General Function Set

The GFS is a set of all mathematical objects – functions, operators and terminals. GFS consists of functions with different numbers of arguments. GFS is user-defined, so the content may differ. We must split the content of GFS into classes based on the numbers of arguments: 0<sub>args</sub> are terminals, 1<sub>args</sub> (sin, cos, tan, ...), 2<sub>args</sub> (+, -, *, /, ...), etc. Choosing the right set may have a large influence on the convergence of AP.

5.3 Evolutionary algorithm

AP was designed to be a very robust method and we can use almost any
evolutionary algorithm. Individual steps, such as mutation, crossover,
etc. are fully handled by the chosen evolutionary algorithm. Operations
which make EAs when the algorithm is run do not have any influence
on performance. Thus the result performance depends mainly on the
correct choice of GFS individuals. For every evolutionary algorithm, the
main goal is generally to reduce the fitness value to below a user-defined
threshold, or to achieve maximum number of migrations for a Self Or-
ganizing Migrating Algorithm (SOMA), or in the case of Differential
Evolution (DE) – generations.

5.4 Mapping operators

Mapping is the phase when an individual is transformed into a useful
mathematical function. It consist of two parts, DSH and security func-
tions, to exclude the creation of pathological individuals. In the evolu-
tionary algorithm the individual is represented by a vector of indexes
and is remapped to mathematical objects from GFS.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Schema of mapping and security principles. Because of measuring distance to end of expression is tan replaced by ω [4]</figcaption></figure>

## 5.5 Reinforced evolution

By running evolution, more or less suitable individuals are created, so one very good idea is to include the best individual as the terminal in GFS. The main idea of reinforcement is based on the addition of a just synthesized and partly successful program in an initial set of terminals [3]. The decision on whether the best value will be added to GFS is ensured by a user defined threshold value. If it is reached, from that moment it is added in GFS as a terminal, best solution and updated whenever a better solution with lower fitness is found.

## 5.6 Security procedures

Evolution can result in an expression, which is not mathematically correct, thus we create a pathological individual (without argument). We can prevent this by distributing GFS into classes ordered by the number of arguments. By mapping from evolution space to GFS, we measure the distance to the end of the expression. When the number of arguments is greater than the distance to the end of the expression, we choose individuals from the lower class. We must also pay attention so as to exclude errors from the fitness function, such as division by zero, functions
with an imaginary or real part (if not expected), frozen functions (an extremely long time to get a cost value), etc. [3].

## 6 Parallel implementation using CUDA – an overview

CUDA (Compute Unified Device Architecture) is a platform for parallel computing, developing by NVIDIA. CUDA is SIMT (Single Instruction Multiple Threads). Its main feature is that one instruction is executed by thousands of threads. However, kernels can effectively perform only basic operations. In the CUDA device various types of memory reside, as you can see in Table 1. Note: registers are the fastest memory on the GPU.

<table><caption>Table 1. CUDA Memory Types and Characteristics [13], (u.c. ≡ unless cached)</caption><thead><tr><th>Memory</th><th>Location</th><th>Cached</th><th>Access</th><th>Scope</th></tr></thead><tbody><tr><td>Register</td><td>On-chip</td><td>No</td><td>Read/Write</td><td>One thread</td></tr><tr><td>Local</td><td>On-chip</td><td>Yes</td><td>Read/Write</td><td>One thread</td></tr><tr><td>Shared</td><td>On-chip</td><td>N/A</td><td>Read/Write</td><td>All threads in a block</td></tr><tr><td>Global</td><td>Off-chip (u.c)</td><td>Yes</td><td>Read/Write</td><td>All threads + host</td></tr><tr><td>Constant</td><td>Off-chip (u.c.)</td><td>Yes</td><td>Read</td><td>All threads + host</td></tr><tr><td>Texture</td><td>Off-chip (u.c.)</td><td>Yes</td><td>Read/Write</td><td>All threads + host</td></tr></tbody></table>

The idea of implementing AP on parallel architecture is based on successful parallel implementation of evolutionary algorithms on CUDA, where a significant speeding up was achieved [10, 11]. We were inspired to implement a parallel version of AP on CUDA, because of the good results achieved by the evolutionary algorithms when implemented on CUDA, although execution of this number of operations and structures may lead to worse performance results, as expected.

## 7 Testing parallel implementation – methods, datasets, results

We cut the spectra from Figure 1 and obtained only the emission line in Figure 5 with 50 equidistant points, to satisfy our requirements. The data on the Be stars spectra come from the archive of the Astronomical Institute of the Academy of Sciences of the Czech Republic<fnref n="5" />.

To be able to carry out a comparison in our tests we used a classical pseudo-random number generator and on the other site we tried to compare it to implemented deterministic chaos, with setting A = 4, which equation is as follows:

$$x_{n-1} = Ax_n(1 - x_n)$$

<footnote n="5">⁵ Available from: http://astropara.projekty.ms.mff.cuni.cz/spectra/newest/</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. The flowchart of AP implementation</figcaption></figure>

We carried out tests with SOMA and DE as the main evolutionary algorithm of AP. For estimating constants we chose the SOMA algorithm. The measured results are shown in Tables 2 and 3.

The settings for our implementation are shown in Table 2.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Extracted emission line (blue) of Be star and our best result (red dots) of the SOMA algorithm with deterministic chaos</figcaption></figure>

<table><caption>**Table 2.** Algorithm settings</caption><thead><tr><th>DE</th><th></th><th>SOMA</th><th></th><th>SOMA constants</th><th></th></tr></thead><tbody><tr><td>NP</td><td>500</td><td>PopSize</td><td>500</td><td>PopSize</td><td>10</td></tr><tr><td>Dimensions</td><td>40</td><td>Dimensions</td><td>20</td><td>Dimensions</td><td>20</td></tr><tr><td>Generations</td><td>100</td><td>Migrations</td><td>10</td><td>Migrations</td><td>5</td></tr><tr><td>F</td><td>0.9</td><td>PRT</td><td>0.1</td><td>PRT</td><td>0.1</td></tr><tr><td>CR</td><td>0.5</td><td>PathLength</td><td>4</td><td>PathLength</td><td>3</td></tr><tr><td></td><td></td><td>Step</td><td>0.21</td><td>Step</td><td>0.25</td></tr></tbody></table>

<table><caption>**Table 3.** Minimum, average and maximum fitness achieved from our tests. Each method was tested 10 times.</caption><thead><tr><th></th><th>DE PRNG</th><th>DE Chaos</th><th>SOMA PRNG</th><th>SOMA Chaos</th></tr></thead><tbody><tr><th>MIN</th><td>4.02327</td><td>4.00432</td><td>0.85327</td><td>0.83326</td></tr><tr><th>AVG</th><td>4.30381</td><td>4.27877</td><td>1.41000</td><td>1.39508</td></tr><tr><th>MAX</th><td>4.57460</td><td>4.51520</td><td>1.84721</td><td>1.91367</td></tr></tbody></table>

<table><caption>**Table 4.** Execution of AP: with SOMA 95238096 evaluations of the cost function in the main EA, with DE 100000000</caption><thead><tr><th>[s]</th><th>DE PRNG</th><th>DE Chaos</th><th>SOMA PRNG</th><th>SOMA Chaos</th></tr></thead><tbody><tr><th>MIN</th><td>2688.59</td><td>1962.53</td><td>3161.47</td><td>2761.57</td></tr><tr><th>AVG</th><td>2805.84</td><td>2296.42</td><td>3305.62</td><td>2942.40</td></tr><tr><th>MAX</th><td>2873.77</td><td>3037.32</td><td>3372.46</td><td>3098.13</td></tr></tbody></table>
Thanks to the use of deterministic chaos against a pseudo-random num-
ber generator (PRNG) we speed up the running time by using simple
mathematical operations, which is good for the CUDA kernel. Conver-
gence of the fitness value is also a bit better.

8 Conclusion

We successfully implemented AP on CUDA and obtained relevant re-
sults. Both PRNG and deterministic chaos are suitable for running with
parallel implementation of AP. The use of deterministic chaos seems to
be a bit better. This is given by the fewer and simpler math operations
needed to get a value. CUDA is not designed to run as long kernels, as
we implemented, but there is a lot of room for future optimization of
these kernels. We use lots of registers in the kernel, so performance falls
from the ideal. In feature CUDA architecture, kernel under kernel will be
probably run, which should bring a large improvement in performance.
We want to compare our results with another parallel approach in the
near future. On graphics hardware it will be Opencl, and OpenPM as a
CPU variant respectively.

**Acknowledgement** The following grants are acknowledged for the financial support provided for this research: Grant Agency of the Czech Republic – GACR P103/13/08195S, is partially supported by Grant of SGS No. SP2014/159, VŠB–Technical University of Ostrava, Czech Republic, by the Development of human resources in research and development of latest soft computing methods and their application in practice project, reg. no. CZ.1.07/2.3.00/20.0072 funded by Operational Programme Education for Competitiveness, co-financed by ESF and state budget of the Czech Republic.

References

1. Wells D.C., Greisen E.W., Harten R.H.: FITS: A Flexible Image Transport System. In Astronomy and Astrophysics Supplement Series. 44:363-370. (1981)

2. Rivinius T., Carciofi A.C., and Martayan C.: Classical Be stars. In The Astronomy and Astrophysics Review 1-86. doi=10.1007/s00159-013-0069-0. (2013)

3. Zelinka I., Davendra D., Senkerik R., Jasek R., and Oplatkova Z.: Analytical Programming – a Novel Approach for Evolutionary Synthesis of Symbolic Structures. Evolutionary Algorithms, Eisuke Kita (Ed.), ISBN: 978-953-307-171-8, InTech, DOI: 10.5772/16166. (2011)

4. Zelinka I., Oplatkova Z., Nolle L.: Analytic Programming – Symbolic Regression by Means of Arbitrary Evolutionary Algorithms. In Special Issue on Intelligent Systems, International Journal of Simulation, Systems, Science and Technology, Volume 6, Issue 9, pp. 44–56, ISSN 1473-8031. (2005)
5. Zelinka I.: Symbolic regression – an overview. Available at http://www.mafy.lut.fi/EcmiNL/older/ecmi35/node70.html

6. Zelinka I., Oplatkova Z.: Analytic programming – Comparative Study. In CIRAS’03: The second International Conference on Computational Intelligence, Robotics and Autonomous Systems. ISSN 0219-6131. (2003)

7. Lampinen J., Zelinka I.: New Ideas in Optimization – Mechanical Engineering Design Optimization by Differential Evolution. Volume 1. London : McGraw-Hill, 20 p. ISBN 007-709506-5. (1999)

8. Zelinka I.: Artificial Intelligence in The Problems of Global Optimization. (in Czech) BEN, 190 p. ISBN 80-7300-069-5. (2002)

9. Thizzy: Be stars. Available at http://www.shelyak.com/contenu.php?id_contenu=30&id_dossier=24 (2008)

10. de Veronese L.P., and Krohling R.A.: Differential evolution algorithm on the GPU with C-CUDA. Evolutionary Computation (CEC), IEEE Congress. (2010)

11. Kralj P. Differential Evolution with parallelised objective functions using CUDA. Available at http://www.codehunter.eu/media/Kralj_Differential_Evolution_with_parallelised_objective_functions_using_CUDA.pdf (2013)

12. Kromer P., Platos J., Snasel V., and Abraham A.: Many-threaded implementation of differential evolution for the CUDA platform. In Proceedings of the 13th annual conference on Genetic and evolutionary computation (GECCO ’11), Natalio Krasnogor (Ed.). ACM, New York, NY, USA, 1595-1602. DOI=10.1145/2001576.2001791 http://doi.acm.org/10.1145/2001576.2001791. (2011)

13. Farber R.: CUDA Application Design and Development. (1st ed.). Morgan Kaufmann Publishers Inc., San Francisco, CA, USA. (2011)

14. Breiman L.: Random Forests. Machine Learning, 45:532. (2001)

15. Breiman L.: Bagging predictors. Machine Learning, 24, 2:123-140. DOI=10.1023/A:1018054314350 http://dx.doi.org/10.1023/A:1018054314350. (1996)

16. Caruana R., Niculescu-Mizil A.: An Empirical Comparison of Supervised Learning Algorithms. In Proceedings of the 23rd International Conference on Machine Learning. ACM, New York, NY, USA, 161-168. DOI=10.1145/1143844.1143865 http://doi.acm.org/10.1145/1143844.1143865. (2006)

17. Ho T. K.: The Random Subspace Method for Constructing Decision Forests. IEEE Trans. Pattern Anal. Mach. Intell. 20, 8, 832-844. DOI=10.1109/34.709601 http://dx.doi.org/10.1109/34.709601. (1998)

18. Ben-Haim Y., Tom-Tov E.: A Streaming Parallel Decision Tree Algorithm. J. Mach. Learn. Res. 11, 849-872. (2010)

19. Li B., Chen X., Li M. J., Huang J. Z., Feng S.: Scalable random forests for massive data. In Proceedings of the 16th Pacific-Asia conference on Advances in Knowledge Discovery and Data Mining - Volume Part I (PAKDD’12), Pang-Ning Tan, Sanjay Chawla, Chin Kuan Ho, and James Bailey (Eds.), Vol. Part I. Springer-Verlag, Berlin, Heidelberg, 135-146. DOI=10.1007/978-3-642-30217-6_12 http://dx.doi.org/10.1007/978-3-642-30217-6_12. (2012)