George Potamias
email: potamias@ics.forth.gr

Sofia Kaforou

Dimitris Kafetzopoulos

Towards Optimal Microarray Universal Reference Sample Designs: An In-Silico Optimization Approach

Keywords: Bioinformatics, microarrays, universal reference sample

Assessment of the reliability of microarray experiments as well as their cross-laboratory/platform reproducibility rise as the major need. A critical challenge concerns the design of optimal Universal Reference RNA (URR) samples in order to maximize detectable spots in two-color/channel microarray experiments, decrease the variability of microarray data, and finally ease the comparison between heterogeneous microarray datasets. Towards this target we devised and present an in-silico (binary) optimization process the solutions of which present optimal URR sample designs. Setting a cut-off threshold value over which a gene is considered as detectably expressed enables the process. Experimental results are quite encouraging and the related discussion highlights the suitability and flexibility of the approach.

Introduction

Scientific experimental science is founded on the 'right' design of experiments where some common-sense principles should apply [START_REF]Design of Microarray Experiments[END_REF]. The observation is more than critical in the set-up and the complex design of microarray experiments. With more than fifteen years of research and experimentation in the field the real challenge as well as the actual tendency is moving towards the comparison of different microarray experiments originated from different platforms, from different laboratories, and from different designs. In this context, assessing the reliability of microarray experiments as well as their cross-laboratory/platform reproducibility rise as the major need.

Measuring reliable fluorescence intensities with microarrays is not a straightforward task. The basic problem rises from the variability in microarray spot geometry and the quantity of DNA deposited at each spot. So, absolute fluorescence intensity cannot be used as a reliable measure of RNA level. However, if two RNA samples are differentially labeled and co-hybridized to spots on the same microarray, the ratio of their signal intensities accurately reports the relative quantity of RNA targets in both samples -the so called, two-color (two-channel) hybridization set-up (2CH). Obtaining reliable and reproducible 2CH gene expression data is critically important for understanding the biological significance of perturbations made on a cellular system. Moreover, and most critical, maximizing detectable spots on the reference image channel also decreases the variability of microarray data allowing for reliable detection of smaller differential gene expression changes. [START_REF] Khan | A universal reference sample derived from clone vector for improved detection of differential gene expression[END_REF].

