# JAMIP: an artificial-intelligence aided data-driven infrastructure for
computational materials informatics

Xin-Gang Zhao  These authors contributed equally  State Key Laboratory of
Integrated Optoelectronics, College of Materials Science and Engineering,
Jilin University, Changchun 130012, China  Kun Zhou  These authors contributed
equally  State Key Laboratory of Integrated Optoelectronics, College of
Materials Science and Engineering, Jilin University, Changchun 130012, China
Bangyu Xing  These authors contributed equally  State Key Laboratory of
Integrated Optoelectronics, College of Materials Science and Engineering,
Jilin University, Changchun 130012, China  Ruoting Zhao  These authors
contributed equally  State Key Laboratory of Integrated Optoelectronics,
College of Materials Science and Engineering, Jilin University, Changchun
130012, China  Shulin Luo  State Key Laboratory of Integrated Optoelectronics,
College of Materials Science and Engineering, Jilin University, Changchun
130012, China  Tianshu Li  State Key Laboratory of Integrated Optoelectronics,
College of Materials Science and Engineering, Jilin University, Changchun
130012, China  Yuanhui Sun  State Key Laboratory of Integrated
Optoelectronics, College of Materials Science and Engineering, Jilin
University, Changchun 130012, China  Guangren Na  State Key Laboratory of
Integrated Optoelectronics, College of Materials Science and Engineering,
Jilin University, Changchun 130012, China  Jiahao Xie  State Key Laboratory of
Integrated Optoelectronics, College of Materials Science and Engineering,
Jilin University, Changchun 130012, China  Xiaoyu Yang  State Key Laboratory
of Integrated Optoelectronics, College of Materials Science and Engineering,
Jilin University, Changchun 130012, China  Xinjiang Wang  State Key Laboratory
of Integrated Optoelectronics, College of Materials Science and Engineering,
Jilin University, Changchun 130012, China  Xiaoyu Wang  State Key Laboratory
of Integrated Optoelectronics, College of Materials Science and Engineering,
Jilin University, Changchun 130012, China  Xin He  State Key Laboratory of
Integrated Optoelectronics, College of Materials Science and Engineering,
Jilin University, Changchun 130012, China  Jian Lv  State Key Laboratory of
Superhard Materials, College of Physics, Jilin University, Changchun 130012,
China  Yuhao Fu<EMAIL_ADDRESS>State Key Laboratory of Superhard
Materials, College of Physics, Jilin University, Changchun 130012, China
Lijun Zhang<EMAIL_ADDRESS>State Key Laboratory of Integrated
Optoelectronics, College of Materials Science and Engineering, Jilin
University, Changchun 130012, China

###### Abstract

Materials informatics has emerged as a promisingly new paradigm for
accelerating materials discovery and design. It exploits the intelligent power
of machine learning methods in massive materials data from experiments or
simulations to seek new materials, functionality, and principles, etc.
Developing specialized facilities to generate, collect, manage, learn, and
mine large-scale materials data is crucial to materials informatics. We herein
developed an artificial-intelligence-aided data-driven infrastructure named
Jilin Artificial-intelligence aided Materials-design Integrated Package
(JAMIP), which is an open-source Python framework to meet the research
requirements of computational materials informatics. It is integrated by
materials production factory, high-throughput first-principles calculations
engine, automatic tasks submission and monitoring progress, data extraction,
management and storage system, and artificial intelligence machine learning
based data mining functions. We have integrated specific features such as an
inorganic crystal structure prototype database to facilitate high-throughput
calculations and essential modules associated with machine learning studies of
functional materials. We demonstrated how our developed code is useful in
exploring materials informatics of optoelectronic semiconductors by taking
halide perovskites as typical case. By obeying the principles of automation,
extensibility, reliability, and intelligence, the JAMIP code is a promisingly
powerful tool contributing to the fast-growing field of computational
materials informatics.

##  I Introduction

Currently, materials science research is entering a new paradigm featuring
materials informatics, which applies machine-learning techniques to massive
materials dataAgrawal and Choudhary (2016); Schleder  _et al._ (2019); Draxl
and Scheffler ; Himanen  _et al._ to accelerate materials discovery and
design. It was first promoted in part by the Materials Genome InitiativeKalil
and Wadia (2011) and emerged in practice as the result of the fast development
of experimental materials synthesis, characterization approaches, and
theoretical materials simulation through available computation resources that
yield massive materials dataAgrawal and Choudhary (2016). By utilizing the
recognized predictive power of materials informatics, encouraging
breakthroughs have been made, including the design of new materials, the
identification of new relationshipsGhiringhelli  _et al._ (2015); Duvenaud
_et al._ (2015), and the prediction of new principlesSchmidt and Lipson
(2009); Rudy  _et al._ (2017). In the fields of both experimental and
computational materials informatics, there is an urgent need to develop
powerful facilities/infrastructures to meet research requirements under this
new paradigm.

Generating massive materials data of the composition-structure-property
relationship is essentially the first step for materials informatics studies.
In the field of computational materials science, there are generally three
materials data generation techniques: The first is through direct calculations
of single or a few materials. This is used frequently in the old paradigm when
theorists want to study some material properties observed experimentally or to
explore new physics in some previously unstudied materials. Such data are
dispersedly distributed in literature and may be gathered manually or by
automatic text extraction approachesFundel  _et al._ (2007); Tshitoyan  _et
al._ (2019). The second is through the high-throughput (HT) calculations where
extremely large numbers of materials spanning different chemical compositions
are automatically calculated based on prototype structures. Along this
direction, a series of computational infrastructures are developed and used in
high-throughput calculationsCurtarolo  _et al._ (2012a); Zhao  _et al._
(2017); Luo  _et al._ (2021); Curtarolo  _et al._ (2012b); Pizzi  _et al._
(2016); Mathew  _et al._ (2016); Jain  _et al._ (2015); Yang  _et al._ (2018);
Zhu  _et al._ (2020); Mathew  _et al._ (2017); Choudhary  _et al._ (2020);
Wang  _et al._ (2021). The accumulated computational materials data lead to
the formation of the databases that are accessible for query and research
purposesPickard and Needs (2011); Hachmann  _et al._ (2011); Kirklin  _et al._
(2015); Draxl and Scheffler (2019); Gorai  _et al._ (2016). The third is
through crystal structure searches. These seek new stable and metastable
materials structures under the typical chemical composition. This method
complements high-throughput calculations. It explores the potential energy
profile of materials, usually by combining artificial intelligence global
optimization methods with energetic calculationsGlass  _et al._ (2006); Lonie
and Zurek (2011); Pickard and Needs (2011); Wang  _et al._ (2012); Tipton and
Hennig (2013). Because in structure search studies, the object of greatest
concern is the lowest-energy ground-state structure, the generated materials
data are usually regarded as intermediate products, and do not accumulate in
materials informatics studies.

Machine learning techniques act as an engine that triggers/accelerates
materials discovery in computational materials informaticsSchmidt  _et al._
(2019); Liu  _et al._ (2017); Ramprasad  _et al._ (2017). Various machine-
learning models are already available for use in computational material
informatics, such as supervised learning, unsupervised learning, and active
learningKusne  _et al._ (2020). Supervised learning requires an artificially
constructed training process, which then understands the mapping relationships
between material descriptors and material properties. Using this model, both
the underlying physical laws between descriptors and the material properties
can be minedOuyang  _et al._ (2018); Bartel  _et al._ (2019). Additionally, it
can be used to reverse-engineer new materials to uncover new materialsLu  _et
al._ (2018). Unsupervised learning, by contrast, is generally used to uncover
differences among materials or chemical systems in terms of unmarked
“descriptors”. Clustering algorithms are the most widely used for this
purpose. For example, Chen et al. used the K-mean clustering algorithm to
analyze the oxygen diffusion pattern, hopping statistics, and site occupation
within crystalsChen  _et al._ (2015). In addition, in the field of materials
discovery, active learning and Bayesian optimization, with the self-optimizing
characteristics of model, allow us to search the potential materials space for
high-quality materials in a limited data range, and contribute to the material
experimental and computational designTran and Ulissi (2018).

The software/infrastructure for computational materials informatics needs to
meet at least the following requirement, but is not limited to them: (1)
Adaptability for diverse materials discovery and design: The software should
flexibly deal with different material systems from inorganic and organic to
organic-inorganic hybrid systems, and different materials
variations/combinations, such as alloy or defect, surface or interface, and
heterostructure or superlattice. To fit the potential complexity of different
materials with different properties, computational workflow modules should
have high scalability and ensure flexible combinations among different types
of computing tasks. (2) Efficiency and reliability in data generation and
management: A high-level automation framework to initialize, run, and analyze
large-scale high-throughput calculations is essential. Monitoring
computational tasks in real-time and reporting/correcting potential errors are
important. The software should have tools to efficiently extract and collect
results and store them in a self-contained database. (3) Orientation by
materials functionality: The software should be oriented by the functionality
of materials (e.g., photovoltaic, thermoelectric, ferroelectric, and
catalytic), which is central to materials discovery. Rational design of
workflow modules needs to be sufficiently considered to effectively calculate
the functionalities of different material systems. (4) Synthetical integration
of data generation, storage, processing, and learning: It is desirable to
integrate data generation, storage, processing, and learning modules with a
unified infrastructure framework, which greatly facilitates the data fluxion
and conversion procedure, as well as the efficiency of data mining. In the
integrated framework, newly calculated data motivated by feedback from data-
learning procedures may accelerate knowledge accumulation.

In this paper, we report on an artificial-intelligence-aided data-driven
infrastructure that we have been developing to fulfill the above requirements
for computational materials informatics. It is an open-source Python-
integrated framework named the Jilin Artificial-intelligence-aided Materials-
design Integrated Package (JAMIP). The code is integrated by intimately
connected units of Data generation (e.g., high-throughput materials
calculations and screening), Data collection (e.g., automatic data extraction,
management, and storage) and Data learning (e.g., integrated descriptors
engineering, and machine learning based data mining functions). Below, we
describe the JAMIP code framework and its usage in carrying out materials-
informatics-related research on optoelectronic semiconductors by taking halide
perovskites as instances. The code and manuals describing its detailed usage
options are freely available for download for academic research.

##  II Code Framework

