Law Discovery using Neural Networks

1 Introduction

The discovery of a numeric law from a set of data is the central part of scientific discovery systems. Such systems, for example, can detect a relationship between the distance r to the sun and the revolution period T of five planets known to Kepler's third law T = kr 3 1 2 (k is a constant).

After the pioneering work of the BACON systems [Langley, 1978; Langley et al. , 1987], several methods [Langley and Zytkow, 1989; Falkenhainer and Michal­ ski, 1990; Langley and Zytkow, 1990; Schaffer, 1993; Sutton and Matheus, 1991] have been proposed. The ba­ sic search strategy employed by these methods is much the same: two variables are recursively combined into a new variable by using multiplication, division, or some predefined prototype functions . BACON and FAHREN­ HEIT [Langley and Zytkow, 1989] use trend detectors to combine variables and employ a heuristic form of depth-first search. ABACUS [Falkenhainer and Michal­ ski, 1990] creates a proportional graph and performs a modified beam search. In IDS [Langley and Zytkow, 1990], correlation analysis is applied, and a beam search is performed. The E* algorithm [Schaffer, 1993] con­ siders only bivariate functions . Also, Sutton-Matheus' algorithm [Sutton and Matheus, 1991] performs a re­ gression, and the correlation between the squared error and the square of the variable's value [Sanger, 1991] is used to combine variables.

These existing methods suffer from the following prob­ lems: first , because combining two variables into a new one must be done in order, a combinatorial explosion may occur when complex laws are sought for data con­ sisting of a large number of variables, or the desired laws will be missed when some heuristic search param­ eters are inappropriate. Second, when some powers ap­ pearing in a law are not restricted to integers, the law may remain unknown unless some appropriate prototype functions such as r 3 1 2 are prepared in advance. How­ ever, a priori information is rarely available. Third, these methods are often criticized for their lack of ro­ bustness; noise tolerance is definitely required since real observed data contain noise [Langley and Zytkow, 1990; Schaffer, 1993].

We believe a connectionist approach has great poten­ tial to solve the above problems. In order to directly learn a generalized polynomial term whose power val­ ues are not restricted to integers, a computational unit called a product unit has been proposed [Durbin and Rumelhart , 1989]; instead of calculating a weighted sum of input values, this unit calculates a weighted prod­ uct, where each input value is raised to a power deter­ mined by a variable weight. However, serious difficulties have been reported when using standard BP [Rumel­ hart et al. , 1986] to train networks containing these units [Leerink et al., 1995]. Although some heuristic strategies such as multiple learning algorithms have been proposed [Leerink et al., 1995], their improvements over BP have been less than remarkable. Moreover, these earlier stud­ ies dealt only with binary data and did not specifically address the problem of numeric law discovery.

In this paper, we propose a connectionist approach called RF5 for discovering numeric laws. Section 2 explains how neural networks are used to discover a class of numeric laws. Section 3 describes a second-order learn­ ing algorithm called BPQ which trains the neural net­ works described in Section 2. Section 4 explains a crite­ rion for selecting the most suitable candidate out of the trained ones. Section 5 evaluates the proposed method RF5 by doing experiments using artificial data, real data, and time series data.

2 Law Discovery using Neural Nets

This section explains a connectionist problem formaliza­ tion for numeric law discovery, first proposed in [Durbin and Rumelhart , 1989] . Let {(x1, Yl), · · ·, (xm, Ym)} be a set of examples, where Xt is an n-dimensional input vector and Yt is a target value corresponding to Xt · In this paper, a class of numeric laws expressed as