With the 2CH microarray technology there are two alternatives for the comparison between different sample conditions, C 1 vs. C 2 : (i) the 'loop' design where, direct comparisons between a series of arrays are performed with a sample of type C 1 in one channel and one of type C 2 in the other channel for all other arrays, and the comparison of samples in circular or multiple pair-wise fashion -it is a useful approach when a small number of samples is to be analyzed; and (ii) the Universal RNA Reference (URR) sample, C -devised by a pool of diverse (in terms of their ability to express different families of genes) cell-lines. With the URR approach a series of hybridizations are carried out with an experimental sample of type C 1 / C 2 in one channel and the URR sample in the other. Then differences in gene-expression between C 1 and C 2 are acquired by comparing the ratios C 1 /C vs. C 2 /C [START_REF] Manduchi | Issues Related to Experimental Design and Normalization[END_REF]. The hybridization of each (experimental) sample with a URR mixture serves as a common denominator between different microarray hybridizations [START_REF] Eisen | DNA arrays for analysis of gene expression[END_REF]. With the employment of a reproducible URR the quantization of gene-expression levels is more reliable and offers a valuable tool for monitoring and controlling intra-and inter-experimental variation. The URR design has several practical advantages over the loop design: (i) it extends easily to other experiments if the common URR is preserved, (ii) is robust to multiple chip failures; and (iii) reduces incidence of laboratory mistakes as each (experimental) sample is handled the same way [START_REF]Design of Microarray Experiments[END_REF]. In addition, the URR design facilitates the normalization and the subsequent comparison between heterogeneous microarray data sets [START_REF] Novoradovskaya | Universal Reference RNA as a standard for microarray experiments[END_REF]. Typically molecular scientists and experimentalists prepare their own array-specific reference samples e.g., genomic DNA [START_REF] Williams | Genomic DNA as a cohybridization standard for mammalian microarray measurements[END_REF][START_REF] Gadgil | An analysis of the use of genomic DNA as a universal reference in two channel DNA microarrays[END_REF], a mixture of clones spotted on the arrays [START_REF] Sterrenburg | A common reference for cDNA microarray hybridizations[END_REF], as well as complementary to every microarray spot short oligomers [START_REF] Dudley | Measuring absolute expression with microarrays with a calibrated reference sample and an extended signal intensity range[END_REF]. Many groups utilize a mixture of cell-lines for their URR (e.g., Stratagene's human URR [START_REF] Novoradovskaya | Universal Reference RNA as a standard for microarray experiments[END_REF]) or, clone-vectors [START_REF] Khan | A universal reference sample derived from clone vector for improved detection of differential gene expression[END_REF]. In most of the cases the utilized URR samples are not reproducible between labs and may provide detectable signal for a low percentage of the microarray spots [START_REF] Yang | Within the fold: assessing differential expression measures and reproducibility in microarray assays[END_REF]. Subsequently, they do not provide good expression coverage, i.e., a high proportion of genes being adequately hybridized and expressed in order to avoid spots with zero denominators in C 1 /C vs. C 2 /C ratio comparisons, a fact that would force discarding those spots from the analyses [START_REF] Manduchi | Issues Related to Experimental Design and Normalization[END_REF].

It is natural to assume, and this is actual the case, that as more cell-lines are utilized to prepare the URR sample mixture then, the biggest the chance to get dilution and saturation effects. So, it is of fundamental importance to determine an optimum number of cell-lines so that such effects are avoided as much as possible. This is the target of the present work. In particular, our aim was to devise a general methodology that guides the determination of an optimal ("the less the better") number of cell-lines the mixture of which is capable to hybridize and express the highest possible number of array genes. The paper is organized as follows: section 2 outlines the state-of-theart research on the field; in section 3 we introduce and formally present our optimization approach, and in section 4 we present and discuss the results of the performed experiments; finally, in section 6 we conclude and point to future research plans.

Background to the URR sample design

The MAQC (MicroArray Quality Control) project2 systematically evaluated the reliability of microarray technology by profiling two human URR samples using different microarray and QPCR (Quantitative PCR) platforms [START_REF] Shi | The MicroArray Quality Control (MAQC) project showsinter-andintraplatformreproducibilityofgeneexpressionmeasurements[END_REF]: Stratagene's (acquired by Agilent Technologies3 ) human Universal Reference RNA™ -derived from a mixture of cell lines [START_REF] Novoradovskaya | Universal Reference RNA as a standard for microarray experiments[END_REF], and Ambion's human brain URR -pooled from multiple donors and several brain regions 4 .

According to Sterrenburg et al., a common reference for DNA microarrays would consist of a mix of the products being spotted on the array. Their polymerase chain reaction (PCR) reference was generated by pooling a fraction of all amplified microarray probes prior to printing. A very low number of spots could not be analyzed because reference, as opposed to target, did not give a significant signal [START_REF] Sterrenburg | A common reference for cDNA microarray hybridizations[END_REF]. Yang et al., suggested the use of a limited number of cell lines, each expressing a large number of diverse genes. They constructed two reference pools from those cell lines with the greatest representation of unique genes, which are also easy to grow and yield high quantities of RNA. They found that adding more cell lines to the pool would not necessarily improve the overall gene representation because some genes were diluted below the detection limit. The first reference sample exhibited similar coverage with Statagene's human URR sample (75%), and the second reference sample had coverage equal to 80%. Thus, according to Yang et al., a simple pool of RNA from diverse cell lines can provide a superior reference [START_REF] Yang | Within the fold: assessing differential expression measures and reproducibility in microarray assays[END_REF]. Human, mouse and rat reference RNA samples were considered by Novoradovskaya et al. [START_REF] Novoradovskaya | Universal Reference RNA as a standard for microarray experiments[END_REF] and were prepared from pools of RNA derived from individual cell lines representing different cell-lines (ten human, eleven mouse and fourteen rat). They evaluated microarray coverage based on a pre-specified threshold equal to the background intensity or twice the background intensity of each channel, using different microarray platforms. Probes with intensities above threshold were characterized as present. They reported microarray coverage greater than 80% for all arrays tested when threshold was equal to background, and greater than 60% when threshold was equal to twice the background. Consequently, they agreed with Yang et al., that pools of RNA derived from a limited but diverse set of cell lines result in an optimal reference sample.

Furthermore, in a recent study we examined whether a small number of cell-lines (in a total of 22), each of which expresses different genes, could outperform more complex reference mixtures [START_REF] Tsiliki | A computational approach to microarray universal reference sample[END_REF]. Following different techniques (exhaustive combination of cell-lines, heuristic search and stochastic simulated annealing) we achieved an optimal cell-line mixture of 11 cell-lines with theoretical gene coverage of 57.84%. The mixture was tested with a wet-lab experiment (i.e., it was hybridized with the mixture of all 22 cell-lines), and showed coverage of 57.12% -quite near to the theoretical one. Testing (with a similar experimental set-up) the standard Stratagene URR mixture an inferior coverage of 54.02% could be reached.

The relatively low gene coverage (around 57%) but also, the need to offer a more natural, universal and robust design approach forced us to consider a straightforward optimization approach. In particular, we utilized the mathematical programming framework and devised a linear/binary programming optimization methodology. Our methodology is flexible enough and could be easily tuned to fit the specifics of different microarray platforms and heterogeneous experimental settings.

Towards Optimal Reference Sample Designs

Assume a set of genes, G, and a set of cell-lines, C. In addition, we assume that the URR sample is to be devised using equal shares of fixed quantity from each cell-line. A series of C two-channel microarray experiments are performed, each time hybridizing an individual cell-line with the reference-mix. After standard, crossexperiment normalization of the acquired intensities, we get the gene-expression profile for each cell-line, and the respective gene-expression matrix is formed -a matrix with G rows and C columns. Our fundamental concern is twofold: which genes are expressed over all cell-lines (i.e., across rows), and which genes are expressed by each cell-line (i.e., across columns). We address the problem by setting a cut-off value, with gene-expression values greater than the cut-off to be considered as detectably expressed (with the remaining considered as not detectable). In most microarray experiments expression cut-off values could be acquired either during the microarray image analysis phase where, the (manual) inspection of the underlying background intensities guides to the determination of a specific cut-off value in order to consider the gene signal as detectable or, by the inspection of the normalized geneexpression matrix and the calculation of different cut-off values for different foldchange gene expression rates. As the determination of the cut-off value is beyond the scope of the current work, we assume an absolute cut-off based on the median over all gene-expression values.

We consider a gene as detectably expressed if an adequate number of selected celllines express it, i.e., the corresponding gene's expression values for these cell-lines are over the cut-off. In the present paper we report results with two different cut-off values: one that equals the median, and one being set to the 70% of the median of all gene-expression values. Elaborating on the previous discussion, we post two criteria in order to consider the number of cell-lines (to be utilized in the URR) as optimal: (i) each array gene should exhibit an adequate number of its corresponding expression-values over the cut-off threshold -in this case we consider the gene as expressed, and (ii) the number of expressed genes is over a (user-defined) pre-specified percentage of all array genes. All genes considered as expressed by (i) are considered as covered, and their percentage over all the array genes presents the coverage figure of the final URR mixture.

The Optimization Methodology

In order to meet and fulfill the aforementioned criteria we devised an optimization process enabled by a linear/binary programming problem formulation. Notation. Cell-line: given a set of m cell-lines, assume, C j , a set of binary variables, C j ∈ {0,1}, 0 ≤ j ≤ m to denote the presence (C j = 1) or, absence (C j = 0) status of a cell-line. We do not cope with the case of multiple (i.e., integer or, real) values for C j . This is a more difficult problem that calls for a more elaborate and complex modeling. In such a case we need more microarray experiments in order to assess exact dilution and saturation effects when the cell-lines are taking part in a reference-mix in nonequal shares; Gene-Expression value: given a set of n genes, G i , 0 ≤ i ≤ n, we denote with G ij the expression value of gene i for cell-line j; Cut-Off: the pre-specified cut-off threshold is denoted with T. Based on this notation we form the following Binary Programming Optimization Problem (BPOP): minimize OF:

€ C j=1 m ∑ j (1)
subject to:

€ C j= 1 m ∑ j * Gij ≥ kT (2) i , 1 ≤ i ≤ n C ∈ {0,1}, G ij ∈ R 0 ≤ k ≤ m
Solutions to the BPOP problem will assign values to the binary variables C j . If C j =1 then, the respective cell-line is to be utilized in the URR sample, otherwise it will not. The objective function (1) seeks for the minimum number of cell-lines to be utilized. Each constraint in (2) states the following requirement: the weighted sum of geneexpression values (over all present cell-lines) should be greater or equal to kT, i.e., k times the cut-off value T. In other words, we try minimize the number of cell-lines, and at the same time to keep the gene-expression values over k times the (prespecified) cut-off threshold. In this way, the constraints guide the minimization process to select those cell-lines for which at least one or more of its gene-expression values is over the cut-off. In the ideal, the corresponding gene-expression values will be greater than the cut-off in all present cell-lines -of course this may happen if all genes exhibit at least one of their expression values (across all present cell-lines) over the cut-off. What we really seek is an adequate number of cell-lines for which each gene exhibits an expression value over the cut-off. Tolerating (by varying the k control number) the required number of cell-lines covering a gene the BPOP process results into different solutions. For example, we may consider a gene as expressed if at-least-one of its expression values is over the cut-off, i.e., at-least-one cell-line covers it. This requirement is guaranteed by the formulation of BPOP -the constraints presented by (2) guides the optimization process to select those cell-lines for which, not only just one value is over the cut-off but also all other values are high enough in order to pass the kT threshold. So, for different and increasing values of k the actual BPOP solutions meet the aforementioned criteria for the design of an optimal reference sample.

Experiments

A set of 22 tissue-specific cell-lines were targeted as candidates for the design of the optimal URR (their code, name and tissue-origin are shown in Table 1). With these cell-lines a set of 22 microarray experiments were conducted on a custom-made microarray platform comprising 34771 gene transcripts (for details of cell-line hybridizations, acquisition of signal-intensities and normalization across experiments the reader may refer to [START_REF] Tsiliki | A computational approach to microarray universal reference sample[END_REF]). As a number of missing gene-expression values were present, and in order to avoid as much as possible noise effects, we conducted a data cleaning procedure: (a) genes having more than half, across all cell-lines, missing values were eliminated, and (b) for the remaining genes we replaced missing values with the average of the respective gene's expression values. The data cleaning process concluded into a final set of 34673 genes. We set the cut-off threshold to be equal to the median (T=368.64) of all geneexpression values, and to a proportion of 70% of the median (T = 258.05). Because of space limitations we do not report the full spectrum of results but just some observations should be reported: (a) a percentage of 28.1% and 5.0% of genes are not covered by any cell-line, (b) the percentage of genes being covered by just-one (k=1) time the respective cut-off threshold is 7.4% and 4.5%, respectively -in this case the respective genes are considered as tissue-specific, and are not taking part in the optimization process, and (c) the percentage of genes being covered by all, (k=22) times the respective cut-off threshold is 31.5% and 46.7%, respectively -these genes are considered as tissue-independent, and again are not taking part in the optimization process (whatever cell-line is to be included in the final URR, these genes will always exhibit gene-expression levels over the respective cut-off threshold; and (iv) the percentage of genes being covered by at-least-one (k≥1) times the respective cut-off threshold is 71.9% and 95.0%, respectively -this observation validates the determination of the median as a rational cut-off value.

Results and Discussion

We run the BPOP process for each of the specified cut-off thresholds using the MOSEK optimization suite 5 . MOSEK optimization software solves large-scale mathematical optimization problems, and provides specialized solvers for linear/binary programming and mixed integer programming problems. In the sequel we present and discuss on the respective results. T = 368.64 cut-off. Varying the k control number from 2 to 21 (times the cut-off threshold value -recall that genes with just one, k=1, and all, k=22 of their expression values over the cut-off are discarded), we came up with the results presented in Table 2. The 'k' row refers to the number of times the sum over all weighted geneexpression values exceeds the cut-off -it controls the lower limit for the formed BPOP constraints. Note that for k>13 no solutions could be found, i.e., there are no genes covered by more than 13 cell-lines. The respective non-zero binary variables in each solution are shown in row 'C j =1', for example, when k=4 the solution includes 7 non-zero, i.e., the sum of genes' expression values is greater or equal to four times the cut-off is satisfied. In other words, all constraints are satisfied at an optimal selection of seven cell-lines. Setting a lower (user specified) coverage limit of 60%, and in the case that all 22 cell-lines are to be included in the final URR mixture, all genes should exhibit gene-expression levels over the cut-off in at-least-three cell-lines. For k=3 the corresponding adequate number of cell-lines is equal to at-least-one. On this level of adequacy, i.e., when we want the genes to exhibit in at-least-one cell-line an expression value over the cut-off, we achieve coverage of 61.9% with the inclusion of just 6 cell-lines. For k=4 the corresponding adequate number of cell-lines is still to at-least-one. On this level of adequacy we achieve coverage of 63.3% with the inclusion of just 7 cell-lines. For k=8 the corresponding adequate number of cell-lines is equal to at-least-two. On this level of adequacy, i.e., when we want the genes to exhibit in at-least-two cell-lines expression values over the cut-off, we achieve coverage of 60.5% in the presence of 13 cell-lines. When the corresponding adequate number of cell-lines is equal to atleast-three, we may achieve a coverage of 60.4% but with a costly reference sample design of 21 cell-lines being present (k=13). With a careful inspection of the results we devised a map of the cell-lines being present in the aforementioned solutions (Table 3).

Table 3. Present cell-lines (marked) in the solutions for the design of an optimal reference sample for cut-off T=368.64

From table 3 we may observe that there are exactly four cell-lines, HL60, Panc1, JAR and HepG2, being present in all solutions. So, it is natural to consider them as the most promising and the only components of the optimal URR sample mixture. With just these cell-lines in the final URR reference mixture we were able to achieve coverage of 60%. Note, that two of these cell-lines (Panc1, HepG2) are also included in the experimentally tested URR mixture that achieves coverage of 57.12% (refer to [START_REF] Tsiliki | A computational approach to microarray universal reference sample[END_REF]). T=258.05 cut-off. The BPOP process was also applied on the same data with a lower cut-off threshold T = 258.05 (70% of the median) but with a higher coverage limit of 75%. Because of space limitations we do not report results analogue to Tables 2 and3. With a careful inspection of the results we again came up with a map of the cell-lines being present in the solutions that meet the aforementioned criteria. There are thirteen of them with just six cell-lines, HL60, Hs578T, Raji, HepG2, LNCap, and WERI1 being present in at least eleven of the thirteen solutions. So (as in the case of the 368.64 threshold) it is natural to consider them as the most promising and the only components of the optimal URR sample mixture. With just these six cell-lines we were able to achieve a coverage of 87.8%. Again, five (Hs578T, Raji, HepG2, LNCap, and WERI1) of the six cell-lines are included in the experimentally tested URR mixture (refer to [START_REF] Tsiliki | A computational approach to microarray universal reference sample[END_REF]). So, and provided that with the presence of less cell-lines saturation and dilution effects are avoided as much as possible, we may hypothesize that the presented optimization approach greatly improves previous URR sample designs.

Concluding Remarks

We have presented an in-silico methodology for the design of optimal URR samples suited for 2CH microarray experiments. The biological problem to overcome relates to the avoidance (as much as possible) of dilution and saturation effects taking place when a big number of cell-lines are present. In the final URR mixture sample. We approached the problem with the careful formation of a binary programming optimization methodology (BPOP). In particular, we tried to minimize a set of binary variables that represent the inclusion or not of a cell-line in the final URR mixture, and maximize the overall expression of genes. Setting different cut-off values, to determine if a gene is detectably expressed or not, we achieved high enough gene coverage percentages. In particular, setting the cut-off to the 70% the median of gene expression values, we were able to achieve a quite high coverage of 87.8% with the presence of just six cell-lines. A figure that outperforms current state-of-the-art URR sample designs.

Our future R&D plans include: (a) testing and assessment of the introduced optimization methodology on other cell-line configurations, and especially on integrated microarray data from different platforms and studies, and (b) elaboration of more advanced optimization techniques, e.g., integer and/or real programming optimization in order to extend our approach to varying cell-line shares in the final URR mixture.

Table 1 .

 1 The target cell-lines

	Cell line	Tissue origin
	1. HL60	Bone marrow
	2. Hs578T	Mammary Gland
	3. McF7	Mammary Gland
	4. OVCAR3	Ovary
	5. Panc1	Pancreas
	6. SKMEL3	Skin
	7. SKMM2	Bone Marrow
	8. T47D	Mammary Gland
	9. TERA1	Testis
	10. U87MG	Brain
	11. Raji	B-lymphoblasts
	12. JAR	Genital
	13. Saos2	Bone
	14. SW872	Liposarcoma
	15. THP1	Peripheral Blood
	16. HCT116	Colon
	17. HUVEC	Umbelical Vein
	18. HepG2	Liver
	19. HeLa	Cervix
	20. LNCap	Prostate
	21. Molt4	T-lymphoblasts
	22. WERI1	Retina

Table 2 .

 2 Gene coverage results of the BPOP process for cut-off T=368.64

http://www.fda.gov/ScienceResearch/BioinformaticsTools/MicroarrayQualityControlProject/default.htm

http://www.genomics.agilent.com/

http://www.ambion.com/catalog/CatNum.php?6050

http://www.mosek.com


Behrouz Minaei-Bidgoli

Maryam Asadi
email: asadi@iust.ac.ir

Hamid Parvin
email: parvin@iust.ac.ir

An Ensemble Based Approach for Feature Selection

Keywords: Feature Selection, Ensemble Methods, Fuzzy Entropy

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

We have to use features of a dataset to classify data points in pattern recognition and data mining. Some datasets have a large number of features. Processing these datasets is not possible or is very difficult. To solve this problem, the dimensionalities of these datasets should be reduced. To do this, some of the redundant or irrelevant features should be eliminated. By eliminating the redundant and irrelevant features, the classification performance over them will be improved. Three different approaches are available for feature selection mechanism [START_REF] Tan | Introduction to Data Mining[END_REF]. The first ones are embedded approaches. In these algorithms, feature selection is done as a part of the data algorithm. The second ones are filter approaches. These algorithm selected features before the data mining algorithm is run. The last ones are wrapper approaches. In these algorithms the target data mining algorithm is used to get the best subset of features.

A lot of methods for feature subset selection have been presented, such as similarity measures [START_REF] Tsang | OFFSS: optimal fuzzyvalued feature subset selection[END_REF], gainentropies [START_REF] Caruana | Greedy attribute selection[END_REF], the relevance of features [START_REF] Baim | A method for attribute selection in inductive learning systems[END_REF], the overall feature evaluation index (OFEI) [START_REF] De | Neuro-fuzzy feature evaluation with theoretical analysis[END_REF], the feature quality index (FQI) [START_REF] De | Neuro-fuzzy feature evaluation with theoretical analysis[END_REF], the mutual information-based feature selector (MIFS) [START_REF] Battiti | Using mutual information for selecting features in supervised neural net learning[END_REF], classifiability measures [START_REF] Dong | Feature subset selection using a new definition of classi-fiability[END_REF], neuro-fuzzy approaches [START_REF] De | Feature analysis: neural network and fuzzy set theoretic approaches[END_REF][START_REF] Platt | Using analytic QP and sparseness to speed training of support vector ma-chines[END_REF], fuzzy entropy measures [START_REF] Shie | Feature subset selection based on fuzzy entropy measures for handling classification problems[END_REF], etc.

In this paper we try to improve Shie-and-Chen's method. We try to solve the drawback of parameter sensitivity. To do this we use ensemble method. We get the results for different threshold values. For each threshold values, we get a subset of features. We give a score to each feature in these subsets. Finally by use of ensemble concept, we select the features which have the highest scores. This method is not a parameter sensitive one, and also it has been shown that using the method based on the fuzzy entropy results in more reliable selected features than the previous methods'.

Proposed Algorithm

Shie-and-Chen's Algorithm which is presented in [START_REF] Shie | Feature subset selection based on fuzzy entropy measures for handling classification problems[END_REF] is parameter sensitive. So if these parameters change, the result of algorithm can be changed significantly. When these parameters are given by the user, the quality of algorithm results will be even weaker. Because user selects the parameters randomly and experimentally, so it is possible that they are not proper values for an exemplary dataset. So the result of algorithm is not trustable. Also the proper values are not available for some datasets which are not used in this algorithm. So to find the best result we need to test the algorithm for a lot of possible threshold values. Then we must select the threshold values which cause the best results. To solve this problem we use ensemble method.

We do not select threshold values experimentally in our algorithm. Our algorithm test different possible values for thresholds and then by doing some steps, it selects the subset of features. This algorithm has 5 steps. We employ Shie-and-Chen's method by a little change in our algorithm. The result of their algorithm is a subset of features. But we get a sequence of features instead of a subset. Actually the order of feature appearance is important in our algorithm. First step runs Shie-and-Chen's method for each pair of (T r , T c ). The result of algorithm at this step is a table of feature sequences which are selected for each pair of threshold values. For example the result of our algorithm for Iris is shown in Table 1. We obtained this result for 5 different values for T c and T r . Each element in this table is a feature sequence selected by the algorithm of Fig. 1 with a different pair of threshold values.

For T r = base_t r : step_t r :1 For T c = base_t c : step_t c :1 AllFSeq(T r , T c ) = Shie-and-Chen's algorithm(T r , T c );

Fig. 1. Pseudo code of the first step of algorithm

It has two loops. One of them slides over T r and the other one slides over T c . Two parameters base_t r and base_t c are the minimum values used for T r and T c respectively. Two parameters step_t r and step_t c determine the distance between two consecutive threshold values of parameters T r and T c respectively. FSeq is a two dimensional matrix whose elements are features sequences obtained by the algorithm of Fig. 2 with each possible tested pair of threshold values. As it is inferred from Table 1, at the first and the last rows of each column we have some similar results for some threshold values. There is a similar discussion about the first and the last columns of each row. The results of algorithm for the first and the last columns of each row and the first and the last rows of each column are not trustable to reach some proper threshold values. Since these results have strongly negative effect on the final evaluation, at the second step we have to remove these repetitions. This step has two parts. The first part removes the repetitions of columns and the second part removes the repetitions of rows. First part keeps only the results at the beginning and ending of each column to reach a dissimilar result at the beginning and ending of each column. And the second part keeps only the results at the beginning and ending of a row to reach a dissimilar result at the beginning or ending of each row. In other words, we use only one of the same results at the beginning and ending parts of each row and each column in final evaluation. The following pseudo code is the first part of second step of the algorithm.

New_AllFSeq = AllFSeq For T r = base_t r : step_t r :1 q = base_ t c While (true) q = q + step_ t c if is_same ( AllFSeq ( T r , base_t c ) , AllFSeq ( T r , q )) New_AllFseq ( Tr , q ) = EmptySeq else break q = last_ t c While (true) q = q -step_ t c ; if is_same ( AllFSeq ( T r , last_t c ) , AllFSeq ( T r , q )) New_AllFseq ( T r , q ) = EmptySeq else break 1is a function that checks the similarity of its inputs. It has two input parameters which can be two sequences of features. If they are similar, the output will be 1 and if they are not similar the output is 0.

݁݉ܽݏ_ݏ݅ ሺܽ , ܾሻ ൌ ൜ 1, ݂݅ ݔ ൌ ݕ 0, ݁ݏ݅ݓݎ݄݁ݐ 1
It checks the similarity between the first sequence of a column and the consecutive sequences of that column. By reaching the first dissimilar sequence at the beginning or ending of a column, this part of algorithm is done for each column. Output for Iris example of doing the first part of the second step of the algorithm is available in Table 2 by horizontal shading (+ sings). The second part of the second step of the algorithm is like Fig. 2. It is like the first part of the second step. It checks the similarity between the first sequence of a row and the other sequences in that column. By reaching the first dissimilar sequence at the beginning or ending of a row, this part of algorithm is done for each row. Result of doing the second part of the second step of the algorithm over the Iris dataset which is obtained from the first step is shown in Table 2 by vertical shading (* sings). Third step uses majority voting to reach the best subset of features. We have to give a score to each feature. There is a subset of selected features for each pair of T r and T c . We change this subset to a sequence of features by their ranks of appearing at the first step. In other words each feature that appears sooner has more effect on output, so it is given a higher score. Then we sum all given scores to features for each pair of threshold values. We define the score of each feature as equation 2.

After obtaining Table 2 for each dataset, we give a score to each of its features according to equation 2. In the equation 2, we give the higher weight to the first feature which appears sooner, and we give the lower weight to the last feature which appears at the end of the sequence. For example if there are 10 features, the weight of the first feature is considered 10, and the weight of the last feature is considered 1.

݁ݎܿܵ ൌ ∑ ∑ ∑ ݈ܽሺݑݍ݁ݏ݅ ݍ݁ܵܨ݈݈ܣ ሺ ܶ , ܶ ሻሺ݅ሻ, ݂ሻ * ሺ |ݍ݁ܵܨ݈݈ܣ| െ ݅ 1ሻ ெ௫ௌி ୀଵ ் ் 2
where MaxSF is obtained by equation 3. ܨܵݔܽܯ ൌ ݔܽ݉ ் ,் , ݍ݁ܵܨ݈݈ܣ|ሺ ሺ ܶ , ܶ ሻሺ݅ሻ|ሻ 3 Finally we sum all the weighted scores obtained by the algorithm for different pairs of threshold values. For example, in the Iris example the MaxFS is 3. In the example we get these results: Score (3) = 21, Score (4) = 21, Score (1) = 7 and Score (2) = 2. Table 3. Comparison between feature subsets selected by our and Shie-and-Chen's methods.

Feature subsets selected by two methods

Data sets

Our method Shie-and-Chen's method {4,3} Then we sort all features by their scores. After that we select the features with maximum scores. We select the same number of features as the Shie-and-Chen's method. In Iris example the subset of {3, 4} features is selected as final selected sub-set, because these features have the highest scores, and Shie-and-Chen's method selected two features for this example.

Experimental Results

We have implemented our feature selection algorithm in Matlab. We use Weka software to evaluate the mapped datasets into the selected features obtained by our feature selection algorithms. We compare the feature subsets selected by our method with those selected by Shie-and-Chen's method in Table 3 for all of datasets which are used to compare in [START_REF] Shie | Feature subset selection based on fuzzy entropy measures for handling classification problems[END_REF]. Also Table 4 shows that the obtained accuracies of different classifiers on the selected features obtained by proposed method are better than the obtained accuracies of the same classifiers on the selected features obtained by Shie-and-Chen's algorithms the most datasets.

Conclusion

This paper improves one of the existing feature selection algorithms, Shie-and-Chen's method. This feature selection algorithm uses fuzzy entropy concept. The problem of Shie-and-Chen's method is that it is a parameter sensitive algorithm. User should select threshold values in that algorithm experimentally. The result of algorithm for some threshold values is very weak and it is not trustable. To solve this problem we use ensemble method. Our paper runs Shie-and-Chen's algorithm for different values as thresholds and then gives a weight to each selected features according its rank. Finally by using one of the ensemble methods, majority voting, it selects the best features which have the highest scores. So this algorithm does not need any input parameter. Also the obtained accuracies of different classifiers on the selected features obtained by proposed method are better that the obtained accuracies of the same classifiers on the selected features obtained by Shie-and-Chen's algorithms.

Fig. 2 .

 2 Fig. 2. Pseudo code of the first part of the second step of the algorithm

Table 1 .

 1 Feature subsets selected for some pairs of threshold values over Iris dataset.

	Tc, Tr	0.01	0.21	0.41	0.61	0.81
	0.01	4, 3	3, 4	3, 4	3, 4	3, 4
	0.21	4, 3	4, 3	3, 4	3, 4	3, 4
	0.41	3, 4	4, 3	3, 4	3, 4	3, 4
	0.61	4, 1	4, 2	3, 1	3, 1	3, 1
	0.81	3, 1	4, 3, 1	3, 4	3, 4	3, 4
	Equation

Table 2 .

 2 Delete repetitions in columns ofTable 1 then delete repetitions in rows of Table 1.

	Tc, Tr	0.01	0.21	0.41	0.61	0.81
	0.01	4, 3	*	*	*	3, 4
	0.21	+	+*	+*	+*	+
	0.41	3, 4	4, 3	+*	+*	+
	0.61	4, 1	4, 2	*	*	3, 1
	0.81	3, 1	4, 3, 1	*	*	3, 4

Table 4 .

 4 Comparing classification accuracy rates of our and Shie-and-Chen's methods

Average classification accuracy rates of different me- thods Classifiers Data sets Shie-and-Chen's me- thod

		C4.5	67.14± 0.61 %	67.14 ± 0.61%
		LMT	97.22 ± 0.47%	99.77 ± 0.10%
	Monk-3 data set	Naive Bayes SMO	97.21 ± 2.71% 97.22 ± 0.47%	97.21 ± 2.71% 100.00 ± 0.00%
		C4.5	100.00 ± 0.00%	100.00 ± 0.00%
			Our method
		LMT	76.30 ±4.84%	77.22 ± 4.52%
	Pima diabetes data set	Naive Bayes SMO	76.30 ±4.84% 75.65 ±5.61%	77.47 ± 4.93% 77.08 ± 5.06%
		C4.5	94.62 ±2.12%	74.88 ± 5.89%
		LMT	82.42 ± 5.34%	82.87 ± 6.23%
	Cleve data set	Naive Bayes	80.41 ± 3.95%	84.48 ± 3.93%
		SMO	80.00 ± 5.99%	83.51 ± 6.09%
		C4.5	76.90 ± 8.40%	76.90 ± 8.40%
		LMT	100.00 ± 0.00%	100.00 ± 0.00%
	Correlated data set	Naive Bayes SMO	86.03 ± 3.75% 89.87 ± 6.88%	86.03 ± 3.75% 89.87 ± 6.88%
		C4.5	94.62 ± 4.54%	94.62 ± 4.54%
		LMT	100.00 ± 0.00%	100.00 ± 0.00%
	M of N-3-7-10 data set	Naive Bayes SMO	89.33 ± 1.56% 100.00 ± 0.00%	89.33 ± 1.56% 100.00 ± 0.00%
		C4.5	100.00 ± 0.00%	100.00 ± 0.00%
		LMT	86.53 ±3.87%	85.22 ± 4.04%
	Crx data set	Naive Bayes	86.53 ±3.87%	85.51 ± 4.25%
		SMO	86.53 ±3.87%	85.80 ± 3.71%
		C4.5	85.36 ± 4.12%	85.51 ± 4.25%
		LMT	100 ± 0.00%	100.00 ± 0.00%
	Monk-1 data	Naive Bayes	72.22 ± 6.33%	74.97 ± 1.95%
	set	SMO	72.22 ± 6.33%	75.02 ± 5.66%
		C4.5	100.00 ± 0.00%	100.00 ± 0.00%
		LMT	67.14 ± 0.61%	67.36 ± 1.17%
	Monk-2 data	Naive Bayes	67.14 ± 0.61%	67.14 ± 0.61%
	set	SMO	67.14 ± 0.61%	67.14 ± 0.61%


Petia Koprinkova-Hristova

Nikolay Tontchev
email: tontchev@vtu.bg

Silviya Popova
email: popova@icsr.bas.bg

Neural Networks Approach to Optimization of Steel Alloys Composition

Keywords: steel alloys, neural networks, modeling, optimization

The paper presents modeling of steels strength characteristics in dependence from their alloying components quantities using neural networks as nonlinear approximation functions. Further, for optimization purpose the neural network models are used. The gradient descent algorithm based on utility function backpropagation through the models is applied. The approach is aimed at synthesis of steel alloys compositions with improved strength characteristics by solving multi-criteria optimization task. The obtained optimal alloying compositions fall into martenzite region of steels. They will be subject of further experimental testing in order to synthesize new steels with desired characteristics.

Introduction

Production of high strenght steel alloys is of big importance for the modern methalurgy. The main aim is to obtain high quality materials reducing quantitiy of used expensive compounds. Thus the production of steel alloys starts with optimization of number and content of used alloying components that improve their quality at reasonable price. Since the dependences between input and output variables in that case are strongly nonlinear, application of nonlinear modelling and optimization techniques must be applied.

The fact that neural networks are universal approximations of complex nonlinear dependences that apply "black-box" modeling approach [START_REF]Neural Networks for Instrumentation, Measurement and Related Industrial Applications[END_REF][START_REF] Cichoski | Neural Networks for Optimization and Signal Processing[END_REF][START_REF] Rumelhart | Parallel Distributed Processing 1[END_REF] is well known. Therefore, they are proper candidates for modeling structure of such MIMO models. Another useful characteristic of neural networks are their training procedures that are in fact optimization algorithms aimed at minimization of error at neural network output with respect to network connection weights [START_REF] Rumelhart | Parallel Distributed Processing 1[END_REF][START_REF] Werbos | Backpropagation Through Time: What It Does and How to Do It[END_REF]. The well-known backpropagation algorithm is procedure for propagating of derivatives of given function of network output backwards to the input [START_REF] Werbos | Backpropagation Through Time: What It Does and How to Do It[END_REF]. Thus neural network training procedures offer a common approach to optimization tasks in process optimization and control applications [START_REF] Nguyen | Neural Networks for Self-Learning Control Systems[END_REF][START_REF] Werbos | Backpropagation Through Time: What It Does and How to Do It[END_REF][START_REF] Werbos | An Overview of Neural Networks for Control[END_REF].

The application of intelligent modelling approaches such as neural networks is relatively new in this area but there are some examples in the literature concerning titanium alloys [START_REF] Malinov | Software Products for Modeling and Simulation in Material Science[END_REF][START_REF] Malinov | Modelling and Correlation between Processing Parameters and Properties of Titanium Alloys using Artificial Neural Network[END_REF] as wll as several for different steel alloys compositions and different strength characteristics accounting [START_REF] Bhattacharyya | Development and Characterization of C-Mn-Al-Si-Nb TRIP Aided Steel[END_REF][START_REF] Dini | Predicting of Mechanical Pproperties of Fe-Mn-(Al, Si) TRIP/TWIP Steels using Neural Network Modelling[END_REF][START_REF] Esfahani | Modeling the Yield Strength of Hot Strip Low Carbon Steels by Artificial Neural Network[END_REF].

In our previous works the approach described in [START_REF] Werbos | Backpropagation Through Time: What It Does and How to Do It[END_REF] was successfully applied to dynamic optimization [START_REF] Koprinkova | Optimal Control of Batch Biotechnological Processes using Neural Network Model[END_REF] as well as to optimization of cultural media composition and initial conditions [START_REF] Koprinkova-Hristova | Neural Network Optimization of Initial Conditions of Milk Starter Culture Cultivation[END_REF] of two kinds of biotechnological process that are other examples of highly nonlinear systems. In [START_REF] Koprinkova-Hristova | Neural Networks for Mechanical Characteristics Modeling and Compositions Optimization of Steel Alloys[END_REF][START_REF] Tontchev | Comparative Study on Intelligent and Classical Modeling and Composition Optimization of Steel Alloys[END_REF] we applied that approach to optimization of steel composition having 11 alloying elements aimed at simultaneous maximization of 6 steel strength characteristics; we also showed that neural networks outperform non-linear regression models. However, it appeared that having that much input and output variables do not allow obtaining neural network model able to fit with the same accuracy all output variables and solving of six criteria optimization task showed some controversial results.

In the present paper, we continue our work with increased data set and accounting for lower number of variables. Moreover, we have created separate neural network model for each one of the considered strength characteristic. The optimization criteria were simplified and multi-criteria optimization task from [START_REF] Koprinkova-Hristova | Neural Networks for Mechanical Characteristics Modeling and Compositions Optimization of Steel Alloys[END_REF][START_REF] Tontchev | Comparative Study on Intelligent and Classical Modeling and Composition Optimization of Steel Alloys[END_REF] was restricted to two simplified sub-tasks. Comparison between optimized steel compositions with respect to their mechanical characteristics was made. The obtained results will be subject of experimental proof further.

Experimental Data Set

Here we used data base containing 91 steel alloys with their 8 alloying elements and the corresponding strenght characterestics tested after thermal treatment of the steels: Rm -tensile strength; Re -yield strength; A -elongation; Z -reduction of area; HB -Brinell hardness. Tables 1 and2 below summarize the minimal and maxinmal values of all considered variables. For the sick of simplicity, names (X i for alloying elements and Y i for strength characteristics) are given to all considered variables in order to be easier to refer to them further in the text.

Neural Network Models of Steel Characteristics

First task of the present study was to train neural network models that approximate dependence between amounts of alloying elements in the steel for each of the five considered strength characteristics. For that purpose, the data are scaled in proper interval [0, 1]. The used neural network structure (presented on Figure 1) is multilayered without feedback connections since the modeled dependences are static. The neurons transfer function is log sigmoid. Training procedure is resilient backpropagation. Different in number of layers and hidden neurons neural network structures were tested and the better one is chosen. Based on previous investigations [START_REF] Koprinkova-Hristova | Neural Networks for Mechanical Characteristics Modeling and Compositions Optimization of Steel Alloys[END_REF][START_REF] Tontchev | Comparative Study on Intelligent and Classical Modeling and Composition Optimization of Steel Alloys[END_REF] We have relatively small database and we do not have guarantee that it includes all possible combinations between alloying elements. That is why for assessment of the generalization error of neural network models we apply k-fold cross-validation. In our case we split database into k=18 data subsets each containing five data samples. Next, we train 18 neural network models for each strength characteristic using all except one data subset. The left subset is further used for generalization error assessment.

The obtained training and testing mean square errors (MSE) for all five neural network models are given in Table 3 below. Since for some of the strength characteristics there were no data for all alloying element combinations (Y 4 and Y 5 ) smaller number of subsets was created for them (17 and 16).

For further optimization purpose, we choose to use the models with smallest testing errors as follows: models number one for Y 1 and Y 2 , model number 13 for Y 3 , model number 9 for Y 4 , and model number 5 for Y 5 .

Optimization Procedure

Since the explored input/output space is multi-dimensional and the modeling function is highly nonlinear in order to find optimal values of input variables with respect to given quality criteria that comprises output variables, there is need to explore whole region of the input space. However, because of big number of possible combinations the exhaustive search on whole variables space will take too much time. Because of this, we applied gradient optimization technique starting from several different points of input variables surface and compare the obtained results.

Figure 2 below presents the optimization procedure scheme adopted from the socalled "backpropagation of utility" method [START_REF] Werbos | Backpropagation Through Time: What It Does and How to Do It[END_REF]. The optimization task here is defined as follows: find values of input vector X that minimize/maximize the utility function:

( ) Y X J J , = (1)
Here Y is vector of output variables that are related to the input once by a given function (model) F as follows:

( ) p X F Y , = (2)
Here, p is model parameters vector. The optimization procedure needs calculation of utility function gradients with respect to the optimized variables as follows:

X F Y J X J dX dJ dX ∂ ∂ * ∂ ∂ + ∂ ∂ = = (3)
In case when J does not depend explicitly on X, the first term in equation ( 3) is zero and thus gradient depends only on function F.

The layered neural networks structure offers a convenient way for calculating derivatives from equation (3) because the backpropagation method [START_REF] Rumelhart | Parallel Distributed Processing 1[END_REF][START_REF] Werbos | Backpropagation Through Time: What It Does and How to Do It[END_REF] was developed initially as procedure for error derivatives calculation and their "propagation" from the output to the input of the network. From a more common point of view it is method for a given function derivative calculations with respect to variables of an ordered system of equations [START_REF] Werbos | Backpropagation Through Time: What It Does and How to Do It[END_REF]. Hence, it could be applied to any optimization problem that can be described in appropriate way. So application of collection of neural network models trained before as model function F within optimization scheme from Figure 2 allows easy gradients calculation using backpropagation method. Next, any gradient iterative optimization procedure could be applied to find optimal values of the variables X using calculated derivatives as follows:

Model

Y=F(X) X Y J(X, Y) dJ/dX

i X i X i X Δ * ± - = α 1 (5)
Here α is parameter called learning speed and is step-change of X for the i i X Δ th iteration calculated as follows:

( ) i dX g i X = Δ (6)
Here g is some function of optimized variable derivative. Usually it is proportional to the derivative dX i but it can also depend on the old values of .

i X Δ

In present study, we used simple gradient optimization procedure with identity function for g, i.e.

. The learning speed α is set to a relatively small value and stopping criteria is very small change of performance function gradient.

i dX i X = Δ
Two different optimization tasks were solved here. First, one is aimed at maximization of Y 2 that is yield strength Re. The second task includes two criteria maximization: maximal yield strength and maximal ratio between Y 1 (tensile strength Rm) and Y 2 . In both cases, a restriction was added to obtain low carbon steels: X 1 ≤0.3. The two tasks optimization results are presented and compared in the following sections.

First Optimization Task

First, we solved one criteria optimization task:

max 2 → = Y J (7)
We explore 15 different initial steel compositions. The obtained after optimization new 15 alloying elements "recipes" are shown in Table 4 below. The corresponding five strength characteristics of the new steel compositions and the needed number of iterations in each optimization run are presented in Table 5. The maximal value of criterion [START_REF] Koprinkova-Hristova | Neural Network Optimization of Initial Conditions of Milk Starter Culture Cultivation[END_REF] was obtained in the case 4 (emphasized in bold in the table). The number of the iterations obviously depends on the starting steel composition position on the search surface; it varies between 160 and 530 iterations. In all cases, the obtained new compositions belong to the martenzite steels subset (see Figure 3 below). All the obtained new steels have relatively close values of the optimization criterion. The final decision, which composition suits better for a given purpose, should be taken considering all strength characteristics from Table 5. For some of these characteristics (e.g. 4 and 5) differences are not so big but for characteristics 1 and 3 there are significant differences observed.

Second Optimization Task

Here in addition to the first criterion [START_REF] Koprinkova-Hristova | Neural Network Optimization of Initial Conditions of Milk Starter Culture Cultivation[END_REF] the second one was added as follows: Here we start again with the same 15 initial steel compositions. The obtained after optimization new 15 "recipes" are shown in Table 6. The corresponding to them strength characteristics iterations needed and two optimization criteria values are presented in Table 7. Again, obtained new steels are in the martenzite region on Figure 3 but there are bigger differences among them.

1 , max 2 2 1 1 , 2 1 2 1 2 2 1 = = → + = = = a a J a J a J Y Y J Y J ( 8
)
Because in this task we have two controversial criteria -one targeted towards increasing of yield strength and the other having Y 2 in denominator, i.e. targeted towards decreasing it -the final decision about the best obtained "recipie" is not definitive. The best steel composition with respect to the first citerion is number 8 while with respect to the second -number 14 (all in bold in the table). This can explain the big variety of iterations number needed for individual optimization runsfrom 3 to 2516. The variety of combinations between all stregth characteristics is also bigger making final decision more complicated too.

Results Analysis

In order to compare two optimization tasks results we present them in the form of steel types diagram accepted in metallurgy -Figure 3. It consists of Cr and Ni equivalent values calculated according to the equation on diagram axis. Here Nb denotes element niobium that is not present in our recipes, i.e. Nb=0. According to the equivalent Cr and Ni values, steels are derived in several classes shown on the diagram. The bigger is the first of the optimized characteristics -yield strength -the bigger is the steel strengh as a holle. From the other hand the ratio tensile strength/ yield strength (the second criterion) is related to local overloading steel resistance. Obtaining the good balance between these two major steel characteristics is difficult in practice. The imposed restiriction in both optimization task (low carbon concentration) was also another characteristic that guarantees obtaining of high strength steels.

From the Figure 3 we conclude that all the obtained new steel compositions belong to the targeted martenzite area. All the results from the first task are concentrated in the middle of that area while those from the second one diverge towards area borders. This result is expected having in mind controversal criteria of the second optimization task.

Conclusions

The application of neural networks for modeling of the complex dependences of steel compositions characteristics from the relatively big number of alloying elements considered demonstrated their power as non-linear MIMO functions approximations. They allow us to carry out further optimization experiments by simulations in order to spent expensive laboratory synthesis and analysis of different steel compositions.

The results of the solved above optimization tasks are meaningful and will help further creation of high strength steel compositions for machine-building industry.

Our intensions for our future work are to include some other optimization criteria as well as to test the obtained steel "recipes" in practical experimental investigations in order to prove and refine their quality.

Fig. 1 .

 1 Fig. 1. Neural network models structure

Fig. 2 .

 2 Fig. 2. Optimization procedure; dashed lines represent gradients calculation

Fig. 3 .

 3 Fig. 3. Optimization results: squares represent the first optimization task solutions; circles represent the second optimization task solutions

Table 1 .

 1 Alloying elements minimal and maximal values.

	var. name X 1	X 2	X 3	X 4	X 5	X 6	X 7	X 8
	element	C	Si	Mn	Ni	S (P)	Cr	Mo	V
	min	0.120 0.100 0.300	0	0.015 0.150	0	0
	max	0.520 1.400 1.750 4.220 0.035 3.250 1.500 0.150

Table 2 .

 2 Steel strength characteristics minimal and maximal values.

	var. name	Y 1	Y 2 Y 3 Y 4 Y 5
	characteristic Rm	Re	A Z HB
	min	500 300 7 30 179
	max

Table 3 .

 3 Training and testing errors for all NN models.

	No	Y 1	Y 1	Y 2	Y 2	Y 3	Y 3	Y 4	Y 4	Y 5	Y 5
		train	test	train	test	train	test	train	test	train	test
		MSE	MSE	MSE	MSE	MSE	MSE	MSE	MSE	MSE	MSE
	1	0.008 0.055 0.033 0.060 0.004 0.083 0.022 0.115 0.003 0.126
	2	0.007 0.445 0.015 0.224 0.006 0.143 0.007 0.200 0.006 0.197
	3	0.013 0.099 0.012 0.361 0.004 0.108 0.016 0.299 0.005 0.063
	4	0.013 0.064 0.010 0.112 0.009 0.079 0.012 0.137 0.006 0.099
	5	0.019 0.171 0.019 0.125 0.005 0.065 0.008 0.166 0.003 0.050
	6	0.010 0.170 0.015 0.217 0.006 0.088 0.009 0.115 0.005 0.135
	7	0.012 0.130 0.010 0.086 0.006 0.112 0.006 0.160 0.004 0.063
	8	0.012 0.104 0.017 0.083 0.005 0.162 0.005 0.382 0.004 0.075
	9	0.014 0.131 0.010 0.114 0.007 0.058 0.008 0.091 0.005 0.101
	10 0.010 0.087 0.012 0.222 0.007 0.135 0.025 0.213 0.005 0.085
	11 0.023 0.127 0.018 0.156 0.009 0.112 0.008 0.168 0.004 0.065
	12 0.013 0.136 0.006 0.085 0.010 0.067 0.013 0.157 0.006 0.140
	13 0.013 0.341 0.010 0.397 0.007 0.045 0.010 0.316 0.006 0.124
	14 0.009 0.195 0.010 0.159 0.006 0.108 0.019 0.289 0.005 0.089
	15 0.008 0.100 0.009 0.199 0.007 0.114 0.011 0.243 0.006 0.107
	16 0.010 0.080 0.008 0.124 0.004 0.157 0.009 0.123 0.004 0.118
	17 0.007 0.354 0.009 0.342 0.005 0.153 0.008 0.484
	18 0.009 0.091 0.019 0.089 0.007 0.104

Table 4 .

 4 Optimized steel compositions from the first task.

	No	X 1	X 2	X 3	X 4	X 5	X 6	X 7	X 8
	1 0.1562 0.9180 0.4517 0.9215 0.0202 0.8730 0.2450 0.0410
	2 0.1562 0.9137 0.4517 0.8663 0.0208 0.8740 0.2757 0.0423
	3 0.1562 0.7997 0.7740 1.4991 0.0184 1.1572 0.1943 0.0087
	4 0.1562 1.0861 0.8087 1.9735 0.0243 1.7111 0.1733 0.0087
	5 0.1584 1.2125 0.8475 2.5041 0.0274 2.0026 0.1815 0.0465
	6 0.1603 0.8917 1.3092 4.1876 0.0381 1.8100 0.6971 0.1025
	7 0.1562 1.0171 1.4663 4.0238 0.0446 2.1378 1.0380 0.1023
	8 0.1562 1.0466 1.4833 4.1469 0.0446 2.2158 1.0497 0.1146
	9 0.1562 1.0566 1.4695 4.1876 0.0446 2.2936 1.0486 0.1263
	10 0.1562 1.0589 1.4319 4.1876 0.0446 2.2193 0.9355 0.1440
	11 0.1562 0.9196 0.4517 0.9113 0.0203 0.8721 0.2501 0.0416
	12 0.1870 0.2352 0.6480 2.5490 0.0446 1.4195 0.9603 0.1583
	13 0.1566 0.2352 0.7037 3.1023 0.0446 1.7221 0.9820 0.1583
	14 0.2714 1.1062 1.0644 2.3624 0.0219 1.0496 0.1572 0.0087
	15 0.2673 1.2441 1.1861 2.7854 0.0269 1.2693 0.2301 0.0218

Table 5 .

 5 Optimized steel characteristics from the first task, optimization criteria and iterations needed

	No	Y 1	J 1 =Y 2 Y 3 Y 4	Y 5 iterations
	1 1643.6 1350.4 8.5 54.7 263.8	432
	2 1640.2 1349.9 8.8 54.8 260.6	445
	3 1407.0 1369.5 7.8 49.5 272.6	261
	4 1638.0 1372.6 8.1 49.8 282.7	437
	5 1621.3 1359.9 9.7 54.5 290.2	530
	6 1122.8 1362.1 13.7 49.9 288.2	276
	7 1394.7 1360.2 16.0 42.5 282.3	378
	8 1294.9 1360.7 15.0 47.2 280.8	427
	9 1227.7 1360.9 14.6 51.2 279.7	362
	10 1246.9 1365.3 11.9 53.1 279.8	523
	11 1642.9 1350.4 8.6 54.7 263.5	442
	12 662.7 1341.0 24.5 55.0 208.9	312
	13 841.0 1336.3 25.0 55.0 222.5	200
	14 1666.9 1370.3 11.4 51.1 281.6	160
	15 1647.5 1370.9 12.2 54.1 287.4	214

Table 6 .

 6 Optimized steel compositions from the second task.

	No	X 1	X 2	X 3	X 4	X 5	X 6	X 7	X 8
	1 0.1562 0.6937 0.5036 0.5673 0.0094 0.4832 0.3116 0.0301
	2 0.1577 0.6915 0.5092 0.5807 0.0135 0.4832 0.2698 0.0225
	3 0.1562 1.2231 0.7684 1.6164 0.0196 0.4832 0.1572 0.1255
	4 0.1562 0.9408 0.7296 1.1316 0.0249 0.4832 0.1602 0.1459
	5 0.1562 0.8331 1.7453 0.4564 0.0228 1.9982 0.5520 0.0632
	6 0.1562 1.2004 1.7453 0.4564 0.0446 0.4832 1.4948 0.0087
	7 0.1562 1.0013 1.6050 2.4654 0.0416 2.0609 1.4905 0.0570
	8 0.2964 1.1273 1.4969 3.2602 0.0375 2.6033 1.2407 0.1245
	9 0.2980 1.2624 1.6303 3.7373 0.0413 2.9228 1.3716 0.1427
	10 0.3000 1.3968 1.7453 4.1874 0.0446 3.2287 1.4948 0.1583
	11 0.1562 0.2352 0.4517 0.4569 0.0094 0.4832 0.1650 0.0093
	12 0.1562 0.7446 1.4706 1.9032 0.0442 0.4832 0.9263 0.1096
	13 0.1562 0.5792 1.1822 1.6692 0.0387 0.4832 1.1307 0.1441
	14 0.3000 1.0140 1.1720 2.1100 0.0210 1.0760 0.2940 0.0087
	15 0.2997 1.1352 1.3479 2.5514 0.0247 1.3766 0.4125 0.0246

Table 7 .

 7 Optimized steel characteristics from the second task, optimization criteria and iterations needed

	No	Y 1	J 1 =Y 2 Y 3 Y 4	Y 5 J 2 =Y 1 /Y 2 iterations
	1 1664.2 1233.3 7.6 52.5 261.0 1.3494	865
	2 1665.5 1285.3 8.1 52.8 239.3 1.2958	157
	3 1659.4 957.9 9.3 54.6 271.2 1.7323	266
	4 1456.2 797.5 7.9 54.3 250.9 1.8260	605
	5 499.1 300.1 11.4 40.7 185.3 1.6631	3
	6 499.1 325.4 22.6 41.5 179.1 1.5338	6
	7 735.5 609.3 23.3 52.5 289.3 1.2071	1353
	8 1463.1 681.1 11.3 42.8 274.6 2.1481	462
	9 1488.5 853.6 14.6 41.0 273.8 1.7438	839
	10 1499.8 829.8 20.3 40.8 271.4 1.8074	2
	11 1571.6 916.0 8.2 54.1 227.0 1.7157	289
	12 1199.5 577.6 13.4 53.2 207.8 2.0767	607
	13 760.0 615.1 12.2 54.7 282.1 1.2356	781
	14 1635.4 1326.0 15.5 47.7 258.9 1.2333	3
	15 1613.4 1308.0 15.6 53.4 267.2 1.2335	2516

Acknowledgments. This work was partially supported by the Bulgarian National Science Fund under the Project No DDVU 02/11 "Characteristic modeling and composition optimization of iron-based alloys used in machine-building industry".


Milan Cisty
email: milan.cisty@stuba.sk

Determining Soil -Water Content by Data Driven Modeling when Relatively Small Data Sets are Available

Keywords: soil water regime, pedotransfer function, neural networks, support vector machines, harmony search

A key physical property used in the description of a soil-water regime is a soil water retention curve, which shows the relationship between the water content and the water potential of the soil. Pedotransfer functions are based on the supposed dependence of the soil water content on the available soil characteristics. In this paper, artificial neural networks (ANNs) and support vector machines (SVMs) were used to estimate a drying branch of a water retention curve. The performance of the models are evaluated and compared in case study for the Zahorska Lowland in the Slovak Republic. The results obtained show that in this study the ANN model performs somewhat better and is easier to handle in determining pedotransfer functions than the SVM models.

Introduction

Modeling water content and transport in soil has become an important tool in simulating agricultural productivity as well as in solving various hydrological tasks. For instance, optimum irrigation management requires a systematic estimation of the soil-water status to determine both the appropriate amounts and timing of irrigation. That is why soil characteristics appear as a important input in the numerical simulation of a soil-water regime. A relatively large number of works have appeared which were devoted to determining the water retention curve which is needed for this purpose from more easily available soil properties such as particle size distribution, dry bulk density, organic C content, etc., e.g. [START_REF] Baker | Optimisation of pedotransfer functions using an artificial neural network ensemble method[END_REF], [START_REF] Bouma | Using Soil Survey Data for Quantitative Land Evaluation[END_REF], [START_REF] Minasny | The Neuro-M Methods for Fitting Neural Network Parametric Pedotransfer Functions[END_REF], etc. Pedotransfer functions (PTF) have become the term for such relationships between soil hydraulic parameters and the more easily measurable properties usually available from a soil survey [START_REF] Baker | Optimisation of pedotransfer functions using an artificial neural network ensemble method[END_REF]. Consequently, the method for the quantification of these relationships uses various types of regression analyses. The aim of this paper is a comparison of three regression models for determining pedotransfer functions.

Besides the standard regression methods, artificial neural networks (ANNs) have become the tool of choice in developing PTFs, e.g., [START_REF] Baker | Optimisation of pedotransfer functions using an artificial neural network ensemble method[END_REF], [START_REF] Mohammadi | Testing an artificial neural network for predicting soil water retention characteristics from soil physical and chemical properties[END_REF], [START_REF] Schaap | Neural network analysis for hierarchical prediction of soil hydraulic properties[END_REF], [START_REF] Tamari | Testing an Artificial Neural Network for Predicting Soil Hydraulic Conductivity[END_REF], etc.). Authors of above works confirm that they received better results from ANN-based pedotransfer functions than from standard linear regression-based PTFs.