The organization of JAMIP abides by the data lifecycle in the research field
of computational materials informatics, from data generation to collection and
learning, as shown in Fig. 1. In the data generation stage, users perform the
large-scale high-throughput materials calculations, which are done through
obtaining structure prototypes from the JAMIP database, generating candidate
structures through the materials production factory, initializing the input
parameters of HT calculation tasks, submitting and monitoring tasks, and post-
processing results. The code contains a specially designed computational tasks
allocation and monitoring system to guarantee efficient and reliable large-
scale calculations. We have designed packaged material property calculation
modules ranging from geometrical, thermodynamic, mechanical, electronic,
phononic, photonic to magnetic properties. There are modules for material
property optimization calculations guided by artificial intelligence
optimization methods (e.g., genetic algorithm and particle swarm optimization
algorithm) to achieve optimal functionality of materials with minimal effort
and global sampling in the materials space. Users can conveniently perform
functionality-direct materials screening at this stage. In the data collection
stage, the code has continuously developed automatic toolkits for effective
extraction, collection, and analysis on target data of functional materials.
There is a self-contained professional database for automatically normalizing
and storing the “composition-structure-property” data of materials. The
database is developed based on the Django framework that facilitates access to
data without SQL query language. In the final data learning stage, the uses
can exploit machine-learning methods to perform data mining researches. Our
data learning module implements the popular general machine learning
workflows, such as properties prediction, descriptors search, and materials
classification. The learned principle knowledge from the data learning stage
can in turn motivate and guide new calculations in the first data generation
stage. The developed machine learning potentials can meanwhile accelerate
high-throughput materials calculations.

Figure 1: Overview of the JAMIP code framework. The program comprises three
major parts based on the material data’s lifecycle: Data generation(blue),
Data collection(yellow), and Data learning(green).

2.1 Data generation

Materials data generation involves the configuration and initialization of the
computation environment and the customization, execution, and monitoring of
the calculation tasks (Fig. 2).

Figure 2: Detailed workflow of the HT calculations.

2.1.1 Materials production factory

Currently, complex material structure manipulation mainly relies on third-
party graphical interactive software. Although such software facilitates the
users’ intuitive operation of the material structure, it also adds some
obstacles when performing structural transformation in large quantities.
Another tricky problem is the lack of available interfaces between graphical
interactive software and HT software. With the increasing complexity of the
scientific issues studied, the demand for complex structural modeling
capabilities is increasing for HT calculations. Thus, powerful structural
modeling capabilities are an important feature of next-generation HT software.
These include building complex twin structures and structures with defects
(e.g., plane or screw dislocations), and even generating new structures by
using some structural units or symmetry operation methods based on group
theory.

JAMIP has two kinds of material structure manipulation operators. One is a set
of basic material structure manipulation operators (e.g., deleting atoms,
adding atoms, scaling lattice parameters, supercells, defining crystal
structure, and adding vacuum layer); another is a set of high-level complex
material structure manipulation operators. They are from a flexible
combination of the above basic manipulation operators (e.g., substituting an
atom with a given molecule and constructing a crystal structure with twin
interfaces). Users can easily further expand and customize the material
structure manipulation operators based on the built-in manipulation operators.

The other important feature of our materials production factory is the
integration of our developed inorganic crystal structure prototype database
(ICSPD) as a sub-database of JAMIP’s database. The structure prototype
database can be beneficial for high-throughput calculations where usually
tremendous materials spanning different chemical compositions are
automatically calculated based on the prototype structuresLuo  _et al._
(2021). Our ICSPD was developed in terms of the local atomic environments
(LAEs) of all experimentally known inorganic crystal structures, with the help
of the unsupervised learning algorithm, i.e., hierarchical clustering. It is
distinct from the common structure prototypes database that was constructed
directly based on the space group information of crystal structures. It
provides a more reasonable and flexible classification based on microscopic
atomic crystal structure features. This self-contained structure prototype
database effectively facilitates ongoing high-throughput calculations in
JAMIP.

2.1.2 High-throughput calculation workflow

HT calculations involve the transfer of a large number of different types of
data. This requires the workflow module to correctly transfer data after
combining different calculation tasks. JAMIP includes a task pool to manage
all calculations tasks in a unified manner and transfer the data with a
function called “calculation”. By using the observer patterns, the task can
monitor the status of tasks on the computing nodes and submit new calculation
tasks in time to avoid lags in information updates and long periods of
idleness. Furthermore, some types of calculation tasks, such as phonons, can
generate multiple parallel calculations concurrently. JAMIP can even split
these calculations into multiple compute nodes to speed up the calculations.

Scalability is one of the essential features in the calculation workflow
module because of complex and diversified computing task requirements. We drew
on the factory method in the software design pattern. We designed many
calculation templates of crucial properties, especially the properties related
to optoelectronic materials. In HT calculations, these property calculations
can be arbitrarily combined by adding their name to the tasks array in the
task control file named “input.py”. Moreover, we designed a self-defined
template to meet the needs of the customized calculation tasks. Users only
need to write a small amount of code and can implement their customized
calculations. To implement one customized calculation, users need to add the
name of this task in the “input.py” file and prepare a workflow file with this
task’s name. The tasks module can automatically detect this customized
calculation located in the “diyflow” directory and perform the calculations by
calling its diy_calculator() method. The HT calculations involve the
configuration of computing resources. In JAMIP, all configuration information
is stored in the YAML files in the “env” dictionary. We have built-in
configuration templates for different job management systems and compiler
environments. Currently, JAMIP supports three job management systems: PBS,
LSF, and SLURM. The code provides currently complete interface support for the
Vienna Ab initio Simulation Package (VASP)Kresse and Furthmüller (1996),
partial interface support for Quantum EspressoGiannozzi  _et al._ (2009), and
for more widely used ab initio codes in future releases.

2.1.3 Material property calculations

JAMIP contains a large number of built-in computational templates related to
the semiconductor’s properties, including thermodynamic, geometrical,
electronic, mechanical, phononic, and photonic properties (listed in Table
II). Some of these calculation tasks interact with the other modules of the
program, such as the JAMIP’s database and materials production factory.
Considering the calculation of the decomposition pathway as an example, the
program can automatically search for crystallographic information of competing
phases related to this material in the JAMIP database, calculate all possible
decomposition paths in batches, and output visual analysis results. Most
importantly, these calculations can be combined arbitrarily into one HT
calculation task, in which it can simultaneously calculate multiple properties
and handle their calculated results. Furthermore, as discussed above, JAMIP
can easily be extended to new customized calculation tasks. Note that the
above calculations need to be combined with VASP. For Quantum Espresso, the
current program only provides some basic calculations, including structural
optimization, self-consistent field (SCF) calculation, and electronic band
structure calculation. We are continuously expanding the items for property
calculations.

Table 1: Built-in property calculations. Properties | Specific terms
---|---
Thermodynamics | Total energy; Enthalpy; Convex hull; Triangle zone; Decomposition of the pathway
Geometry | Atomic packing factor; Index of structural connectivity; Tolerance factor;
| X-ray diffraction (XRD) pattern; Scanning tunneling microscope (STM) image;
| Radial distribution function (RDF)
Electronic structure | Band; DOS/PDOS; Electrostatic potential; Effective mass; Born effective charge;
| Deformation potential; Dielectric constant; Electron localization function
(ELF);
| Charge density; Crystal orbital Hamilton population (COHP/COOP); Bond
valence sum
Mechanics | Stress tensor; Elastic constants; Bulk modulus; Poisson,s ratio; Hardness
Phonon | Phonon dispersion; PHDOS; Raman $\&$ IR; Vibration entropy; Zero-point energy
| (ZPE); Thermal coefficient; Gruneisen parameters; Heat capacity;
| Thermal conductivity; Soft mode transition
Photon | Optical bandgap; Absorption spectrum; Jointed density of states (JDOS); Dipole
| matrix elements (DME); Spectroscopic limited maximum efficiency (SLME);
| Exciton binding energy
Magnetism | Magnetic moment, Magnetism type
Ferroelectrics | Electronic polarization magnitude, Polarization pattern

2.1.4 Calculations guided by artificial intelligence optimization methods

In certain material systems with complex degrees of freedom in composition or
structure variation, it is challenging for the HT calculations to exhaust the
materials to be explored and exhibit a combinatorial explosion in numbers. To
achieve optimal functionality of materials with minimal effort and global
sampling in the materials space, we use artificial intelligence optimization
methods to guide the material property calculationsFranceschetti and Zunger
(1999); Zhang  _et al._ (2013, 2012). This approach in nature is the same as
in crystal structure search studies that attempted to find the lowest-energy
material structures through the intelligent exploration of the potential
energy landscapeGlass  _et al._ (2006); Lonie and Zurek (2011); Pickard and
Needs (2011); Wang  _et al._ (2012); Tipton and Hennig (2013). Moreover,
fitness is the material property rather than the free energy in crystal
structure search studies. The artificial intelligence optimization methods
implemented include genetic algorithm and particle swarm optimization
algorithm. This allows the users to optimize the target material property in
the materials configuration space that needs to be prescribed. Such studies
are meaningful in exploring the extreme property values that the studied
materials can reach, or in studying the material systems with masses of
metastable states that can be experimentally realized. Otherwise, the
stability of the materials needs to be properly considered by taking the free
energy as an equally important criterion.

2.1.5 Real-time monitoring and calculation-error management

Adaptive ability relating to the diversified calculation tasks is one of the
important indicators for evaluating the quality of HT computational software.
The default calculation parameters may not be suitable for every sub-
calculation task in large-scale calculations owing to the differences in the
structure and composition. Testing and modifying is a very cumbersome job in
the parameter debugging of large-scale calculations. It requires a lot of
manual intervention and even leads to the failure of tasks in some serious
cases. In JAMIP, we append the monitoring and self-correction mechanism to
deal with run-time errors in the calculations. The back-end monitor process is
triggered after finishing a calculation task. This process judges whether or
not the calculation is completed normally. For abnormally stopped
calculations, the program’s control flow automatically switches to the
correction step of the error. In the current version, we collected some common
solutions and constructed an “error-solution” database. The sketched flow
chart is shown in the monitor module of Fig. 3. By retrieving this database,
JAMIP can automatically try to correct errors encountered in large-scale
calculations. More importantly, the “error-solution” database can be expanded
by adding customized solutions of the new error to the YAML-typed file named
“error.yaml”, which comprises three parts: name, keywords, and solution for
the error, as shown in Table I. With the enrichment of error solutions, the
ability of error correction will become more powerful and intelligent and
significantly reduce the intensity of human intervention.

