Law Discovery using Neural Networks
Kazumi Saito and Ryohei Nakano
NTT Communication Science Laboratories
2 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-02 Japan

Abstract
This paper proposes a new connectionist approach to numeric law discovery; i.e., neuralnetworks (law-candidates) are trained by using
a newly invented second-order learning algorithm based on a quasi-Newton method, called
BPQ, and the Minimum Description Length
criterion selects the most suitable from lawcandidates. The main advantage of our method
over previous work of symbolic or connectionist approach is that it can eﬃciently discover
numeric laws whose power values are not restricted to integers. Experiments showed thatthe proposed method works well in discoveringsuch laws even from data containing irrelevant
variables or a small amount of noise.

1

Introduction

The discovery of a numeric law from a set of data is
the central part of scientiﬁc discovery systems. Suchsystems, for example, can detect a relationship between
the distance r to the sun and the revolution period T ofﬁve planets known to Kepler’s third law T = kr3/2 (k isa constant).
After the pioneering work of the BACON systems[Langley, 1978; Langley et al., 1987], several methods[Langley and Zytkow, 1989; Falkenhainer and Michalski, 1990; Langley and Zytkow, 1990; Schaﬀer, 1993;
Sutton and Matheus, 1991] have been proposed. The basic search strategy employed by these methods is muchthe same: two variables are recursively combined into a
new variable by using multiplication, division, or somepredeﬁned prototype functions. BACON and FAHRENHEIT [Langley and Zytkow, 1989] use trend detectors
to combine variables and employ a heuristic form ofdepth-ﬁrst search. ABACUS [Falkenhainer and Michalski, 1990] creates a proportional graph and performs a
modiﬁed beam search. In IDS [Langley and Zytkow,
1990], correlation analysis is applied, and a beam searchis performed. The E* algorithm [Schaﬀer, 1993] considers only bivariate functions. Also, Sutton-Matheus’
algorithm [Sutton and Matheus, 1991] performs a regression, and the correlation between the squared error
and the square of the variable’s value [Sanger, 1991] isused to combine variables.
These existing methods suﬀer from the following problems: ﬁrst, because combining two variables into a newone must be done in order, a combinatorial explosion
may occur when complex laws are sought for data consisting of a large number of variables, or the desiredlaws will be missed when some heuristic search parameters are inappropriate. Second, when some powers appearing in a law are not restricted to integers, the law
may remain unknown unless some appropriate prototype
functions such as r3/2 are prepared in advance. However, a priori information is rarely available. Third,these methods are often criticized for their lack of robustness; noise tolerance is deﬁnitely required since real
observed data contain noise [Langley and Zytkow, 1990;
Schaﬀer, 1993].
We believe a connectionist approach has great potential to solve the above problems. In order to directly
learn a generalized polynomial term whose power values are not restricted to integers, a computational unit
called a product unit has been proposed [Durbin and
Rumelhart, 1989]; instead of calculating a weighted sumof input values, this unit calculates a weighted product, where each input value is raised to a power determined by a variable weight. However, serious diﬃculties
have been reported when using standard BP [Rumelhart 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 havebeen less than remarkable. Moreover, these earlier studies dealt only with binary data and did not speciﬁcallyaddress the problem of numeric law discovery.
In this paper, we propose a connectionist approachcalled RF5 for discovering numeric laws. Section 2 ex-

plains how neural networks are used to discover a class ofnumeric laws. Section 3 describes a second-order learning algorithm called BPQ which trains the neural networks described in Section 2. Section 4 explains a criterion for selecting the most suitable candidate out of the
trained ones. Section 5 evaluates the proposed method
RF5 by doing experiments using artiﬁcial data, real data,
and time series data.

2

Law Discovery using Neural Nets

This section explains a connectionist problem formalization for numeric law discovery, ﬁrst proposed in [Durbin
and Rumelhart, 1989]. Let {(x1 , y1 ), · · · , (xm , ym )} bea 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 asyt= c0 +

h


win
i1
ci xw
t1 · · · xtn

(1)

i=1

is considered, where each parameter ci or wij is an unknown real number and h is an unknown integer. If atarget law consists of periodic function or discontinuousfunction, we cannot exactly discover it when Eq. (1) isassumed. However, if the range of each input variable
is bounded, such a law can be closely approximated toa multivariate polynomial function with a ﬁnite numberof terms. Moreover, since the power values are not restricted to integer, we can expect that the approximatedpolynomial function has a fewer number of terms. Hereafter, (c0 , · · · , ch )T and (wi1 , · · · , win )T are expressed asc and wi , respectively, where aT means the transposed
vector of a. In addition, a vector consisting of all parameters, (cT , w1T , · · · , whT )T is simply expressed as Φ, andN (= nh + h + 1) denotes the dimension of Φ.
By adding an adequate value to each element of the
input vectors, if necessary, without losing generality wecan assume xti > 0; then, Eq. (1) is equivalent to


h
n


yt = c0 +ci exp 
wij ln(xtj ) .
(2)
i=1

j=1

Equation (2) can be regarded as the feedforward computation of a three-layer neural network where the activation function of each hidden unit is exp(s) = es . Hereh, wi , and c denote the number of hidden units, theweights between the input units and hidden unit i, andthe weights between the hidden units and the output
unit, respectively. Hereafter, the output value of hidden
unit i is described asn
vi (xt ; wi ) = exp 
wij ln(xtj )
j=1

and then the output value of the output unit is describedasz(xt ; Φ) = c0 +h

ci vi (xt ; wi ).

i=1

The hidden units deﬁned by v(x; w) are called productunits [Durbin and Rumelhart, 1989]. The discovery ofnumeric laws subject to Eq. (1) can thus be deﬁned as thefollowing learning problem in neural networks. That is,
the problem is to ﬁnd the Φ that minimizes the following
objective function:mf (Φ) =

31(yt − z(xt ; Φ))2 .
2 t=1

(3)

BPQ Algorithm

In our early experiments and as reported in earlier studies [Leerink et al., 1995], the problem of minimizing
Eq. (3) turned out to be quite tough. Thus, in order toeﬃciently and constantly obtain good results, this paper
employs a new second-order learning algorithm called
BPQ [Saito and Nakano, 1997]; by adopting a quasiNewton method [Gill et al., 1981; Luenberger, 1984] as abasic framework, the descent direction, ∆Φ, is calculatedon the basis of a partial BFGS update and a reasonably
accurate step-length, λ, is eﬃciently calculated as theminimal point of a second-order approximation. In ﬁrstorder learning algorithms which calculate the search direction as the gradient direction, a large number of iterations are often required until convergence. On the other
hand, in existing second-order methods [Gill et al., 1981;
Luenberger, 1984] which converge more quickly by using
both gradient and curvature information, it is diﬃcult tosuitably scale up for large problems, and a large amountof computation is required for calculating the optimal
step-length. BPQ can be reasonably scaled up by introducing a storage space parameter, and the computational complexity for calculating the optimal step-length
is reasonably small, almost equivalent to that of gradientvector evaluation.
For the problem of minimizing Eq. (3), the partial
BFGS update can be directly applied, while the basic procedure for calculating the step-length λ must be
slightly modiﬁed. In the step-length calculation, since λis the only variable, we can express f (Φ + λ∆Φ) simply as ζ(λ). Its second-order Taylor approximation is
given as ζ(0) + ζ  (0)λ + 12 ζ  (0)λ2 . When ζ  (0) < 0 andζ  (0) > 0, the minimal point of this approximation is
given byζ  (0)λ = −  .
ζ (0)Here, the method for coping with the other cases is exactly the same as described in [Saito and Nakano, 1997].

For three-layer neural networks deﬁned by Eq. (3), we
can eﬃciently calculate both ζ  (0) and ζ  (0) as follows.
By diﬀerentiating ζ(λ) and substituting 0 for λ, we obtain

where MSE represents the value of the mean squared
error deﬁned by
m

MSE =

m

ζ (0) = − (yt −z(xt ; Φ))z  (xt ; Φ),


t=1

m


 
ζ  (0) =(z (xt ; Φ))2 −(yt −z(xt ; Φ))z  (xt ; Φ) .
t=1

Now that the derivative of z(xt ; Φ) is deﬁned asz  (xt ; Φ) =dz(xt ; Φ + λ∆Φ)|λ=0 ,dλz  (xt ; Φ) = ∆c0 +

h


5.1
(∆ci vi (xt ; wi )+ci vi (xt ; wi )) ,i=1z  (xt ; Φ) =h(2∆ci vi (xt ; wi )+civi (xt ; wi )) ,Evaluation by Experiments
Artiﬁcial data

The rule discovery method, RF5, was evaluated by using
an artiﬁcial problem proposed by Sutton and Matheus[Sutton and Matheus, 1991] and our modiﬁed version.
The original problem is to restore a law described asi=1vi (xt ; wi ) = vi (xt ; wi )

n

j=1
n


(6)

while a law for the modiﬁed problem is described as∆wij ln(xtj ),y∆wij ln(xtj ),j=1

and ∆ci and ∆wij are modiﬁcation values corresponding
to ci and wij , respectively.
Incidentally, we can employ other second-order learning algorithms such as SCG [Møller, 1993] or OSS [Battiti, 1992], but BPQ worked the most eﬃciently among
them in our own experience.

4= 2 + 3x1 x2 + 4x3 x4 x5 ,

y

where
vi (xt ; wi ) = vi (xt ; wi )(5)Here, Φ̂ is a set of weights which minimizes Eq. (3), N
is the number of parameters in Φ, and m is the numberof examples. Hereafter, our discovery method employing the connectionist problem formalization, the BPQ
algorithm and the MDL criterion, is called RF5 (Ruleextraction from Facts version 5).

5

we obtain1 (yt − z(xt ; Φ̂))2 .
m t=1

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 lawcandidate 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 thenumber of hidden units. In this paper, by assumingthat the target output values are corrupted by Gaussian
noise with a mean of 0 and an unknown standard deviation of σ, ﬁnding an adequate number of hidden units
is formalized as a model selection problem of the maximum likelihood estimation problem. Thus, we adopt
the MDL (Minimum Description Length) criterion [Rissanen, 1989] for this purpose. The MDL ﬁtness value isdeﬁned by
MDL = 0.5m log(MSE) + 0.5N log(m),(4)

1/2 −1/3

3= 2 + 3x−1
1 x2 + 4x3 x4 x5.

(7)Each example is generated as follows: each value of variables x1 , · · · , x5 is randomly generated in the range of[0, 1], and the corresponding value of y is calculated using Eq. (6) or (7). In these problem, the total number of variables is 9 (n = 9); each value of irrelevant
variables x6 , · · · , x9 is also randomly generated in therange 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 amean of 0 and a standard deviation of 1; the initial values 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 suﬃciently small, i.e.,m1 (yt − z(xt ; Φ))2
m t=1

< 10−8 ,(8)the gradient vector was suﬃciently small, i.e.,1∇f (Φ)2
N< 10−8 ,or the total processing time exceeded 100 seconds.

(9)

Table 1: Learning statistics (artiﬁcial data)

original problem

modiﬁed problem

hidden
unit
h=1
h=2
h=3
h=1
h=2
h=3

MSE value
best
avg.
s.d.
0.126 0.130 0.000
0.000 0.000 0.000
0.000 0.000 0.000
1.317 1.321 0.000
0.000 0.033 0.172
0.000 0.000 0.000

MDL valuebest
avg.
-177.9-178-1786.4 -1786-1759.9 -1717
56.7
57-1786.4 -1731-1760.0 -1727

s.d.
0.0
0.0
64
0.0
314
47

iteration
avg. s.d.
71
1.0
81
8.4
130 2970
5116 27240 104time
avg. s.d.
0.32 0.01
0.68 0.07
1.65 0.36
0.32 0.02
0.97 0.23
3.02 1.31

iteration
avg. s.d.
68
4939
784 8290
12
228 134
753 127time
avg. s.d.
0.31 0.02
0.77 0.07
9.80 1.02
0.41 0.05
1.90 1.11
9.41 1.58

Table 2: Learning statistics (noisy artiﬁcial data)

original problem

modiﬁed problem

hidden
unit
h=1
h=2
h=3
h=1
h=2
h=3

MSE value
best
avg.
s.d.
0.160 0.160 0.000
0.009 0.009 0.000
0.008 0.008 0.000
2.326 2.326 0.000
0.010 0.031 0.207
0.009 0.009 0.000

In the experiments, we changed the number of hidden
units from 1 to 3 (h = 1, 2, 3) and performed 100 trials 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 numberof hidden units, h = 2, was successfully found because
the best MDL value was minimized. When h = 2, all 100
trials converged to the global minimum for the original
problem; however, a few of the 100 trials converged toundesirable local minima for the modiﬁed problem. The
original and modiﬁed laws discovered by RF5 wereyy

= 2.00 + 3.00x1.00
x1.00+ 4.00x1.00
x1.00
x1.00
1
2
3
4
5
−1.00 3.00
1.00 0.50 −0.33= 2.00 + 3.00x1
x2 + 4.00x3 x4 x5

where the weight values were rounded oﬀ to the second decimal place. Namely, the laws discovered by RF5
were almost perfect to Eq. (6) or (7). Note that without
preparing some appropriate prototype functions, existing numeric discovery methods cannot ﬁnd such laws as
described in Eq. (7). This point is an important advantage of RF5 over existing methods.
To evaluate RF5’s noise tolerance, we corrupted eachvalue of y calculated from Eq. (6) or (7) by adding
noise generated according to a normal distribution witha mean of 0 and a standard deviation of 0.1. The otherexperimental conditions were exactly the same as before. Table 2 shows the results. The best MSE values
were minimized when h = 3, while the best MDL values
1

Our experiments were done on HP/9000/735 computer.

MDL value
best
avg.
s.d.
-154.0 -154.0 0.0-416.0 -416.0 0.0-405.5 -405.3 0.8
113.6
113.6
0.0
-403.9 -398.6 53-388.9 -384.9 1.2

were minimized when h = 2: this indicates the correct
number of hidden units was found again for the bothproblems. The original and modiﬁed laws discovered by
RF5 wereyy= 1.97 + 3.03x1.00
x0.97
x−0.01x−0.01x0.01
x−0.0112
7458+3.88x−0.03
x−0.01
x1.03
x1.00
x1.05
x−0.01x−0.02
3
4
5
1
2
6
7
3.00
= 2.01 + 3.00x−1.00
x2
10.50 −0.33 −0.01 −0.01+3.98x1.02
x
x
x6
x9
.
3
4
5

When the weight values were rounded oﬀ to the ﬁrst
decimal place, these laws becomeyy

1.0
1.0 1.0 1.0= 2.0 + 3.0x1.0
1 x2 + 3.9x3 x4 x5−1.0 3.0
0.5 −0.3= 2.0 + 3.0x1 x2 + 4.0x1.0.
3 x4 x5

Although some weight values were slightly diﬀerent, laws
almost equivalent to the true ones were found. Thisshows that RF5 is robust and noise tolerant to somedegree.

5.2

Real data

For real data, we used three data sets supporting HagenRubens’ law, Kepler’s third law, and Boyle’s law.2 In
this experiment, since the number of examples for eachdata set was small, the number of hidden units was ﬁxed
at 1. Note that since we consider a constant term c0 ,the problem cannot be reduced to a simple regression
2

We obtained the sets of Kepler’s data and Boyle’s datafrom the UCI repository of machine learning databases.

problem. The trials were performed 10 times for eachdata set, and the weight values in the discovered laws
were rounded oﬀ to the second decimal place. The otherexperimental conditions were exactly the same as before.
Hagen-Rubens’ law is a relation among the frequency
ν of incident light, the electrical conductivity σ of metal,
and the reﬂectance R. The original law is described as
= 1−2

R

ν

1
2

σ

50

20
10

= 0.41r1.5 .

= 0.41r1.50 + 0.19.

= 29.30/p.

= 29.05p−1.08 − 0.61.

If time series data {x1 , · · · , xt , · · ·} are generated according to a nonlinear diﬀerence equation such as kxt =xt−1 (1 − xt−1 ) (k is a constant and x0 ∈ [0, 1]), then, it
is known that we can observe a chaotic behavior. Since
our target law deﬁned by Eq. (1) is a generalized polynomial function, it is expected that the law of chaotic timeseries can be identiﬁed by RF5. Hereafter, a nonlinear
autoregressive function deﬁned by= z(xt−1 , · · · , xt−τ ; Φ),0time series data, we used a data set recorded from a
Far-Infrared-Laser in a chaotic state [Weigend and Gershenfeld, 1994] 3 , where we used a series of 1,000 training
points (Data Set A) and τ was set to 8.
In the experiments, we changed the number of hidden
units from 1 to 7 and performed 10 trials for each ofthem. 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, thelaw discovered by RF5 was not able to outperform thebest result reported in the competition [Weigend andGershenfeld, 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 agood short-term prediction performance for the unknown
data. Figure 2 shows a part of the predictions by thediscovered law, where only ﬁrst 8 points were given to
the law. We can see that the predicted points are very
close to the true ones.

6

Time series data

xt

400

Figure 1: learning results

The law discovered by RF5 is

5.3

MSE
h=1 h=2 h=3 h=4 h=5 h=6 h=7

An undesired constant term appeared here, but the law
is very similar to the reference relation. The value ofthe constant term is small in the discovered law; as such
we can perform another trial using the model where theconstant term is omitted.
Boyle’s law is a relation between the pressure p andvolume V of a quantity of enclosed air. The original law
is described as

V

800

0

The law discovered by RF5 is

V

1200
MDL

.

Although the data has some outstanding outlier examples, the discovered law is very similar to the reference
relation.
Kepler’s third law is a relation between the distancer to the sun and the revolution period T of ﬁve planets.
The original law is described as

T

2000
1600

30

= 1.00 − 2.08ν0.57 σ−0.57 .

T

MDL

40

The law discovered by RF5 is
R

MSE

(10)

is considered, where z(·) denotes the right hand-side ofEq. (2) and τ is the length of delayed inputs. For real

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 formalization, RF5 adopts a second-order learning algorithm
BPQ for training and the MDL criterion for model selection. Experiments showed that RF5 successfully discovered underlying laws whose power values are not restricted to integers, even if the data contained irrelevant
variables or a small amount of noise. In the future, weplan to do further experiments to evaluate the proposed
method using a wider variety of problems.
3This data set was used in the Santa Fe Institute Time
Series Prediction and Analysis Competition.

true value

predictied value

200

100

0
1200

[Leerink et al., 1995] L. Leerink, C. Giles, B. Horne, andM. Jabri. Learning with product units. In G. Tesauro,
D.S. Touretzky, and T.K. Lee, editors, Advances in
Neural Information Processing Systems 7. MIT Press,
Cambridge, MA, 1995.
[Luenberger, 1984] D. Luenberger. Linear and nonlinear
programming. Addison-Wesley, Reading, MA, 1984.
[Møller, 1993] M. Møller. Supervised learning on large
redundant training sets. International Journal of Neural Systems, 4(1):15–25, 1993.

1300

1400

Figure 2: prediction by discovered law

References[Battiti, 1992] R. Battiti. First- and second-order methods for learning between steepest descent and newton’s method. Neural Computation, 4(2):141–166,
1992.
[Durbin and Rumelhart, 1989]
R. Durbin and D.
Rumelhart. Product units: a computationally powerful and biologically plausible extension. Neural Computation, 1(1):133–142, 1989.
[Falkenhainer and Michalski, 1990] B. Falkenhainer and
R. Michalski. Integrating quantitative and qualitative discovery in the abacus system. In Y. Kodratoﬀ
and R. Michalski, editors, Machine learning: an artificial intelligence approach, Volume III, pages 153–190.
Morgan Kaufmann, San Mateo, CA, 1990.
[Gill et al., 1981] P. Gill, W. Murray, and M. Wrigh.
Practical optimization. Academic Press, London,
1981.
[Langley and Zytkow, 1989] P. Langley and J. Zytkow.
Data-driven approaches to empirical discovery. Artificial Intelligence, 40:283–312, 1989.
[Langley and Zytkow, 1990] P. Langley and J. Zytkow.
A robust approach to numeric discovery. In Proc.
seventh International Machine Learning Conference,
pages 411–418, Austin, Texas, 1990.
[Langley et al., 1987] P. Langley, H. Simon, G. Bradshaw, and J. Zytkow. Scientific discovery: computational explorations of the creative process. MIT Press,
Cambridge, MA, 1987.
[Langley, 1978] P. Langley. Bacon.1: a general discovery
system. In Proeeedings of the Seeond National Conferenee of the Canadian Soeiety for Computational Studies of Intelligenee, pages 173–180, 1978.

[Rissanen, 1989] J. Rissanen. Stochastic complexity instatistical inquiry. World Scientiﬁc, Singapore, 1989.
[Rumelhart et al., 1986] D. Rumelhart, G. Hinton, and
R. Williams. Learning internal representations by error propagation. In D.E. Rumelhart and J.L. McClelland, editors, Parallel Distributed Processing, pages
318–362. MIT Press, Cambridge, MA, 1986.
[Saito and Nakano, 1997] K. Saito and R. Nakano. Partial BFGS update and eﬃcient step-length calculation
for three-layer neural networks. Neural Computation,
9(1):239–257, 1997.
[Sanger, 1991] T. Sanger. Basis-function trees as a generalization of local variable selection method for function approximation. In D.S. Touretzky, editor, Neural Information Processing Systems 3, pages 707–713.
Morgan Kaufmann, San Mateo, CA, 1991.
[Schaﬀer, 1993] C. Schaﬀer. Bivariate scientiﬁc function ﬁnding in a sampled, real-data testbed. Machine
Learning, 12(1/2/3):167–183, 1993.
[Sutton and Matheus, 1991] R. Sutton and C. Matheus.
Learning polynomial functions by feature construction. In Proceedings of the Eighth International Machine Learning Workshop, pages 208–212, Evanston,IL, 1991.
[Weigend and Gershenfeld, 1994]
A. Weigend and
N. Gershenfeld. Time series prediction: Forecasting
the futire and undestanding past. Addison-Wesley,
Reading, MA, 1994.