Artificial neural networks include the ability to learn and generalize from examples with the aim of providing meaningful solutions to the problems to be solved. This process is called "training". When the training of an ANN is accomplished with a set of input and output data, the aim is to adjust the parameters of the ANN and make the ANN also provide corresponding outputs for other sets of input data (for which the outputs are not known).

Also second data driven method was used in this study-support vector machines (SVMs), which were developed by Vapnik [START_REF] Vapnik | The Nature of Statistical Learning Theory[END_REF] and are gaining in popularity due to their attractive features and promising empirical performance. The formulation embodies the structural risk minimization principle in addition to the traditional empirical risk minimization principle employed by conventional neural networks. It is this difference which gives SVMs a greater ability to generalize, which is the goal of statistical learning.

The objective of this work is to compare abovementioned methods while developing PTFs for the Zahorska Lowland in Slovakia, which was selected as a representative region for the investigation (e.g., while solving the regression task of determining the water retention curve from easily available soil properties).

In the following part of the paper ("Methodology") the three methods used in this study -ANN, SVM and multiple linear regression are briefly explained. Then the data acquisition and preparation is presented. In the "Results" part, the settings of the experimental computations are described in detail, and the "Conclusions" of the paper evaluates these experiments on the basis of the statistical indicators.

Methodology