Table 2: Summary of control files. Calculation params configuration | Running environment configuration
---|---
base: |  |  manager:PBS |
ismear:0 |  |  queue:defaut |
npar:4 |  |  cores:16 |
sigma:0.05 |  |  nodes:1 |
algo:fast |  |  maximum:10 |
band: |  |  env: module load python3 |
lcharg:11 |  | … |
…: |  |  |
Errors and Corrections configuration
ISYM: |  |  |
\- VERY BAD NEWS! internal error in subroutine IBZKPT |  |  |
\- isym: 0 |  |  |
… |  |  |
Job submission script
1 Input: Structures, calculation params, … |  |  |
2 Output: Task Pool |  |  |
3 Function jump_input(): |  |  |
4 vasp = SetVasp() |  |  |
5 vasp.program =‘/executable/vasp/program’ |  |  |
6 vasp.potential=‘/directory/of/vasp/pseudopotentials/folders’ |  |  |
7 vasp.task = ‘scf band partchg emass’ |  |  |
8 vasp.force = 1e-3 | #force convergence |  |
9 vasp.energy = 1e-6 | #energy convergence |  |
10 vasp.cutoff = 1.3 | #cutoff energy of ware function |  |
11 vasp.kpoints = 0.25 | #k-mesh density |  |
12 pool=Prepare.pool(vasp) |  |  |
13 pool.set_structure(‘input’,‘output’) | #read structures |  |
14 pool.save(‘silicon.dat’,overwrite=True) | #save pool of tasks |  |
15 Prepare.cluster(‘pbs’) | #check environment configuration |  |
16 Prepare.incar(‘vasp.tasks’) | #check calculation configuration |  |

2.2 Data collection

JAMIP materials database is based on the Django framework that was developed
for database-driven Python web applications. Django is an excellent middleware
that link the database and application-layered code. It can avoid users
directly accessing the database by using obscure standard query language (SQL)
and reduce the difficulty of development. All development activities are in a
pure Python environment. JAMIP can easily switch to other relational
databases. Currently, JAMIP supports MySQL, MariaDB, PostgreSQL, Oracle, and
SQLite. As shown in Fig. 3, the database engine can read crystallographic
information from different types of raw structure data files. Our developed
material production factory can be used as another important source of
structures. In addition, the database engine can batch-extract calculation
results in the data generation stage by executing the command “jp –db” in the
terminal. All the material data are decomposed into some basic data units by
the database engine and stored in the underlying relational database.

Figure 3: The overall perspective of the material database.

2.2.1 Data extraction

Once the HT calculations are completed, users need to extract information from
different calculation files, summarize/analyze the calculation data, and
intuitively obtain the performance differences between materials. For
instance, the results analysis of band-alignment calculations involves the
following: (1) finishing HT calculations of the heterojunction, including
structural optimization, self-consistent field and band structure calculation;
(2) extracting the vacuum energy level from the electrostatic potential
calculated by self-consistent field and the electronic structure from the
output files of band structure calculations, and searching the VBM and CBM of
materials by the electron occupation; (3) plotting the band offsets diagram
with the aligned band sides, and arranging it in VBM order from low to high.
We packaged the above data processing flow as a separate post-processor. Using
this module, users can complete the process of HT calculations, data
extraction, and intuitively obtain the relative band edge position and
heterojunction type of materials from the band offset diagram.

JAMIP has a complete built-in data extraction module. Users can extract the
required information from HT output files with the post-processing program and
store it as images, or into a comma-separated values (CSV) file and the
database. The post-processing program currently supports the following
functions: (1) Total energy, enthalpy, decomposition of the pathway, convex
hull, and triangle zone; (2) Band structure, electron density of states,
charge density, effective mass, born effective charge, deformation potential,
dielectric constant, electron localization function, and bond valence sum; (3)
Stress tensor, elastic constants, bulk modulus, and Poisson’s ratio; (4)
Phonon dispersion, Raman $\&$ IR, vibration entropy, zero-point energy (ZPE),
thermal expansion coefficient, thermal conductivity, and soft mode transition;
(5) Optical bandgap, absorption spectrum, jointed density of states (JDOS),
dipole matrix elements (DME), spectroscopic limited maximum efficiency (SLME),
and exciton binding energy; etc.

2.2.2 Data formatting and normalization

In JAMIP, the calculated raw materials data are extracted and decomposed into
a series of Python objects, including entries, structures, compositions, and
elements. These objects are designed with deep optimizations based on the
Django framework. They are mainly composed of attributes, associative
relationships with other objects, and their auxiliary methods. Each object
maps to a single database table in our relational database. Among these
objects, entries are used to save calculation-related information, such as
computational parameters and calculated results (e.g., bandgaps and
electron/hole effective masses). Moreover, each entry object also records its
corresponding associated structure. Similarly, other objects contain a large
number of attributes and association relationships. For example, an element
object stores a large amount of element-related information (e.g., element
symbol, atomic number, and atomic radius) and associated objects such as
structures, compositions, and atoms. The establishment of these relationships
forms a complex relationship network, which is the key to material
informatics. The formation of these relationship networks also allows users to
easily query and filter out the required data by using some simple Django
query syntax.

2.2.3 Data storage in the self-contained database

The SQL query is time-consuming because of data filtering and merging
operations that are required for large amounts of data entries. Optimizing the
efficiency of the SQL query is the key to improving the database’s
performance. For example, material structure manipulation involves a large
number of queries and operations on the structural information (e.g.,
extracting atoms within a given region and removing these atoms). Frequently
performing the SQL query operations significantly slows down the response
speed of the program. Thus, we designed an additional run-time data structure
in memory for material structure manipulation. Inside a series of Python
objects such as structures, compositions, and elements, we maintain a series
of attribute arrays synchronously. These attribute arrays store all objects
associated with them separately. When adding or deleting some atoms, all
attribute arrays can be refreshed by calling the update() method in the file
“structure.py”. By directly accessing these data in memory, the time spent
querying the same data is reduced by two orders of magnitude compared to
querying the data via the SQL query language. However, it should be noted that
the query syntax in the run-time data structure is different from the one in
the database. The key is that only structure objects are allowed to access the
run-time data structure. For instance, Table III shows the query syntax and
efficiency when extracting all atoms of the specified element (Na) in a
structure object (s0).

Table 3: Query efficiency comparison between the two query methods. Type |  Query statement |  Time-consuming ($\times$10-5s))
---|---|---
Run-time data | s0.get_element(‘Na’).atoms |  3.4
structure database |  Element.objects.get(symbol=‘Na’).atom_set.filter(atom_strucutre=s0) |  290
Figure 4: Detailed workflow in the Machine-Learning.

2.3 Data learning

After completing the HT computing tasks, materials data mining can be
performed by studying and learning the data stored in the JAMIP database. This
version of the JAMIP’s data mining module is constructed under the supervised
learning model. It includes data preprocessing, material feature engineering,
machine learning model construction, and performance evaluation (Fig. 4).

2.3.1 Data filtering and preprocessing

Based on the JAMIP structure classes and database, users can retrieve
information on chemical composition (e.g., ionic radius, electronegativity,
and ionic potential), crystal structure (e.g., packing factor, radial
distribution, and smooth overlap of atomic position), diagrammatic function
(e.g., XRD diagram), electronic structure (e.g., band structure and density of
state), and optoelectronic properties from structures that have completed the
HT tasks. Before inputting the data in the machine learning model, the data
need to be converted into a language that can be identified by the model for
sufficient completion. However, there are often missing values or outliers in
the building material datasets. For example, different structures may have
inconsistent dimensions when converted to structural features, or the data
distribution may have a “long tail effect”, JAMIP includes a data cleaning
function to deal with missing values and outliers. In addition, the original
dataset generally has category features, such as crystal systems and space
groups. It also includes a category feature encoding (e.g., one-hot and
binary) function to ensure that the data can be identified.

2.3.2 Material feature engineering

The descriptors in the original dataset are mostly constructed from intuition.
However, most descriptor sets appear highly dimensional or lack features,
which often affects the learning of material properties. Thus, a robust method
for identifying descriptors is needed. JAMIP’s material feature engineering is
divided into the following steps: 1) Feature scaling: This is a normalized
data method to ensure that the role of each feature is close to the material
target properties (e.g., normalization, non-linear normalization, and
regularization); 2) Feature construction: As the relationship between
descriptors and target properties is not simple and is often an exponential or
logarithmic relationship, users can customize features by setting the
“expression_of_feature” parameter in the “feature_crosses” method; 3) Feature
selection: After the feature construction is completed and the number of
features is added, it also imposes redundant features. This increases the
model complexity and enhances the overfitting possibility. JAMIP offers
multiple methods to accomplish feature dimensionality reduction to help
uncover the features most closely associated with the material target
properties and identify the physical mechanisms behind them, including chi-
squared feature selection, tree-based model feature selection, and recursive
feature selection.

2.3.3 Machine learning model and performance evaluation

We bulid the machine-learning models by using the scikit-learn machine
learning packageFranceschetti and Zunger (1999). Currently, JAMIP provides
four types of machine learning algorithms: regularization linear models,
clustering models, kernel models, and tree-based models, and provides post-
fitting functions for the above algorithms. For different learning tasks
(classification or regression), the user can select different evaluation
metrics (e.g., root mean square error and confusion matrix) to test the fitted
model. For randomness in the integrated model, the program takes multiple
fittings and averages them to eliminate as much randomness as possible.

##  III Demonstration of usage

In principle, our JAMIP code can be used to perform materials-informatics-
related research on the functional material system of which the functionality-
involved properties can be accurately described by first-principles
calculations. In this section, we briefly describe its use in studying halide
perovskites-based semiconductors for optoelectronic applications. Halide
perovskites (HPs, formula as ABX3), as promising optoelectronic materials
(e.g., photovoltaic solar cells, light-emitting diodes, lasers, and
photodetectors), have attracted tremendous research interest. We chose this
family of optoelectronic semiconductors as instances.