is considered, where each parameter c; or Wij is an un­ known real number and h is an unknown integer. If a target law consists of periodic function or discontinuous function, we cannot exactly discover it when Eq. (1) is assumed. However, if the range of each input variable is bounded, such a law can be closely approximated to a multivariate polynomial function with a finite number of terms. Moreover, since the power values are not re­ stricted to integer, we can expect that the approximated polynomial function has a fewer number of terms. Here­ after, (eo,···, ch)T and (wil, · · ·, Win)T are expressed as c and wi , respectively, where aT means the transposed vector of a. In addition, a vector consisting of all param­ eters, (cT, wf, · · ·, wnT is simply expressed as~' and N(= nh + h + 1) denotes the dimension of~.

By adding an adequate value· to each element of the input vectors, if necessary, without losing generality we can assume Xti > 0; then, Eq. (1) is equivalent to

Equation (2) can be regarded as the feedforward com­ putation of a three-layer neural network where the acti­ vation function· of each hidden unit is exp( s) = e 8 · Here h, wi, and c denote the number of hidden units, the weights between the input units and hidden unit i, and the weights between the hidden units and the output unit, respectively. Hereafter, the output value of hidden unit i is described as

and then the output value of the output unit is described as

The hidden units defined by v(x; w) are called product units [Durbin and Rumelhart, 1989]. The discovery of numeric laws subject to Eq. (1) can thus be defined as the following learning problem in neural networks. That is, the problem is to find the ~ that minimizes the following objective function:

3 BPQ Algorithm

In our early experiments and as reported in earlier stud­ ies [Leerink et al., 1995], the problem of minimizing Eq. (3) turned out to be quite tough. Thus, in order to efficiently and constantly obtain good results, this paper employs a new second-order learning algorithm called BPQ [Saito and Nakano, 1997]; by adopting a quasi­ Newton method [Gillet al., 1981; Luenberger, 1984] as a basic framework, the descent direction, D.~, is calculated on the basis of a partial BFGS update and a reasonably accurate step-length, A, is efficiently calculated as the minimal point of a second-order approximation. In first­ order learning algorithms which calculate the search di­ rection as the gradient direction, a large number of itera­ tions are often required until convergence. On the other hand, in existing second-order methods [Gillet al., 1981; Luenberger, 1984] which converge more quickly by using both gradient and curvature information, it is difficult to suitably scale up for large problems, and a large amount of computation is required for calculating the optimal step-length. BPQ can be reasonably scaled up by in­ troducing a storage space parameter, and the computa­ tional complexity for calculating the optimal step-length is reasonably small, almost equivalent to that of gradient vector evaluation.

For the problem of minimizing Eq. (3), the partial BFGS update can be directly applied, while the ba­ sic procedure for calculating the step-length A must be slightly modified. In the step-length calculation, since A is the only variable, we can express f(~ +All~) sim­ ply as ((A). Its second-order Taylor approximation is given as ((0) + ('(O)A + ~("(O)A 2 · When ('(0) < 0 and ("(0) > 0, the minimal point of this approximation is given by

Here, the method for coping with the other cases is ex­ actly the same as described in [Saito and Nakano, 1997].

For three-layer neural networks defined by Eq. (3), we can efficiently calculate both (' ( 0) and (" ( 0) as follows. By differentiating ((>.) and substituting 0 for >., we ob­ tain

Now that the derivative of z(xt; 4>) is defined as

we obtain

where

and ~Ci and ~Wij are modification values corresponding to Ci and Wij, respectively.

Incidentally, we can employ other second-order learn­ ing algorithms such as SCG [M0ller, 1993] or OSS [Bat­ titi, 1992], but BPQ worked the most efficiently among them in our own experience.

4 Criterion for Selection

In general, for a given set of data, we cannot know the optimal number of hidden units in advance. Moreover, since the data is usually corrupted by noise, the law­ candidate which minimizes Eq. (3) is not always the best one. We must thus consider a criterion to adequately evaluate the law-candidates discovered by changing the number of hidden units. In this paper, by assuming that the target output values are corrupted by Gaussian noise with a mean of 0 and an unknown standard devi­ ation of a, finding an adequate number of hidden units is formalized as a model selection problem of the max­ imum likelihood estimation problem. Thus, we adopt the MDL (Minimum Description Length) criterion [Ris­ sanen, 1989] for this purpose. The MDL fitness value is defined by

where MSE represents the value of the mean squared error defined by

Here, ~ is a set of weights which minimizes Eq. (3), N is the number of parameters in 4>, and m is the number of examples. Hereafter, our discovery method employ­ ing the connectionist problem formalization, the BPQ algorithm and the MDL criterion, is called RF5 (Rule extraction from Facts version 5).

5 Evaluation by Experiments

5.1 Artificial data

The rule discovery method, RF5, was evaluated by using an artificial problem proposed by Sutton and Matheus [Sutton and Matheus, 1991] and our modified version. The original problem is to restore a law described as

while a law for the modified problem is described as

Each example is generated as follows: each value of vari­ ables x1, · · · , X5 is randomly generated in the range of [0, 1], and the corresponding value of y is calculated us­ ing Eq. (6) or (7). In these problem, the total num­ ber of variables is 9 (n = 9); each value of irrelevant variables X6, · · · , Xg is also randomly generated in the range of [0, 1], and the number of examples is set to 200 (m = 200).

In the experiments, the initial values for the weights between the input and hidden units were independently generated according to a normal distribution with a mean of 0 and a standard deviation of 1; the initial val­ ues for the weights between the hidden and output units were set to 0, but the bias value at the output unit was initially set to the average output value of all training examples. The iteration was terminated when any of the following three conditions was met: the MSE value was sufficiently small, i.e.,

the gradient vector was sufficiently small, i.e.,

or the total processing time exceeded 100 seconds.

In the experiments, we changed the number of hidden units from 1 to 3 (h = 1, 2, 3) and performed 100 tri­ als for each of them. Table 1 shows the basic statistics of MSE values, MDL values, iterations, and processing times (sec.). 1 This table shows that the correct number of hidden units, h = 2, was successfully found because the best MDL value was minimized. When h = 2, all100 trials converged to the global minimum for the original problem; however, a few of the 100 trials converged to undesirable local minima for the modified problem. The original and modified laws discovered by RF5 were

where the weight values were rounded off to the sec­ ond decimal place. Namely, the laws discovered by RF5 were almost perfect to Eq. (6) or (7) . Note that without preparing some appropriate prototype functions, exist­ ing numeric discovery methods cannot find such laws as described in Eq. (7) . This point is an important advan­ tage of RF5 over existing methods.

To evaluate RF5's noise tolerance, we corrupted each value of y calculated from Eq. (6) or (7) by adding noise generated according to a normal distribution with a mean of 0 and a standard deviation of 0.1. The other experimental conditions were exactly the same as be­ fore. Table 2 shows the results. The best MSE values were minimized when h = 3, while the best MDL values were minimized when h = 2: this indicates the correct number of hidden units was found again for the both problems. The original and modified laws discovered by RF5 were

1 0ur experiments were done on HP /9000/735 computer.

When the weight values were rounded off to the first decimal place, these laws become

Although some weight values were slightly different, laws almost equivalent to the true ones were found. This shows that RF5 is robust and noise tolerant to some degree.

5.2 Real data

For real data, we used three data sets supporting Ragen­ Rubens' law, Kepler's third law, and Boyle's law.2 In this experiment, since the number of examples for each data set was small, the number of hidden units was fixed at 1. Note that since we consider a constant term eo, the problem cannot be reduced to a simple regression problem. The trials were performed 10 times for each data set, and the weight values in the discovered laws were rounded off to the second decimal place. The other experimental conditions were exactly the same as before.

2 We obtained the sets of Kepler's data and Boyle's data from the UCI repository of machine learning databases.

Hagen-Rubens' law is a relation among the frequency v of incident light, the electrical conductivity <7 of metal, and the reflectance R. The original law is described as

The law discovered by RF5 is

Although the data has some outstanding outlier exam­ ples, the discovered law is very similar to the reference relation.

Kepler's third law is a relation between the distance r to the sun and the revolution period T of five planets. The original law is described as

The law discovered by RF5 is

An undesired constant term appeared here, but the law is very similar to the reference relation. The value of the constant term is small in the discovered law; as such we can perform another trial using the model where the constant term is omitted.

Boyle's law is a relation between the pressure p and volume V of a quantity of enclosed air. The original law is described as

The law discovered by RF5 is

5.3 Time series data

If time series data { x1 , · · · , Xt, · · ·} are generated accord­ ing to a nonlinear difference equation such as kxt = Xt-1(1- Xt - 1 ) (k is a constant and x0 E [0, 1]), then, it is known that we can observe a chaotic behavior. Since our target law defined by Eq. (1) is a generalized polyno­ mial function, it is expected that the law of chaotic time series can be identified by RF5. Hereafter, a nonlinear autoregressive function defined by

is considered, where z(·) denotes the right hand-side of Eq. (2) and T is the length of delayed inputs. For real time series data, we used a data set recorded from a Far-Infrared-Laser in a chaotic state [Weigend and Ger­ shenfeld, 1994] 3 , where we used a series of 1,000 training points (Data Set A) and Twas set to 8.

In the experiments, we changed the number of hidden units from 1 to 7 and performed 10 trials for each of them. The other experimental conditions were exactly the same as before. Figure 1 shows the best MSE values and MDL values. This indicates the best MDL values were minimized when h = 6. In this experiments, the law discovered by RF5 was not able to outperform the best result reported in the competition [Weigend and Gershenfeld, 1994], and we must go further to improve our preliminary result. Actually, for the successive 100 test points, the prediction performance of the discovered law was the 9th position among the 14 results submitted to the competition. However, the discovered law shows a good short-term prediction performance for the unknown data. Figure 2 shows a part of the predictions by the discovered law, where only first 8 points were given to the law. We can see that the predicted points are very close to the true ones.

6 Conclusion

To discover an underlying law from a set of numeric data, we have proposed a new connectionist method called RF5. After employing the connectionist problem formal­ ization, RF5 adopts a second-order learning algorithm BPQ for training and the MDL criterion for model se­ lection. Experiments showed that RF5 successfully dis­ covered underlying laws whose power values are not re­ stricted to integers, even if the data contained irrelevant variables or a small amount of noise. In the future, we plan to do further experiments to evaluate the proposed method using a wider variety of problems.

3 This data set was used in the Santa Fe Institute Time Series Prediction and Analysis Competition.