The first approach for modeling the PTFs used in this paper is the application of artificial neural networks (ANNs). This approach has been described e.g. in [START_REF] Minasny | The Neuro-M Methods for Fitting Neural Network Parametric Pedotransfer Functions[END_REF] or [START_REF] Schaap | Neural network analysis for hierarchical prediction of soil hydraulic properties[END_REF]. Briefly summarized, a neural network consists of input, hidden and output layers, which contains processing elements. The number of processing elements in the input layer and output layer correspond to the number of input (e.g., the soil's bulk density, the soil's particle size data, etc.) and output variables of the model. So-called "learning" involves adjustment of the synaptic connections that exist between the neurons or weights in hidden layer, which are used for the transformation of the inputs to the outputs. A type of ANN known as a multi-layer perceptron (MLP), which uses a back-propagation training algorithm, was used for generating the PTFs in this study. The training process was performed by back propagation training algorithm of an MLP. The basic information about the application of an ANN to regression problems is available in the literature and is well known, so we will not provide a more detailed explanation here.

A second approach called support vector machines (SVM) for estimating the pedotransfer functions used in this study is explained hereinafter, with brief explanations of its principles. A more detailed description of the methodology can also be found in available sources, e.g., in [START_REF] Smola | A Tutorial on Support Vector Regression[END_REF].