Figure 5: Summary of halide perovskites and derivatives in HT calculations, in
which ABX3 and A2B1B2X6 (FA=NH2CHNH${{}_{2}}{{}^{+}}$;
MA=CH3NH${{}_{3}}{{}^{+}}$; HA=NH3OH+; DA=NH2NH${{}_{3}}{{}^{+}}$;
EA=CH3CH2NH${{}_{3}}{{}^{+}}$; FM=CHONH${{}_{3}}{{}^{+}}$;
M=NH${{}_{4}}{{}^{+}}$; DEA=CH3NH2CH${{}_{3}}{{}^{+}}$) for three-dimensional
materials, A2BX6 (A=Cs; B=Pd, Sn, Ti; X= Cl, Br, I) and A3C2X9 (A=Cs; C=Bi,
In, Sb; X=Cl, Br, I) for 0-dimentional materials, ABX3 (A=Cs, B=Pb, Sn; X=Cl,
Br, I), A2BCX6 (A=Cs; B=Ag; C=Bi, In; X=Cl, Br, I), and A3C2X8 (A=Cs; C=Bi,
In; X=Cl, Br, I) for one-dimensional materials; A3C2X9 (A=Cs; C=Bi, In; X=Cl,
Br, I) and A4BC2X12 (A=Cs; B=Cu; C=Bi, In; X=Cl, Br, I) for two-dimensional
materials.

Fig. 5 illustrates the studied materials space of the candidate structures of
halide perovskites. Besides extensively studied single and double perovskites,
we also included some important low-dimensional perovskite-like structures.
Nearly 700 structures were explored by employing the JAMIP code (Fig. 6a). We
constructed inputs of three-dimensional candidate structures based on the
crystal structure prototypes of single, double, and low-dimensional
perovskites with given chemical formula. Note that the lattice parameters can
be automatically adjusted according to the composition elements when
generating the structures. The orientation of organic molecules in the
candidate structures is generated by using the structure factory method to
rotate the molecular orientation appropriately. The way to generate the two-
dimensional perovskite structure is as follows: (i) we first redefined the
crystal lattice of single perovskite structure, rotated the halogen atoms
layer by layer to form a twisted octahedron, and (ii) then deleted the outward
atoms in the cell and established two-dimensional perovskites.

Next step, by using generated structures as inputs, we carried out the HT
calculations of the electronic and optical properties of these halide
perovskites with a series of different compositions and structural
dimensionality. For each candidate structure, we set the automatic calculation
flow to work as following: (i) firstly perform structural optimization on the
perovskite structure and self-consistent field calculation to obtain stable
crystal structure and ground-state charge density; (ii) calculate the
properties of perovskites, including band structure, effective mass,
absorption spectrum, dielectric constant, exciton binding energy, and dipole
matrix element, etc. To facilitate analysis and obtain knowledge, we extracted
structure-, energetic-, electronic-, and optical-related information from the
HT calculation output files, and imported these information into the JAMIP
database. Some of the calculated results were shown in Fig. 6 as the function
of the packing factor (PF) that is an important descriptor to judge the
structure stability of perovskites. Fig. 6b shows non-monotonic trends of
calculated electronic properties in 3D and 2D perovskites (including bandgap,
carrier masses, electron-contributed dielectric constant, and exciton binding
energy) as the function of the PF, indicating the weak coupling between
electronic properties and PF. Interestingly, the magnitudes of electron
effective masses for different structures distributes in a narrow range, while
the magnitudes of hole effective mass distribute a broad range, as well as
bandgap, electron-contributed dielectric constant, and exciton binding energy
results.

Fig. 6c illustrates the bandgap values of different halide perovskites as the
function of electronegativity difference between anions and cations (e.g., the
difference between X and B site in ABX3 halide perovskites). For the ionic
crystals, the electronegativity difference between anions and cations was
regarded as a good descriptor of the bandgap. Taking rock salt structures
(NaX, X=F, Cl, Br, I) as an example, the bandgap decreases monotonically with
the electronegativity difference. Indeed the bandgap values of halide
perovskites monotonically depend on the electronegativity difference, though
the magnitudes of bandgap for single, double, and low-dimensional perovskites
distributes in a broad range. Fig. 6d shows the optical transition intensity
for different candidate structures as function of band gap, where no obvious
tendency could be observed. The above analyses indicate that understanding the
compounds-structure-property relationship is a huge challenge that need
researchers’ sharp intuition and rich experience.

Figure 6: Calculated electronic and optical properties of halide perovskites
(HPs) and derivatives. (a) The number of calculated compounds, including
hybrid organic-inorganic perovskites, all-inorganic perovskites, and
perovskites-like structures. Orange represents the singe perovskites. Cyan
indicates the double perovskites, and yellow represents the low-dimensional
perovskite-like halides. (b) Calculated electronic properties (bandgap (Eg in
eV), electron and hole effective mass (me and mh), dielectric constant
($\varepsilon$), and exciton binding energy (Eeb in eV)) versus the atomic
packing factor. The inorganic atomic radii values were taken from SlaterSlater
and KhanKhan and Tuck (1985), the equivalent atomic radii of organic molecules
were taken from Ref. Yang  _et al._ (2017). (c) Calculated bandgap (Eg) versus
the electronegative difference between halogen anions and B-site cations. The
electronegativity values were taken from Pauling’s scaleLide (2004). (d)
Inter-band optical transition intensity versus optical bandgap
(E${{}_{g}}^{optical}$).

To facilitate the next machine-learning study and future usage of the data, we
stored our calculated data in the database. The stored calculation data
contain the relaxed structure and calculation instance. The calculation
instance represents a series of information, including input parameters (such
as initial structure and calculation parameters) and extracted output
information, such as relaxation structure, calculation time, calculated output
energy band data, effective mass, dielectric constant, optical absorption
spectrum, and exciton binding energy. In the database, there is a one-to-many
relationship between the structure model and the calculation instance. Users
can access the entry properties of structure class to obtain past calculation
data.

Machine-learning algorithms combined with big data of materials allow
computers to automatically learn the structure-property relationship without
human intervention or assistance. Here, considering bandgap as an example, we
call the gradient boosting regression tree (GBRT) model via the model
construction module of JAMIP, which is a tree-based integrated supervised
learning algorithm that has many applications in the materials fieldLu  _et
al._ (2018); Im  _et al._ (2019). Here, based on calculated halide
perovskites, a total of 60 features (Table S2) were generated by calling the
JAMIP database to produce the initial dataset. Of these features, the crystal
system was encoded by calling the JAMIP one-hot encoding method. Using
statistics of the target values, we found the bandgaps calculated by the
density functional theory (DFT) methods obey Gaussian distribution
($\mu$=1.89, $\sigma$=0.91). To evaluate the effectiveness of the training
model, we divided the dataset into 80% for training and 20% for testing. In
the training process, we used Bayesian optimization combined with the cross-
validation approach to finding the optimal hyperparameters. The coefficient of
determination (R2), mean absolute error (MAE), and root mean square error
(RMSE) were then used to evaluate the fitting quality of the hyperparameters
by the JAMIP’s model evaluation module. After several iterations, the GBRT
model achieved high fitting accuracy, with 0.9256 of R2, 0.1750 of MAE, and
0.2526 of RMSE (Fig. 7b). This indicates by choosing an appropriate machine
learning model and enough features, one can achieve a rather accurate
description of the bandgap values in the current halide perovskites materials.
Most importantly, we then used the JAMIP tree-based model feature selection
method to mine the features that are most relevant to the target property
bandgap, as shown in Fig. 7c. We found that the packing factor plays the most
important role in determining the values of bandgaps. The element information
related to B-site is more important than A/X-site. This is consistent with the
previously reported results. We also used the Pearson correlation coefficient
matrix to identify the correlations between these features (Fig. 7d). We found
that most features exhibit low correlation character. This implicates that
alternative effective features may be needed to reach a good description of
the electronic structures of the current complex materials system.

Figure 7: Learning the physical rules related to bandgaps. (a) Histogram of
bandgaps of 696 structures in HT calculations. (b) Comparison of the bandgap
predicted by the GBRT model with respect to the bandgap by DFT calculations.
(c) Top 30 important features ranked by GBRT algorithm. (d) Heat map of
Pearson correlation coefficient matrix among the 30 most important features.
The blue colors indicate a negative correlation, while the orange colors
represent a positive correlation.

##  IV Conclusion and discussion

To summarize, we have reported the development of an artificial-intelligence-
aided data-driven infrastructure named Jilin Artificial-intelligence aided
Materials-design Integrated Package (JAMIP), which is designed purposely to
meet the requirements of the studies of computational materials informatics.
It is an open-source software implemented mainly in Python language. With the
emphasis on automation, extensibility, reliability, and intelligence, it is
integrated by intimately connected units of Data generation (e.g., high-
throughput materials calculations and screening), Data collection (e.g.,
automatic data extraction, management, and storage), and Data learning (e.g.,
integrated descriptors engineering, and machine learning-based data mining
functions). We integrated the essential modules into the software for
computational materials informatics including materials production factory,
high-throughput first-principles calculations engine, computational tasks
monitoring progress, data extraction, management and storage system, and
artificial intelligence data mining and machine learning function. The
materials production factory consists of specially designed material structure
manipulation operators and an inorganic crystal structure prototype database,
which facilitates batch processing of materials production for high-throughput
calculations. In the high-throughput first-principles calculations engine we
have packaged various material property calculation functions ranging from
thermodynamic, geometrical, electronic, mechanical, phononic, photonic to
magnetic properties. This is expected to be beneficial to diverse functional
materials discovery and design. We have designed a computational task
allocation and monitoring system to guarantee efficient and reliable large-
scale calculations. We have included in the software the computational data
extraction and analysis toolkits, as well as a deliberately developed database
for data storage based on the Django framework. We have integrated the
artificial intelligence machine-learning based data mining functions that
cover data preprocessing, material feature engineering, machine learning model
construction, and performance evaluation. We demonstrated how our software can
be used in performing materials informatics studies of optoelectronic
semiconductors by taking halide perovskites as an instance.

