The search for new chemical materials, e.g., catalytic materials for a plethora of chemical reactions, produces large amounts of data. To discover useful knowledge from those data, statistical as well as machinelearning data mining methods have been used in this area since the late 1990s, the former represented in particular by the analysis of variance, decision trees and support vector regression, the latter by main variants of feed-forward neural networks.
This paper summarizes experience from nearly ten years using and developing neural-networks based data mining methods for catalytic data. Artificial neural networks are the most popular regression model in this application area. In the survey [10], more than 20 published applications of multilayer perceptrons (MLPs) to catalytic data have been listed, as well as several applications of radial basis function networks. The role of feed-forward neural nets as a regression model predicting catalytic performance of materials (such as yield, conversion, selectivity) is due partially to their preceding success in other areas, but mostly to their ability ⋆ The research reported in this paper has been supported by the grant No. 201/08/1744 of the Grant Agency of the Czech Republic and partially supported by the Institutional Research Plan AV0Z10300504.
to serve as universal approximators in very general function spaces [12,14,18]. This ability is particularly valuable in the context of the highly nonlinear nature of the dependencies encountered in catalysis (cf. Figure 1). However, it seems to be little awareness, among researchers using artificial neural networks in catalysis, of the difference among the symbolic nature of the knowledge obtained from data by analysis of variance and decision trees, and the numeric nature of the knowledge obtained by neural network regression. Incited by the situation just outlined, the paper presents two strategies for the application of artificial neural networks to data about chemical materials. The first strategy relies on the numeric knowledge from neural network regression. Although numeric knowledge is much less understandable to humans than symbolic knowledge (in terms of [4], it has a high "data fit", but a low "mental fit"), we show that in this application area, it can be very useful if directly integrated with the optimization of materials performance in an approach called surrogate modelling. In that context, also the possibility to increase the accuracy of neural network regression by means of boosting is mentioned. The other strategy, on the other hand, relies on employing rules extraction methods to obtain, from trained neural networks, symbolic knowledge.
These two strategies determine also the structure of the paper. In Section 2, the surrogate modelling approach is described. Section 3 then explains a method for the extraction of logical rules from trained neural networks. Both strategies are in the respective sections illustrated using real-world examples.
From the point of view of theoretical computer science, the search for most suitable chemical materials entails complex optimization tasks. As objective functions, those tasks use various properties of the materials, e.g. in the case of catalytic materials, properties quantifying their catalytic performance, such as yield, conversion, or selectivity. A crucial feature of such objective functions is that they cannot be expressed analytically, their values must be obtained empirically.
For their optimization, it is not possible to employ most common optimization methods, such as steepest descent, conjugate gradient methods or the Levenberg-Marquardt method. Indeed, to obtain sufficiently precise numerical estimates of gradients or second order derivatives of the empirical objective function, those methods need to evaluate the function in points some of which would have a smaller distance than is the empirical error of catalytic measurements. That is why methods not requiring any derivatives have been used to solve such optimization tasks, such as the simplex method, and most frequently genetic and other evolutionary algorithms [2]. To compensate for missing information about derivatives, these methods need quite large number of objective function evaluations. In the context of catalysis, this is quite disadvantageous because the evaluation of the empirical objective functions used in the search for optimal catalysts is often costly and time-consuming. Testing a generation of catalytic materials proposed by an evolutionary algorithm typically needs several days of time and costs thousands of euros.
The usual approach to decreasing the cost and time of optimization of empirical objective functions is to evaluate the function only in points considered to be most important for the progress of the employed optimization method, and to evaluate its suitable regression model otherwise. That model is termed surrogate model of the function, and the approach is referred to as surrogate modelling [17,20,23,27]. Needless to say, the time and costs needed to evaluate a regression model are negligible compared to time and costs needed to evaluate empirical functions such as yield or conversion. However, it must not be forgotten that the agreement between the results obtained with a surrogate model and those obtained with the original function depends on the accuracy of the model.
The fact that feed-forward neural networks are the most frequent regression models in catalysis suggests them as the most natural candidate for surrogate models in this area. Indeed, several nice examples of the application of neural-network based surrogate modelling to the optimization of performance of catalytic materials have been published during the last five years [3,6,21,25]. Within the overall context of the application of artificial neural networks to mining catalytic data, however, they are still rare.
Although surrogate modelling has been also applied to conventional optimization [5], it is most frequently encountered in connection with evolutionary algorithms because for them, the approach leads to the approximation of the fitness function, whose usefulness in evolutionary computation is already known [13,19]. For the progress of evolutionary optimization, most important criteria are on the one hand points that indicate closeness to the global optimum (through highest values of the fitness function), on the other hand points that most contribute to the diversity of the population.
In the literature, various possibilities of combining evolutionary optimization with surrogate modelling have been discussed [17,24,27]. Nevertheless, all of them are controlled by one of two basic approaches:
A. The individual-based-control consists in choosing between the evaluation of the empirical objective function and the evaluation of its surrogate model individual-wise, basically in the following steps: (i) An initial set E of individuals is collected, in which the considered empirical fitness η was evaluated (for example, the population of several first generations of the evolutionary algorithm). (ii) The surrogate model is constructed using the set of pairs {(x, η(x)) : x ∈ E}. (iii) The evolutionary algorithm is run with the fitness η replaced by the model for one generation with a population Q of size qP , where P is the desired population size for the optimization of η, and q is a prescribed ratio (e.g., q = 10 or q = 100). (iv) A subset P ⊂ Q of size P is selected so as to contain those individuals from Q that are most important according to the considered criteria for the progress of optimization.
(v) For x ∈ P, the empirical fitness is evaluated.
The set E is replaced by E ∪ P and the algorithm returns to the step (ii).
B. The generation-based-control consists in choosing between both kinds of evaluation generation-wise, basically in the following steps: The agreement between the results that are obtained with a surrogate model and those that would be obtained if the empirical objective function were evaluated depends on the accuracy of the model. A popular approach to increasing the accuracy of learning methods is boosting, i.e., construction of a strong learner through combining weak learners. It is important to realize that boosted surrogate models are only particular kinds of surrogate models and their interaction with optimization algorithms in optimization tasks follows the same rules as the interaction of surrogate models in general. In particular in the above outlines of individual-based and generation-based control, boosting is always performed in the step (ii), which has to be replaced with:
where | | denotes the cardinality of a set, ⌊ ⌋ the lower integer bound of a real number, and ⌈ ⌉ its upper integer bound.
(ii'b) For each j = 1, . . . , k, a surrogate model F j 1 is constructed, using only data not belonging to the j-th subset.
(ii'c) A k-fold crossvalidation of regression boosting is performed, and the error of the boosting approximation is in each iteration measured with the prescribed error measure on the validation data. (ii'd) The first iteration i in which the average error of the boosting approximation on the validation data is lower than in the i + 1-th iteration is taken as the final iteration of boosting. (ii'e) Boosting using the complete set {(x, η(x)) :
x ∈ E} is performed up to the final iteration found in step (ii'd), and the result of the application of the employed boosting method in each such iteration of boosting is taken as the boosted surrogate model in that iteration.
A particular method for MLP boosting has been presented in [11]. That method will now be employed in surrogate modelling with data from the investigation of catalytic materials for the high-temperature synthesis of hydrocyanic acid (HCN) [16]. The composition of most of those materials was designed by means of a specific genetic algorithm (GA) for heterogeneous catalysis [26]. As usually in evolutionary optimization of catalytic materials, the GA configuration was determined by the experimental conditions in which the optimization was performed: number of channels of the reactor in which the materials were tested, as well as time and financial resources available for those expensive tests. In the reported investigation, the algorithm was running for 7 generations of population size 92, and in addition 52 other catalysts with manually designed composition were investigated. Consequently, data about 696 catalytic materials were available. The considered MLPs had 14 input neurons: 4 of them coding catalyst support, the other 10 corresponding to the proportions of 10 metal additives forming the active shell, and 3 output neurons, corresponding to 3 kinds of catalytic activity considered as fitness functions.
For boosting, only data about catalysts from the 1.-6. generation of the GA and about the 52 catalysts with manually designed composition were used, thus altogether data about 604 catalytic materials. Data about catalysts from the 7. generation were completely excluded and left out for testing. The set of architectures to which boosting was applied was restricted to MLPS with 1 and 2 hidden layers and was delimited by means of the heuristic pyramidal condition: the number of neurons in a subsequent layer must not exceed the number of neurons in a previous layer. Let n I , n H and n O denote the numbers of input, hidden and output neurons, respectively, and n H1 and n H2 denote the numbers of neurons in the first and second hid-den layer, respectively. Then the pyramidal condition entails the following 90 architectures: (i) one hidden layer and 3 ≤ n H ≤ 14 (12 architectures); (ii) two hidden layers and 3 ≤ n H2 ≤ n H1 ≤ 14 (78 architectures). As was mentioned above, boosting can be combined both with the individaul-based and with the generation-based control of surrogate modelling. In the reported investigation of catalytic materials for HCN synthesis, the indiviual-based control was employed.
The error measure employed in the crossvalidation in the step (ii'c) was MSE. The distribution of the final iterations of boosting, found for MLPs with the 90 considered architectures in the step (ii'd), is depicted in Figure 2. We can see that only for 16 MLPs, already the 1st iteration was the final. For the remaining 74 MLPs, boosting improved the average MSE on the validation data for at least 1 iteration. The mean and median of the distribution of the final iterations were 6.6 and 5, respectively. For testing with the data from the 7th generation of the evolutionary algorithm, we used only the five MLPs most promising from the point of view of the average MSE on the validation data in the final iteration of boosting. These were the following MLPs:
a 1-hidden-layer MLP, with n H = 11 and the 3rd iteration of boosting being the final iteration, four 2-hidden-layers MLPs, with (n H1 , n H2 ) = = (10, 4), (10,6), (13,5), (14,8) and the final iterations of boosting 19, 32, 31 and 29, respectively.
For each of them, the validation proceeded as follows:
1. In each iteration up to the final, a single MLP was trained with data about all the 604 catalytic materials used for boosting.
2. In each iteration up the final iteration of boosting, the boosted surrogate model was constructed for the trained MLP, according to the step (ii'e). 3. From the values predicted by the boosted surrogate model for the 92 materials from the 7. generation of the GA, and from the measured values, the boosting MSE was calculated.
The results are summarized in Figure 3, decomposed to the properties corresponding to the MLP outputs -conversions of CH 4 and NH 3 and yield of HCN. They clearly confirm the usefulness of boosting for the five considered architectures. For each of them, boosting leads to an overall decrease of MSE of the conversion of CH 4 and HCN yield, on new data from the 7th generation of the GA, which is uninterrupted or nearly uninterrupted till the final boosting iteration. On the other hand, boosting did not lead to any decrease of the error of the conversion of NH3, which on the other hand is already from the beginning much lower than the two other performance measures (notice that the scale of the y-axis is 10-times finer for the conversion of NH 3 than for the conversion of CH 4 and HCN yield). The explanation for the different behavior of the conversion of NH 3 is the substantially lower variability of its values in the seventh generation of the GA, used for validating the usefulness of boosting (standard deviation, SD: 2.8, interquartile range, IQR: 1.6), compared to the conversion of CH 4 (SD: 26.1, IQR: 45.0) and HCN yield (SD: 20.1, IQR: 35.9). Due to so low variability, the conversion of NH 3 appears effectively as nearly constant during the validation of boosting, which in turn accounts for a nearly constant MSE. In addition, also implications and equivalences of important kinds of fuzzy logic are frequently extracted [8,15]. In general, extracted formulas of a formal logic are called rules. Over the last two decades, various rules extraction methods have been proposed for neural networks, but so far none of them has become a common standard (cf. the survey papers [1,15,22] and the monograph [7]). Here, a method for the extraction of Boolean implications from multilayer perceptrons with n inputs and m outputs will be sketched that finds to each output condition of the form:
one or more input conditions of the form C I : the value x of the input variables lies in a polyhedron P ⊂ R n (3)
Hence, this method extracts rules of the form:
A detailed explanation of the method can be found in [9]. Its main principles can be summarized as follows:
-An m-dimensional rectangular area R with borders perpendicular to the m coordinate axes has to be chosen in advance in the output space of a trained MLP with sigmoid activation functions. The reason for choosing such an area is that in the space of evaluations of m free variables, each m-dimensional rectangular area is the validity set of the conjunction of some m univariate Boolean predicates. That conjunction then serves as the consequent of the rule to extract. -The activation functions in the hidden neurons are approximated with piecewise-linear sigmoid activation functions. This can be done with an arbitrary precision. -The products of individual linearity intervals of all the activation functions determine areas in the input space in which the final approximating mapping computed by the multilayer perceptron is linear.
-In each such area, all points mapped to R form a polyhedron, which may eventually be empty or may be concatenated with polyhedra from some of the neighboring areas to a larger polyhedron. -The union of all the nonempty concatenated polyhedra P 1 , . . . , P q defines the antecedent of a rule in a combined form
which is equivalent to a logical disjunction of q rules of the simple form (4):
To increase the comprehensibility of the extracted rules, visualization by means of 2-or 3-dimensional cuts of the set P 1 ∪ • • • ∪ P q can be used (Figure 4).
Usually, logical rules of the form (4) are the final results of this rule-extraction method. Nonetheless, there is one exception -when the polyhedron P is also rectangular with borders perpendicular to axes, or more generally, when P can be approximately replaced with such a rectangular area R I in the input space. Then the above rule (4) can be approximately expressed in the conjunctive form
Here, I 1 , . . . , I nI are intervals that constitute the projections of R I into the n I input dimensions. Each such interval can be restricted both from below and from above, restricted only from below or only from above, or finally can be even the complete set of real numbers. However, dimensions for which the corresponding projection of R I equals the complete real axis are usually not included in (7), since they would not provide any new knowledge. Finally, observe that due to ( 5) and ( 7), the final extracted rule is in the disjunctive normal form.
In the rule-extraction method outlined above, the possibility of replacing a polyhedron P with a rectangular area R I is assessed according to the following principles:
1. The resulting dissatisfaction with points that either belong to P but do not belong to R I , or belong to R I but do not belong to P (i.e., with points from the symmetric difference R I ∆P ), has to remain within a prescribed tolerance ε and R I has to be minimal in the input space among rectangular areas of some specified kind with dissatisfacion within that tolerance. 2. The dissatisfaction with points from R I ∆P depends solely on those points and is increasing with respect to inclusion. Consequently, it can be measured using some monotone measure on the input space, possibly depending on P . 3. To be eligible for replacement, P has to cover at least one point of the available data.
For 2., the most attractive monotone measures, due to their straightforward interpretability, are:
-The joint empirical distribution of the input variables in the available data.
-The conditional empirical distribution of the input variables in the available data, conditioned by P .
Rules of the form (7) are also very convenient from the visualization point of view: Since cuts of rectangular areas coincide with the corresponding projections of those areas, the values of no variables need to be fixed.
As an example, Figure 5 shows three-dimensional cuts determining the antecedents of conjunctive-form rules extracted from a trained MLP with 5 input neurons and 1 output neuron such that: (i) the input neurons correspond to variables that record the molar proportions of the oxides of Fe, Ga, Mg, Mn and Mo in the catalytic material; (ii) the output neuron corresponds to a variable recording propene yield. The extracted rules are listed in Table 1.
The paper dealt with employing feed-forward neural networks for knowledge discovery from data about chemical materials. It has shown that in this application area, obtaining numeric knowledge by neural-network regression is justified, in spite of the fact that numeric knowledge is substantially less human-understandable Table 1. Antecedents of the rules of the form (7) extracted using the method described in this section for the consequent "propene yield > 8%" from a trained MLP with 5 input neurons and 1 output neuron, assuming that the above interpretation of the variables to which those neurons correspond is described by (i) and (ii).
Fig. 5. A three-dimensional projection of the union of rectangular areas that replace, following the method described in this section, the union of of polyhedra from the antecedent of a combined form rule extracted from a trained MLP. The projection corresponds to input variables recording the molar proportion of oxides of Ga, Mg and Mo in a catalytic material. The numbers 1, 2, 3 refer to the antecedents of the rules in Table 1.
than symbolic knowledge. Its justification consists in the possibility to use such knowledge in the optimization tasks entailed by search for new materials in the surrogate modelling approach.
In addition to justifying the specific need for numeric knowledge from neural network regression in this application area, the paper recalled the possibility to obtain symbolic knowledge in the form of logical rules from trained neural networks. It explained a recently proposed method for the extraction of Boolean rules in disjunctive normal form, and illustrated it on data about catalytic materials.