The architecture of a SVM is similar to that of an ANN, but the training algorithm is significantly different. The basic idea is to project the input data by means of kernel functions into a higher dimensional space called the feature space, where a linear regression can be performed for an originally nonlinear problem, the results of which are then mapped back to the original input-output space. The linear regression is maintained by quadratic programming, which ensures a global optimum and an optimal generalization. The important idea is to fully ignore small errors (by introducing the "tube" variable ε, which defines what the "small" error is) to make the regression sparse, that is, dependent on a smaller number of inputs (called the support vectors), which makes the methodology much more computationally treatable. The uniqueness of solution produced by SVMs is often emphasized, but the actual truth is that this solution is only unique for a given set of performance parameters, which should be chosen and process of selection of them will be described later.

Objective function of mentioned linear regression in feature space SVM simultaneously minimizes both the empirical risk and the model's complexity; the tradeoff between these two goals is controlled by parameter C. An important characteristic of SVMs as a consequence of this form of objective function is that a better ability to generalize could be expected (by choosing appropriate parameter C), compared, e.g., with ANNs, because unnecessarily complex models usually suffer from over fitting.

The radial basis function was chosen on a trial and error basis as the kernel function for this work. This function has the following form:

K(x i , x j ) = exp(-g|| x i -x j || 2 ), g > 0 . (1)
The parameter γ of this kernel function, the tube size ε for the ε-insensitive loss function, and parameter C should be found, which the basic task is when SVMs are applied to particular task. The harmony search methodology (HS) was used for this purpose instead of the usual trial-and-error principle. A harmony search is a metaheuristic search algorithm introduced by Geem [START_REF] Geem | Various Continuous Harmony Search Algorithms for Web-Based Hydrologic Parameter Optimisation[END_REF] and is inspired by the improvisational process of musicians. In an HS algorithm, each musician corresponds to one decision variable. A musical instrument's pitch range corresponds to a decision variable's value range; the musical harmony at a certain time corresponds to a solution vector at certain iteration; and the audience's aesthetics corresponds to an objective function. Just as a musical harmony is improved time after time, a solution vector is improved iteration by iteration by the application of the improvisation's operators (the random selection of a tone, a musician's memory considerations or a pitch adjustment).

As a criterion for selecting the appropriate combinations of the parameters, the correlation coefficient could be used in regression task of determining PTF as the value of the objective function of the harmony search methodology.

Study Area and Data Collection

The data used in this study were obtained from a previous work [START_REF] Skalova | Pedotransfer Functions of the Zahorska Lowland Soils and Their Application to Soil-Water Regime Modeling[END_REF]. An area of the Zahorska Lowland in Slovakia was selected for testing the methods described. A total of 140 soil samples were taken from various localities in this area.

The soil samples were air-dried and sieved for a physical analysis. A particle size analysis was performed utilizing Cassagrande's methods. The dry bulk density, particle density, porosity and saturated hydraulic conductivity were also measured on the soil samples. The points of the drying branches of the WRCs for the pressure head values of -2.5, -56, -209, -558, -976, -3060 and -15300 cm were estimated using overpressure equipment.

A full database of the 140 samples and their properties was used for creating the input data for the modeling from which the three subsets of the data were produced: training data (88 data samples), validation data (22 data samples), testing data (30 data samples).

A practical way to find a better generalization data driven model is to set aside a small percentage of the training set and use it for the cross validation. When the error in the validation set increases, the training should be stopped. The five divisions of the data used to the training and validation data set were used.

An ensemble of data-driven modeling was used in the present work, which means a collection of a finite number of data-driven models that are trained for the same task. This is meant as a simple variant of bootstrapping (the bootstrap scheme involves generating subsets of the data on the basis of random sampling with replacements as the data are sampled). Five data-driven models are trained independently, and their predictions are combined for the sake of obtaining a better generalization (average value was taken as result). For this reason the mentioned training fraction of the data (110 samples) was divided into the training and validation data sets alternatively in five different versions.

Results

The first approach used in determining the water retention curves in the presented work was an application of ensemble neural networks. The same network architecture for every ANN in the network was determined. In this work the multilayer perceptron (MLP) with 2, 3, and 4 neurons in the hidden layer was tested; an MLP with 3 neurons in the hidden layer was finally chosen for the ensemble neural network model. A neuron with a bias and tanh activation function was used. The Levenberg-Maquardt method was used in the context of the back propagation method. The networks were trained for computing the water content at the pressure head value h w = -2.5, -56, -209, -558, -976, -3060, -15300 cm. Then the testing dataset was computed with the ensemble ANN. The results with five regression coefficients are summarized in Table 1 also with final (average) regression coefficient R.

Given regression problem was also solved by using ensemble of support vector machines. The estimation of the steps of the SVM regression (described in the methodology part of this paper) are the following: 1) the selection of a suitable kernel and the appropriate kernel's parameter (g in eq.1); 2) specifying the ε parameter and specifying the capacity C.