The challenges still faced in the computational materials informatics field
are motivating us to further develop our software. Firstly, the insufficiency
in computational materials data is usually one of the most severe challenges
in data mining studies. To relieve this challenge, we are developing more
material structure manipulation operators into the materials production
factory, and meanwhile are working on an efficient structure search algorithm
based on machine learning approaches. Secondly, the huge computational cost is
a serious problem during data generation through high-throughput first-
principles calculations, especially for complex systems with hundreds or even
thousands of atoms. On this aspect, we are launching in the infrastructure our
specially developed machine-learning potential for energetic calculations with
good transferability. We have also been developing efficient and meanwhile
accurate electronic structure calculation approach based on machine learning.
Thirdly, seeking effective descriptors in machine-learning based studies is
one of the most important subjects in computational materials informatics. We
are developing alternative descriptors to accurately represent the structural,
electronic structure, and spectral, etc. properties of materials. Fourthly,
evaluating the stability of designed materials is an essential step when new
materials were designed theoretically and call for experimental verification.
In addition to the thermodynamic criteria implemented in the software, we are
developing other approaches associated with kinetic processes. This involves
finding effective algorithms to calculate the transition barriers of chemical
reaction or phase transformation, as well as the algorithms to simulate the
crystal growth procedure. We believe that with continuous development and
improvement, our open-source JAMIP software represents an important addition
to the facilities/infrastructures needed by computational materials
informatics, and will contribute to accelerating new materials discovery and
design from the computational side.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(Grants No. 61722403, 92061113 and 12004131). Calculations were performed in
part at the high-performance computing center of Jilin University.

## References

  * Agrawal and Choudhary (2016) A. Agrawal and A. Choudhary, APL Materials 4, 053208 (2016).
  * Schleder  _et al._ (2019) G. R. Schleder, A. C. M. Padilha, C. M. Acosta, M. Costa, and A. Fazzio, Journal of Physics: Materials 2, 032001 (2019).
  * (3) C. Draxl and M. Scheffler, in  _Handbook of Materials Modeling : Methods: Theory and Modeling_, edited by W. Andreoni and S. Yip, pp. 1–25.
  * (4) L. Himanen, A. Geurts, A. S. Foster, and P. Rinke, 6, 1900808.
  * Kalil and Wadia (2011) T. Kalil and C. Wadia, Retrieved from the World Wide Web April 21, 2014 (2011).
  * Ghiringhelli  _et al._ (2015) L. M. Ghiringhelli, J. Vybiral, S. V. Levchenko, C. Draxl, and M. Scheffler, Physical Review Letters 114, 105503 (2015).
  * Duvenaud  _et al._ (2015) D. Duvenaud, D. Maclaurin, J. Aguilera-Iparraguirre, R. Gómez-Bombarelli, T. Hirzel, A. Aspuru-Guzik, and R. P. Adams, arXiv e-prints , 1509.09292 (2015).
  * Schmidt and Lipson (2009) M. Schmidt and H. Lipson, Science 324, 81 (2009).
  * Rudy  _et al._ (2017) S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, Science Advances 3, e1602614 (2017).
  * Fundel  _et al._ (2007) K. Fundel, R. Küffner, and R. Zimmer, Bioinformatics 23, 365 (2007).
  * Tshitoyan  _et al._ (2019) V. Tshitoyan, J. Dagdelen, L. Weston, A. Dunn, Z. Rong, O. Kononova, K. A. Persson, G. Ceder, and A. Jain, Nature 571, 95 (2019).
  * Curtarolo  _et al._ (2012a) S. Curtarolo, W. Setyawan, S. Wang, J. Xue, K. Yang, R. H. Taylor, L. J. Nelson, G. L. W. Hart, S. Sanvito, M. Buongiorno-Nardelli, N. Mingo, and O. Levy, Computational Materials Science 58, 227 (2012a).
  * Zhao  _et al._ (2017) X.-G. Zhao, J.-H. Yang, Y. Fu, D. Yang, Q. Xu, L. Yu, S.-H. Wei, and L. Zhang, Journal of the American Chemical Society 139, 2630 (2017).
  * Luo  _et al._ (2021) S. Luo, T. Li, X. Wang, M. Faizan, and L. Zhang, WIREs Computational Molecular Science 11, e1489 (2021).
  * Curtarolo  _et al._ (2012b) S. Curtarolo, W. Setyawan, G. L. Hart, M. Jahnatek, R. V. Chepulskii, R. H. Taylor, S. Wang, J. Xue, K. Yang, O. Levy, M. J. Mehl, H. T. Stokes, D. O. Demchenko, and D. Morgan, Computational Materials Science 58, 218 (2012b).
  * Pizzi  _et al._ (2016) G. Pizzi, A. Cepellotti, R. Sabatini, N. Marzari, and B. Kozinsky, Computational Materials Science 111, 218 (2016).
  * Mathew  _et al._ (2016) K. Mathew, A. K. Singh, J. J. Gabriel, K. Choudhary, S. B. Sinnott, A. V. Davydov, F. Tavazza, and R. G. Hennig, Computational Materials Science 122, 183 (2016).
  * Jain  _et al._ (2015) A. Jain, S. P. Ong, W. Chen, B. Medasani, X. Qu, M. Kocher, M. Brafman, G. Petretto, G.-M. Rignanese, G. Hautier, D. Gunter, and K. A. Persson, Concurrency and Computation: Practice and Experience 27, 5037 (2015).
  * Yang  _et al._ (2018) X. Yang, Z. Wang, X. Zhao, J. Song, M. Zhang, and H. Liu, Computational Materials Science 146, 319 (2018).
  * Zhu  _et al._ (2020) W. Zhu, Y. Xu, J. Ni, G. Hu, X. Wang, and W. Zhang, Materials Science and Engineering: B 252, 114474 (2020).
  * Mathew  _et al._ (2017) K. Mathew, J. H. Montoya, A. Faghaninia, S. Dwarakanath, M. Aykol, H. Tang, I.-h. Chu, T. Smidt, B. Bocklund, M. Horton, J. Dagdelen, B. Wood, Z.-K. Liu, J. Neaton, S. P. Ong, K. Persson, and A. Jain, Computational Materials Science 139, 140 (2017).
  * Choudhary  _et al._ (2020) K. Choudhary, K. F. Garrity, A. C. E. Reid, B. DeCost, A. J. Biacchi, A. R. Hight Walker, Z. Trautt, J. Hattrick-Simpers, A. G. Kusne, A. Centrone, A. Davydov, J. Jiang, R. Pachter, G. Cheon, E. Reed, A. Agrawal, X. Qian, V. Sharma, H. Zhuang, S. V. Kalinin, B. G. Sumpter, G. Pilania, P. Acar, S. Mandal, K. Haule, D. Vanderbilt, K. Rabe, and F. Tavazza, npj Computational Materials 6, 173 (2020).
  * Wang  _et al._ (2021) G. Wang, L. Peng, K. Li, L. Zhu, J. Zhou, N. Miao, and Z. Sun, Computational Materials Science 186, 110064 (2021).
  * Pickard and Needs (2011) C. J. Pickard and R. J. Needs, Journal of Physics: Condensed Matter 23, 053201 (2011).
  * Hachmann  _et al._ (2011) J. Hachmann, R. Olivares-Amaya, S. Atahan-Evrenk, C. Amador-Bedolla, R. S. Sánchez-Carrera, A. Gold-Parker, L. Vogt, A. M. Brockway, and A. Aspuru-Guzik, The Journal of Physical Chemistry Letters 2, 2241 (2011).
  * Kirklin  _et al._ (2015) S. Kirklin, J. E. Saal, B. Meredig, A. Thompson, J. W. Doak, M. Aykol, S. Rühl, and C. Wolverton, npj Computational Materials 1, 15010 (2015).
  * Draxl and Scheffler (2019) C. Draxl and M. Scheffler, Journal of Physics: Materials 2, 036001 (2019).
  * Gorai  _et al._ (2016) P. Gorai, D. Gao, B. Ortiz, S. Miller, S. A. Barnett, T. Mason, Q. Lv, V. Stevanovic, and E. S. Toberer, Computational Materials Science 112, 368 (2016).
  * Glass  _et al._ (2006) C. W. Glass, A. R. Oganov, and N. Hansen, Computer Physics Communications 175, 713 (2006).
  * Lonie and Zurek (2011) D. C. Lonie and E. Zurek, Computer Physics Communications 182, 372 (2011).
  * Wang  _et al._ (2012) Y. Wang, J. Lv, L. Zhu, and Y. Ma, Computer Physics Communications 183, 2063 (2012).
  * Tipton and Hennig (2013) W. W. Tipton and R. G. Hennig, Journal of Physics: Condensed Matter 25, 495401 (2013).
  * Schmidt  _et al._ (2019) J. Schmidt, M. R. G. Marques, S. Botti, and M. A. L. Marques, npj Computational Materials 5, 1 (2019).
  * Liu  _et al._ (2017) Y. Liu, T. Zhao, W. Ju, and S. Shi, Journal of Materiomics High-throughput Experimental and Modeling Research toward Advanced Batteries, 3, 159 (2017).
  * Ramprasad  _et al._ (2017) R. Ramprasad, R. Batra, G. Pilania, A. Mannodi-Kanakkithodi, and C. Kim, npj Computational Materials 3, 1 (2017).
  * Kusne  _et al._ (2020) A. G. Kusne, H. Yu, C. Wu, H. Zhang, J. Hattrick-Simpers, B. DeCost, S. Sarker, C. Oses, C. Toher, S. Curtarolo, A. V. Davydov, R. Agarwal, L. A. Bendersky, M. Li, A. Mehta, and I. Takeuchi, Nature Communications 11, 5966 (2020).
  * Ouyang  _et al._ (2018) R. Ouyang, S. Curtarolo, E. Ahmetcik, M. Scheffler, and L. M. Ghiringhelli, Physical Review Materials 2, 083802 (2018).
  * Bartel  _et al._ (2019) C. J. Bartel, C. Sutton, B. R. Goldsmith, R. Ouyang, C. B. Musgrave, L. M. Ghiringhelli, and M. Scheffler, Science Advances 5, eaav0693 (2019).
  * Lu  _et al._ (2018) S. Lu, Q. Zhou, Y. Ouyang, Y. Guo, Q. Li, and J. Wang, Nature Communications 9, 3405 (2018).
  * Chen  _et al._ (2015) C. Chen, D. Chen, and F. Ciucci, Physical Chemistry Chemical Physics 17, 7831 (2015).
  * Tran and Ulissi (2018) K. Tran and Z. W. Ulissi, Nature Catalysis 1, 696 (2018).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. d. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Franceschetti and Zunger (1999) A. Franceschetti and A. Zunger, Nature 402, 60 (1999).
  * Zhang  _et al._ (2013) L. Zhang, J.-W. Luo, A. Saraiva, B. Koiller, and A. Zunger, Nature Communications 4, 2396 (2013).
  * Zhang  _et al._ (2012) L. Zhang, M. d’Avezac, J.-W. Luo, and A. Zunger, Nano Letters 12, 984 (2012).
  * (47) J. C. Slater,  41, 3199, publisher: American Institute of Physics.
  * Khan and Tuck (1985) M. A. Khan and D. G. Tuck, Inorganica Chimica Acta 97, 73 (1985).
  * Yang  _et al._ (2017) D. Yang, J. Lv, X. Zhao, Q. Xu, Y. Fu, Y. Zhan, A. Zunger, and L. Zhang, Chemistry of Materials 29, 524 (2017).
  * Lide (2004) D. R. Lide,  _CRC handbook of chemistry and physics_ , Vol. 85 (2004).
  * Im  _et al._ (2019) J. Im, S. Lee, T.-W. Ko, H. W. Kim, Y. Hyon, and H. Chang, npj Computational Materials 5, 1 (2019).


# Algorithms for Batch Hierarchical Reinforcement Learning

Tiancheng Zhao and Mohammad Gowayyed
Language Technologies Institute, Carnegie Mellon University
{tianchez<EMAIL_ADDRESS>
###### Abstract

Hierarchical Reinforcement Learning (HRL) exploits temporal abstraction to
solve large Markov Decision Processes (MDP) and provide transferable subtask
policies. In this paper, we introduce an off-policy HRL algorithm:
Hierarchical Q-value Iteration (HQI). We show that it is possible to
effectively learn recursive optimal policies for any valid hierarchical
decomposition of the original MDP, given a fixed dataset collected from a flat
stochastic behavioral policy. We first formally prove the convergence of the
algorithm for tabular MDP. Then our experiments on the Taxi domain show that
HQI converges faster than a flat Q-value Iteration and enjoys easy state
abstraction. Also, we demonstrate that our algorithm is able to learn optimal
policies for different hierarchical structures from the same fixed dataset,
which enables model comparison without recollecting data.

##  1 Introduction

Conventional tabular reinforcement learning is bottle-necked by the curse of
dimensionality for practical applications. The number of parameters that needs
to be trained grows exponentially with respect to the size of states and
actions. In order to make reinforcement learning practically tractable, one
line of research is hierarchical reinforcement learning (HRL), which develops
principled ways of temporal and state abstraction to reduce the dimensionality
for sequential decision making.

The basic idea of temporal abstraction is to develop macro-actions that take
several steps to terminate before returning. Usually good macro-actions aim to
solve sub-goals, so that multiple macro-actions divide difficult tasks into
several simpler ones. In addition, state abstraction tries to reduce the
dimensionality by removing irrelevant state variables for decision making,
reducing the cardinality of state space and helping in tackling over-fitting.
These two techniques lead to natural hierarchical control architecture, which
intuitively resembles how humans solve complex tasks.

Another area of research closely related to our work is batch reinforcement
learning. Batch reinforcement learning aims to learn the best policy from a
fixed set of prior-known samples. Compared to on-policy algorithms, batch
reinforcement learning enjoys stability and data-efficiency. More importantly,
it allows to apply reinforcement learning in a practical problem that is
expensive in collecting new samples, such as education, spoken dialog system
and medical systems. Well-known algorithms in batch reinforcement learning
include Least Square Policy Iteration (LSPI) lagoudakis2003least , Fitted Q
iteration (FQI) ernst2005tree , Neural Fitted Q Iteration (NFQ)
riedmiller2005neural  and etc.

In this paper, we combine batch learning and hierarchical reinforcement
learning, in order to achieve faster learning speed, data efficiency and model
comparison.

##  2 Related Work

There are three major approaches developed relatively independently
barto2003recent , aiming to formalize the idea of abstraction into
reinforcement learning. The three approaches are: 1) the option framework
sutton1999between , 2) Hierarchies of Abstract Machines (HAMs)
parr1998reinforcement  and 3) MAXQ framework dietterich2000overview .

Under the option framework, the developers augment the original action set by
options, which are macro actions that have their own predefined policy,
termination state and active state. Sutton et al have shown that such a system
is a semi-Markov Decision Process (SMDP), which converges to a unique
hierarchical optimal solution using a modified Q-learning algorithm. For HAM
framework, rather than giving out the entire policy of these macro actions,
developers just need to specify a partial program that specifies a part of the
policy. Using HAMQ learning parr1998reinforcement , HAM can also converge to a
hierarchical optimal solution.

At last, MAXQ framework provides an elegant formulation that decomposes the
original MDP into several subroutines in a hierarchy and the algorithm can
learn policies recursively for all the subroutines. Therefore, in the MAXQ
framework, there is no need to specify the policy for any macro-actions.
However, Dietterich shows that it can only achieve recursive optimal solution,
which in the extreme case, can be arbitrarily worse than the hierarchical
optimal solution.

All of the above work assumes that the agent can interact with the world while
learning. However, in real-world applications that needs HRL, it is usually
very expensive to collect data and terrible failures are not allowed on
operation. This forbids the usage of online learning algorithms that could
potentially preform horribly in the early learning stage. To our best
knowledge, there is little prior work bloch2011reducing  in developing batch
learning algorithms that allow a hierarchical SMDP to be trained from an
existing data set collected from a stochastic behavior policy. We believe that
such algorithms are valuable for applying HRL in complex practical domains.

##  3 Batch Learning for HSMDP

###  3.1 Definitions

Mostly, we follow the definitions in the MAXQ framework. However, for notation
simplicity, we also borrow some notations from the $option$ framework.

###  3.2 Markov Decision Process

An MDP $M$ is described by $(S,A,P,R,P_{0})$

  1. 1.

$S$ is the state space of $M$

  2. 2.

$A$ is a set of primitive actions that are available

  3. 3.

$P(s^{\prime}|s,a)$ defines the transition probability of executing primitive
action $a$ in state, $s$

  4. 4.

$R(s^{\prime}|s,a)$ is the reward function defined over $S$ and $A$

###  3.3 Hierarchical Decomposition

An MDP, $M$, can be decomposed into a finite set of subtasks
$O=\\{O_{0},O_{1}...O_{n}\\}$ with the convention that $O_{0}$ is the root
subtask, i.e. solving $O_{0}$ solves the entire original MDP, $M$. $O_{i}$ is
then a Semi-Markov Decision Process (SMDP) that shares the same $S$, $R$, $P$
with $M$, and has an extra tuple $<\beta_{i},U_{i}>$, where:

  1. 1.

$\beta_{i}(s)$ is the termination predicate of subtask $O_{i}$ that partition
$S$ into a set of active states, $S_{i}$ and a set of terminal states $T_{i}$.
If $O_{i}$ enters a state in $T_{i}$, $O_{i}$ and its subtasks exit
immediately, i.e. $\beta_{i}(s)=1$ if $s\in T_{i}$, otherwise
$\beta_{i}(s)=0$.

  2. 2.

$U_{i}$ is a nonempty set of actions that can be performed by $O_{i}$. The
actions can be either primitive actions from $A$ or other subtask, $O_{j}$,
where $i\neq j$. We will refer to $U_{i}$ as the children of subtask $O_{i}$.

It is evident that a valid hierarchical decomposition forms a direct acyclic
graph (DAG) where each non-terminal node corresponds to a subtask, and each
terminal node corresponds to a primitive action. For later discussion, we will
use hierarchical decomposition and DAG interchangeably.

###  3.4 Hierarchical Policy

A hierarchical policy, $\pi$, is a set of policies for each subtask, $O_{i}$,
$\pi=\\{\pi_{0},\pi_{1}...\pi_{n}\\}$. In the terminology of $option$
framework, a subtask policy is a deterministic $option$, with $\beta_{i}(s)=1$
for $s\in T_{i}$, and $0$ otherwise.

###  3.5 Recursive Optimality

A recursive optimal policy for MDP $M$ with hierarchical decomposition is a
hierarchical policy $\pi=\\{\pi_{0}...\pi_{n}\\}$, such that for each subtask,
$O_{i}$, the corresponding policy $\pi_{i}$ is optimal for the SMDP defined by
the set of states, $S_{i}$, the set of actions $U_{i}$, the state transition
probability $P^{\pi}(s^{\prime},N|s,a)$, and the rewards function
$R(s^{\prime}|s,a)$.

##  4 Algorithm

The problem formulation is as following: given any finite set of samples,
$D=\big{\\{}(s_{m},a_{m},r_{m},s^{{}^{\prime}}_{m})|m=1,2,...,M\big{\\}}$ and
any valid hierarchical decomposition $O$ of the original MDP $M$, we wish to
learn the recursive optimal hierarchical policy $\pi^{*}$.

We now propose Hierarchical Q-value Iteration, HQI, and we prove that it
converges to the recursive optimal solution for any hierarchical decomposition
given that the batch sample distribution has sufficient state action
exploration. The basic idea is to train every subtask using Subtask Q-value
Iteration (SQI) in a bottom up fashion. The training prerequisite of SQI for a
specific subtask $O_{i}$ is that all of its children $U_{i}$ have converged to
their greedy optimal policies. In order to fulfil this constraint, HQI first
topologically sorts the DAG and running SQI from subtasks whose children have
only primitive actions. After those subtasks converge to their optimal policy,
the algorithm continues to other subtasks whose children are either converged
or primitive actions. We will show that there always exist an ordering of
training every subtask in a valid DAG that fulfills the prerequisite of SQI.

One challenge of training a subtask with subtask children is that we cannot
use the optimal SMDP Bellman equation described in the MAXQ framework
dietterich2000overview , $Q_{i}(s,u)$, which is the Q-value function for
subtask, $O_{i}$, at state, $s$ and action $u$:

$Q_{i}(s,u)=V(s,u)+\sum_{s^{\prime},N}P_{i}^{\pi}(s^{\prime},N|s,u)\gamma^{N}Q_{i}^{\pi}(s^{\prime},\pi_{i}(s^{\prime}))\\\
$ (1) $V(s,u)=\begin{cases}max_{u^{\prime}}(Q_{u}(s,u^{\prime}))&\text{$u$ is
subtask}\\\ \sum_{s^{\prime}}P(s^{\prime}|s,u)R(s^{\prime}|s,u)&\text{$u$ is
primitive}\end{cases}$ (2)