As a criterion for selecting the appropriate combinations of the parameters, the correlation coefficient for the training and cross-validation data is calculated within the objective function of the harmony search, where the correlation coefficient of the cross-validation data was weighted by the coefficient 1.2 for the sake of a better generalization.

In the training phase, SVM models for a pressure head value of h w = -2.5, -56, -209, -558, -976, -3060, -15300 cm were created. This was repeated five times because of the five divisions of the data used to train the model on the training and validation data set. A total of 35 computations were run because there is seven variables computed. Then the testing dataset was computed five times with the models obtained, and the final result is the average of the outputs from these five models; the results are summarized with the regression coefficients in Table 2. (

) 2
where θ hw is the water content [cm 3 .cm -3 ] for the particular pressure head value h w [cm]; 1 st cat., 2 nd cat., 3 th cat. and 4 th cat. are the percentages of the clay (d< 0.01 mm), silt (0.01-0.05 mm) and sand (0.1-2.0 mm); r d is the dry bulk density [g.cm -3 ]; and a, b, c, d, e, f are the parameters determined by the regression analysis. In the case of the multi-linear regression it was not possible to use ensemble models, because in this case no iterative process is applied which involves cross validation, so all 110 training samples were used as a whole for the development of the model and 30 samples for testing.

The PTFs designed were evaluated on a testing dataset. The results of the multilinear regression are listed in Table 3. As can be seen, the results from both data driven techniques are clearly better compared with the multi-linear regression and from SVM are somewhat worse compared with the ANN. From these results, it seems that ANNs are more resistant to an insufficient amount of data (which is the case in this work), because, on the other hand, better results with the application of the SVM than with the ANN for the PTF evaluation were reported in the literature [START_REF] Twarakavi | Development of Pedotransfer Functions for Estimation of Soil Hydraulic Parameters Using Support Vector Machines[END_REF]. It should be mentioned that the authors of mentioned paper worked with larger data sets (2134 soil samples). For this reason the authors of the present paper hypothesize that it is advisable to use combined SVM/MLP models, because of the variability of an adequate methodology, but this should be verified in future work.