The main problem of this equation(1) is that in order to estimate the Q-value
for a subtask children, $u$, the parent $O_{i}$ needs to have an estimate
about the transition probability $P_{i}^{\pi}(s^{\prime},N|s,u)$, which is the
distribution of $u$’s exit state and number of primitive steps needed to reach
its termination. Although the termination states of the child $u$ are given by
$T_{u}$, it is difficult to estimate the joint distribution of termination
steps $N$ and $s$’ if $u$ follows an policy that is different from the
behavior policy without recollecting new samples. This is because since the
behavior policy is usually random and poor in performance, the collected
samples do not provide information about how many steps the subtask $u$ would
take to terminate if following a different (optimal) policy.

Therefore, instead of using the above Bellman equation that updates the Q
table of the parent when a child exits, we use the intra-option Bellman
equation proposed in the $option$ framework sutton1999between :

$\displaystyle Q_{i}(s,u)$ $\displaystyle=\sum_{a\in
A}{\pi_{u}(a|s)E[r(s,a)+\gamma V_{i}(s^{\prime},u)]}$ (3)
$\displaystyle=\sum_{a\in
A}{\pi_{u}(a|s)\bigg{[}r(s,a)+\gamma\sum_{s^{\prime}}P(s^{\prime}|s,a)V_{i}(s^{\prime},u)\bigg{]}}$
(4)

Where

$V_{i}(s,u)=(1-\beta_{i}(s))Q_{i}(s,u)+\beta_{i}(s)\max_{u^{\prime}\in
U_{i}}Q_{i}(s,u^{\prime})$ (5)

Equation (3) also yeilds a contraction in the max norm and is able to learn
the Q table after observing every new reward, which eliminates the need to
estimate $P_{i}^{\pi}(s^{\prime},N|s,u)$. Another key benefit is that we can
use flat samples to estimate the one step transition probability and rewards
in equation (3), which makes the algorithm independent of the hierarchical
decomposition and is able to learn optimal polices for different structures
from the same dataset. Specifically, we can estimate the above two terms by
$\sum_{s^{\prime}}P(s^{\prime}|s,a)V(i,s^{\prime},u)\approx\frac{1}{c}\sum_{m=1}^{c}V(i,s^{\prime}_{m}=s^{\prime},u)$
and $r(s,a)\approx\frac{1}{c}\sum_{m=1}^{c}r(s_{m}=s,a_{m}=a)$, where $c$ is
the number of experiences that has $s^{\prime}$ and ($s$, $a$), respectively.
At last, since we assume converged subtasks follow deterministic greedy
policy, $\pi_{u}(a|s)=1$ if $a$ is the greedy primitive action that subtask
$u$ would take at state $s$, and $\pi_{u}(a|s)=0$ otherwise. This step is in
fact crucial for HQI to learn the optimal policy because it allows a subtask
to discard those samples that are not following the optimal behavior of its
children.

The HQI algorithm is summarized in Algorithm 1 and SQI is summarized in
Algorithm 2. Any dataset $D$ can be used at every iteration of SQI. If the
initial data is sufficient to cover important state-action space, the same
dataset is able to train all subtasks of the DAG.

###  4.1 Extension to Function Approximation and State Abstraction

We note that it is trivial to extend SQI to Fitted-SQI which uses a function
approximator to model the Q-value function for a subtask at the end of each
iteration. The direct advantage of using function approximation is that it can
incorporate powerful supervised regression methods, such as Gaussian Processes
or Neural Networks to scale up to large-scale and continuous MDPs. Although
using function approximations usually compromises the theoretical convergence
guarantee for tabular MDP, our experiments shows that Fitted-HQI is able to
converge to the unique optimal solution. The Fitted-SQI is summarized in
Algorithm 3.

Furthermore, state abstraction here means finding a subset of state variables
that are most informative for each subtask. A good hierarchical decomposition
decomposes the original MDP into several simpler ones, such that the agent
only needs to care about a small set of features in each task. Therefore, a
good structure should create easy opportunity for state abstraction at each
level of its hierarchy. Many techniques have been explored painter2012greedy
qin2014sparse geramifard2013batch  in non-hierarchical batch reinforcement
learning to achieve state abstractions. These methods can be directly applied
in fitting step of Fitted-SQI and each subtask can learn its own sparse state
representation. Due to the space limit, we conduct simple manual state
abstraction for each subtask in this paper, and leave the study of analyzing
automatic feature selection techniques into future works.

Algorithm 1 Hierarchical Q-value iteration (HQI)

0: $O$, $D$

$Train\leftarrow O_{i}\in O$ with only primitive children

$Done\leftarrow\\{A\\}$

while $Train\neq empty$ do

for $O_{i}\in Train$ do

SQI$(O_{i},D)$

$Done.add(O_{i})$

end for

$Train\leftarrow O_{i}\in(O-Done)$ AND $U_{i}\in Done$

end while

Algorithm 2 Subtask Q-value Iteration (SQI)

0: $O_{i},D$

while $k<maxIter$ do

for $(s,a,r,s^{\prime})\in D$ do

for $u\in U_{i}$ do

if $s^{\prime}\in T_{i}$ then

$y\leftarrow r$

$Q^{k}_{i}(s,u)\leftarrow(1-\alpha)Q_{i}^{k-1}(s,u)+\alpha y$

else

if GreedyPolicy($u$, $s$) $=a$ then

$y\leftarrow r+\gamma((1-\beta_{u}(s^{\prime}))Q_{i}^{k-1}(s^{\prime},u)$

$\quad\quad\quad+\beta_{u}(s^{\prime})max_{u^{\prime}\in
U_{i}}Q_{i}^{k-1}(s^{\prime},u^{\prime}))$

$Q_{i}^{k}(s,u)\leftarrow(1-\alpha)Q_{i}^{k-1}(s,u)+\alpha y$

end if

end if

end for

end for

end while

Algorithm 3 Fitted Subtask Q-value Iteration (Fitted SQI)

0: $O_{i},D$

while $k<maxIter$ do

$X\leftarrow[]$, $Y\leftarrow[]$

for $(s,a,r,s^{\prime})\in D$ do

for $u\in U_{i}$ do

if $s^{\prime}\in T_{i}$ then

$y\leftarrow r$

else

if GreedyPolicy($u$, $s$) $=a$ then

$y\leftarrow r+\gamma((1-\beta_{u}(s^{\prime}))Q_{i}^{k-1}(s^{\prime},u)$

$\quad\quad\quad+\beta_{u}(s^{\prime})max_{u^{\prime}\in
U_{i}}Q_{i}^{k-1}(s^{\prime},u^{\prime}))$

end if

end if

$X.add((s,u))$, $Y.add(y)$

end for

end for

$Q_{i}^{k}\leftarrow fit(X,Y)$

end while

Algorithm 4 GreedyPolicy

0: $u,s$

if $u\in A$ then

return u

else

$u^{*}=\arg\\!\max_{u^{\prime}\in U_{u}}Q_{u}(s,u^{\prime})$

return GreedyPolicy($u^{*}$, $s$)

end if

##  5 Proof of Convergence

In this section, we prove that our HQI (in the tabular case) converges to the
recursive optimal policy. Assume that the policy at each subtask $M_{i}$ is
ordered, such that it break ties deterministically (e.g favor left to right),
it defines a unique recursive optimal hierarchical policy, $\pi^{*}_{r}$, and
a corresponding recursive optimal Q function $Q^{*}_{r}$. We then show that
HQI converge to $\pi^{*}_{r}$ and $Q^{*}_{r}$. The $r$ subscript refers to
recursive optimality.

###  5.1 Proof

We want to prove that for an MDP $M=(S,A,P,R,P_{0},\gamma)$ with hierarchical
decomposition $O=\\{O_{0},..O_{n}\\}$, HQI converges to recursive optimal
policy for the hierarchical policy of $M$, $\pi^{*}_{r}$.

  1. 1.

We first prove that: For a subtask $O_{i}$, with all of its children converged
to their recursive optimal policies and infinity amount of batch data,
algorithm SQI converge to the optimal Q-value function after infinity number
of iterations, $\lim_{k->\infty}Q^{k}_{i}=Q^{*}_{i}$

  2. 2.

We then show that HQI provides an order of training all the subtasks in the
DAG graph, such that when training a subtask, $O_{i}$, all of its children,
$U_{i}$ already converged to their optimal recursive policies.

###  5.2 Proof Sketch

Step 1 We begin with the base case for subtask whose children are all
primitive actions. We can notice that equation (3) falls back to traditional
Bellman operator for flat MDP, because a primitive action always terminate
after one step:

$Q(i,s,u)=r(s,a)+\gamma\sum_{s^{\prime}}P(s^{\prime}|s,a)\max_{u^{\prime}\in
U_{i}}Q(i,s,u^{\prime})$ (6)

Therefore, for subtask with only primitive children, SQI is equivalent to flat
Q-value iteration, which is guaranteed to converge to optimal policy given
sufficient data.

Then for subtasks with other subtask children, by definition, when we run SQI,
the children of $O_{i}$ ($U_{i}$) have converged to their unique deterministic
optimal recursive policy. This means that every action $u\in U_{i}$, is a
deterministic deterministic Markov option as defined in the $option$ framework
sutton1999between . sutton1999between  proved that ”for any set of
deterministic Markov options one step intra-option Q-learning converges w.p. 1
to the optimal Q-values, for every option regardless of what options are
executed during learning provided every primitive action gets executed in
every state infinitely often”. Refer to the sutton1999between , for the
detailed proof.

Step 2 By definition, a hierarchical decomposition is a Directed Acyclic Graph
(DAG) with edges from parents to their children. In this proof, we first
reverse the edges so that they are from children to their parents. Also we
know from Graph Theory that any Directed Acyclic Graph has at least one
topological sort, such that every edge $uv$, $u$ comes before $v$ in the
ordering cormen2001section . Therefore, we can to topologically sort the
hierarchical decomposition with reversed edges such that SQI can always train
the children before parents.

Also, the definition of topological sort ensures the initial condition that
there is at least one subtask that only has primitive children. Therefore, we
can then conclude that for any DAG, HQI can traverse the subtasks such that
the conditions of SQI convergence are met. Then, HQI converges for all
subtasks.

##  6 Experiments

###  6.1 Experimental Setup

We applied our algorithm to the Taxi domain described in sutton1999between .
This is a simple grid world that contains a taxi, a passenger, and four
specially-designated locations labeled R, G, B, and Y. In the starting state,
the taxi is in a randomly-chosen cell of the grid, and the passenger is at one
of the four special locations. The passenger has a desired destination that
he/she wishes to reach, and the job of the taxi is to go to the passenger,
pick him/her up, go to the passenger’s destination, and drop the passenger.
The taxi has six primitive actions available to it: move one step to one of
the four directions (north, south, east and west), pick up the passenger and
put down the passenger. To make the task more difficult, the move actions are
not deterministic, so that it has $20\%$ chance of moving in one of the other
directions. Also, every move in the grid will cost $-1$ reward. Attempting to
pick up or drop passenger at wrong location will cause $-10$ reward. At last,
successfully finish the task has $20$ reward. The grid is described in figure
1. Therefore, there are 4 possible state for the destination, 5 possible state
for the passenger (4 location and 5 is on the car), 25 possible locations,
which results into $500*6=3000$ parameters in the Q-table that needs to be
learned. We denote the state variable as $[dest,pass,x,y]$ for later
discussion.

The dataset for each run were collected in advance by choosing actions
uniformly at random with different sizes. We evaluate the performance of
algorithms by running greedy execution for $100$ times to obtain average
discounted return at every $5000$ new samples and up to $60,000$ samples. We
repeat the experiments for $5$ times to evaluate the influence of different
sample distribution. The discounting factor is set to be $0.99$.

Figure 1: The Taxi Domain Figure 2: DAG 1

###  6.2 Results

We conducted three sets of experiments: 1) comparison of HQI with flat Q-value
Iteration and the effect of state abstraction. 2) learning polices for
different DAGs from the same dataset and 3) learning policy using Fitted-HQI
with Random Forest as the function approximator.

The first experiment compares HQI against flat Q-value Iteration (FQI). Also,
as pointed out in dietterich2000hierarchical , state abstraction is essential
for MAXQ to have fast learning speed compared to flat Q learning. As a result,
we manually conduct state abstraction for each subtask in DAG $1$. However,
different from the aggressive state abstraction described in
dietterich2000hierarchical , where every subtask and child pair has a
different set of state variables, we only conduct a simple state abstraction
at subtask level, i.e. all children of a subtask has the same state
abstraction. The final state abstraction is listed in Table 1. As described
above, we run $5$ independent runs with different random samples of different
sizes, we report the mean average discounted return over five runs in Figure
3, as well as the best average discounted reward of the five runs in Figure 4.

Table 1: DAG-1 State Abstraction subtask | active states
---|---
root | [pass]
get | [pass x y]
put | [dest x y]
navi_get | [pass x y]
navi_put | [dest x y]

Results show that both HQI with and without state abstraction consistently
outperforms the FQI when there is limited training data. When the dataset is
large enough, they all converge to the same optimal performance, which is
around $1.0$. We also notice that, occasionally, HQI with state abstraction
can learn the optimal performance state abstraction with very limited samples,
i.e $5000$ samples. This demonstrates that with proper hierarchy constraints
and good behavioral policy, HQI can generalize much faster than FQI. Moreover,
even the HQI without state abstraction consistently outperforms FQI in terms
of sample efficiency. This is different from the behavior of the on-policy
MAXQ-Q algorithm reported in dietterich2000hierarchical , which needs state
abstraction in order to learn faster than Q-learning. We argue that HQI
without state abstraction is more sample efficient than FQI for the following
reasons: 1) HQI uses all applicable primitive samples to update the Q-table
for every subtask while MAXQ-Q only updates for the subtask that executes that
particular action. 2) Upper level subtask in MAXQ-Q needs to wait for its
children gradually converges to their greedy optimal policy before it can have
have a good estimate of $P(s^{\prime},N|s,u)$ while HQI does not have this
limitation.

Figure 3: Average discounted reward. The SA postfix means state abstraction.
The error bar is one standard deviation over 5 runs. Figure 4: Best
performance comparison

The second experiment is running HQI on different variations of hierarchical
decomposition of the original MDP. Figure 5 and Figure 6 show two different
valid DAGs that could also solve the original MDP. Figure 7 demonstrates that
with sufficient data all three DAG converge to their recursive optimal
solution, which confirms that HQI is able converge for different hierarchies.
In terms of sample efficiency, three structures demonstrate slight different
behavior. We can notice that DAG $2$ learns particularly slower than the other
two. We argue that this is because of poor decomposition of the original MDP.
Based on the problem settings, pick and drop are all risky actions (illegal
execution lead to $-10$ reward), while in DAG $2$ these two actions are mixed
with low-cost move actions while the other two DAGs isolated them in a higher
level of decision making. Therefore, designing good hierarchy is crucial to
obtain performance gain versus flat RL approaches. This emphasizes the
importance of the off-policy nature of HQI, which allows developers to
experiment with different DAG structures without collecting new samples. How
to effectively evaluate the performance of particular hierarchical
decomposition without using a simulator is a part of our future research.

Figure 5: DAG 2 Figure 6: DAG 3

The last experiment utilizes Random Forests as the function approximator to
model the Q-value function in DAG 1. The main purpose is to demonstrate the
convergence of Fitted HQI. For each subtask $O_{i}$ the Q-value function
$Q_{i}(s,u)$ is modelled by a random forest with $[dest,pass,x,y]$ as the
input feature. Since $dest$ and $pass$ are categorical variables, we represent
them as a one-hot vector, which transforms the state variable into a $11$
dimension vector (4d for destination, 5d for passenger and 2d for the $x,y$
coordinate). We report the mean average discounted rewards over 5 independent
runs with different random samples of different sizes. Figure 8 shows that
Fitted-HQI achieves similar performance compared to Tabular HQI.

Figure 7: Comparison of different DAGs. The error bar is one standard
deviation over 5 runs.

###  6.3 Comparison with MAXQ-Q and Intra-Option Learning

Compared to MAXQ-Q, HQI enjoys sample efficiency and the ability to be off-
policy. The advantage of off-policy is that it does not require hyper-
parameter tuning such as exploration rate. Since high level subtask in MAXQ-Q
needs to wait for its children converge first, developers usually set a faster
exploration decay rate for lower level subtasks, which is an extra
hyperparameter that needs tuning. The limitation of HQI is that it maintains
an independent Q table for each subtask, while MAXQ-Q allows a part of the
parent value function recursively retrieved from its children, a technique
known as value function decomposition dietterich2000hierarchical . This allows
more compact memory usage and accelerates the learning of the parents. How to
share value function in the off-policy setting is a future research topic.

For intra-option learning in option framework, the main advantage of HQI is
that it does not require developers to fully define the policy of every
options. Instead, one only needs to define a DAG with terminal predicate for
each node in the graph. We argue that in general it is easier to define a task
hierarchy than giving a full policy for macro-actions. Therefore, HQI combines
the strength of intra-option off-policy learning with MAXQ. bloch2011reducing
provides a method of training options in an off-policy fashion. Compared to
it, HQI has the advantage of learning all subtasks from flat batch dataset, so
that our algorithm does not require a task DAG priory to collecting data and a
manual definition of option policies.

Figure 8: Comparison between fitted-HQ and HQI on DAG-1. The error bar is one
standard deviation over 5 runs.

##  7 Conclusion and Future Work

In this paper, we introduced an off-policy batch learning algorithm for
hierarchical RL. We showed that it is possible to blindly collect data using a
random flat policy. Then, we use this data to learn different structures that
the data collection were not aware of. Our experiments on the Taxi domain show
that it converges faster than FQI to the optimal policy. It also shows that
different DAG structures are able to learn from this flat data, with different
speeds. Every DAG structure has its own number of parameters, which suggests a
possible line for research to try to minimize the number of parameters in the
hierarchy. Other future work include comparing different feature selections
techniques for Fitted-SQI and applying the algorithm to large-scale and
complex domains.

## References

  * (1) Andrew G Barto and Sridhar Mahadevan.  Recent advances in hierarchical reinforcement learning.  Discrete Event Dynamic Systems, 13(1-2):41–77, 2003.
  * (2) Mitchell Keith Bloch.  Reducing commitment to tasks with off-policy hierarchical reinforcement learning.  arXiv preprint arXiv:1104.5059, 2011.
  * (3) Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, and Clifford Stein.  Section 22.4: Topological sort.  Introduction to Algorithms (2nd ed.), MIT Press and McGraw-Hill, pages 549–552, 2001.
  * (4) Thomas G Dietterich.  Hierarchical reinforcement learning with the maxq value function decomposition.  J. Artif. Intell. Res.(JAIR), 13:227–303, 2000.
  * (5) Thomas G Dietterich.  An overview of maxq hierarchical reinforcement learning.  In Abstraction, Reformulation, and Approximation, pages 26–44. Springer, 2000.
  * (6) Damien Ernst, Pierre Geurts, and Louis Wehenkel.  Tree-based batch mode reinforcement learning.  In Journal of Machine Learning Research, pages 503–556, 2005.
  * (7) Alborz Geramifard, Thomas J Walsh, Nicholas Roy, and Jonathan How.  Batch-ifdd for representation expansion in large mdps.  arXiv preprint arXiv:1309.6831, 2013.
  * (8) Michail G Lagoudakis and Ronald Parr.  Least-squares policy iteration.  The Journal of Machine Learning Research, 4:1107–1149, 2003.
  * (9) Christopher Painter-Wakefield and Ronald Parr.  Greedy algorithms for sparse reinforcement learning.  arXiv preprint arXiv:1206.6485, 2012.
  * (10) Ronald Parr and Stuart Russell.  Reinforcement learning with hierarchies of machines.  Advances in neural information processing systems, pages 1043–1049, 1998.
  * (11) Zhiwei Qin, Weichang Li, and Firdaus Janoos.  Sparse reinforcement learning via convex optimization.  In Proceedings of the 31st International Conference on Machine Learning (ICML-14), pages 424–432, 2014.
  * (12) Martin Riedmiller.  Neural fitted q iteration–first experiences with a data efficient neural reinforcement learning method.  In Machine Learning: ECML 2005, pages 317–328. Springer, 2005.
  * (13) Richard S Sutton, Doina Precup, and Satinder Singh.  Between mdps and semi-mdps: A framework for temporal abstraction in reinforcement learning.  Artificial intelligence, 112(1):181–211, 1999.