Conclusions

The results of this paper contain a description and evaluation of the models of an ensemble of multi-layer perceptrons and an ensemble of support vector machines for the development of pedotransfer functions for the point estimation of the soil-water content for the seven pressure head values h w from the basic soil properties (particlesize distribution, bulk density). Both ensemble data-driven models were compared to a multiple linear regression methodology.

• The accuracy of the predictions was evaluated by the correlation coefficient (R) between the measured and predicted parameter values. The R varied from 0.631 to 0.888 for the multi-linear regression, from 0.864 to 0.930 for the MLP, and from 0.816 to 0.907 for the SVM. The MLP models perform somewhat better than the SVM models. Nevertheless, the results from both data-driven models are quite close, and the results show that they provide a significantly more precise outcome than traditional multi-linear regression. • Although SVM training is faster, the whole process of ANN training for evaluating PTFs is accomplished in less time, because of the ability of ANNs to produce more outputs in one run, which is the advantage versus SVMs.

Because other authors have reported the better regression ability of SVMs compared with ANNs [START_REF] Twarakavi | Development of Pedotransfer Functions for Estimation of Soil Hydraulic Parameters Using Support Vector Machines[END_REF], the authors of the present paper hypothesize that it is advisable to use combined SVM/MLP models, because of this variability in suitable methodology. This should be verified in future work. The authors of the mentioned paper worked with larger data sets (they used 2134 soil samples; 140 samples were used in our work), and the influence of the amount of data or other statistical data set properties on the choice of the methodology suitable to use should be evaluated. Also other types of data driven models should be tested, e.g., generalized regression neural network or radial basis network.

θ

  hw = a*1 st cat. + b*2 nd cat. + c*3 th cat. + d*4 th cat. + e*r d + f .

Table 1 .

 1 Regression coefficients (R1-5) of five ANN ensemble models and resulted R

	h w [cm]	R1	R2	R3	R4	R5	R
	-2.5	0.914	0.921	0.937	0.934	0.888	0.930
	-56	0.905	0.922	0.934	0.897	0.872	0.916
	-209	0.886	0.92	0.934	0.897	0.871	0.912
	-558	0.883	0.927	0.943	0.879	0.877	0.912
	-976	0.872	0.923	0.937	0.871	0.865	0.905
	-3060	0.858	0.92	0.93	0.846	0.848	0.892
	-15300	0.192	-0.270	-0.281	-14.718	54.796	0.864

Table 2 .

 2 Regression coefficients (R1-5) of five SVM ensemble models and resulted R

	h w [cm]	R1	R2	R3	R4	R5	R
	-2.5	0.910	0.908	0.902	0.910	0.891	0.907
	-56	0.865	0.859	0.863	0.848	0.852	0.861
	-209	0.860	0.857	0.867	0.862	0.855	0.863
	-558	0.856	0.854	0.857	0.857	0.848	0.857
	-976	0.833	0.842	0.846	0.840	0.834	0.846
	-3060	0.845	0.847	0.852	0.846	0.837	0.851
	-15300	0.799	0.822	0.819	0.815	0.817	0.816
	A multi-linear regression for assessing the PTFs was accomplished for comparison
	and it was used in the form:

Table 3 .

 3 Results of the multi-linear regression (coefficients of Eq. 2 and regression coeficient)

	hw [cm]	a	b	c	d	e	f	R
	-2.5	-0.224	-0.427	-0.2728	-0.403	-37.251	133.466	0.888
	-56	-0.803	-1.157	-0.905	-1.223	-25.193	180.218	0.798
	-209	-0.794	-1.163	-1.060	-1.2781	-17.761	167.584	0.819
	-558	-0.531	-0.932	-0.833	-1.024	-19.859	143.738	0.856
	-976	-1.7	-2.091	-1.984	-2.174	-18.405	255.177	0.685
	-3060	-1.166	-1.624	-1.499	-1.655	-1.166	199.374	0.770
	-15300	-2.007	-2.459	-2.288	-2.487	-14.882	275.784	0.631

Acknowledgement

This work was supported by the Slovak Research and Development Agency under Contract No. LPP-0319-09 and APVV-0496-10 and by the Scientific Grant Agency of the Ministry of Education of the Slovak Republic and the Slovak Academy of Sciences, Grants No. 1/1044/11.