Numeric Law Discovery using Neural Networks
Kazumi Saito1 and Ryohei Nakano2
1

NTT R&D Management Department, 3-19-2, Nishi-shinjuku, Shinjuku-ku, Tokyo 163-19 Japansaito@rdh.ecl.ntt.co.jp
2
NTT Communication Science Laboratories, 2 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-02 Japannakano@cslab.kecl.ntt.co.jp

Abstract
This paper proposes a new connectionist approach to numeric law discovery; i.e., neural networks (law-candidates) aretrained by using a newly invented second-order learning algorithm based on a quasi-Newton method, called BPQ, and the
MDL criterion selects the most suitable from law-candidates.
The main advantage of our method over previous work ofsymbolic or connectionist approach is that it can eﬃciently
discover numeric laws whose power values are not restricted
to integers. Experiments showed that the proposed method
works well in discovering such 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. Such systems, forexample, 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 = kr 3/2 (k is a constant).
After the pioneering work of the BACON systems [5, 6], several methods [7, 3, 8, 16, 17] have been proposed. The basic
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 predeﬁned prototype functions. BACON and FAHRENHEIT [7] use trend
detectors to combine variables and employ a heuristic form ofdepth-ﬁrst search. ABACUS [3] creates a proportional graphand performs a modiﬁed beam search. In IDS [8], correlation analysis is applied, and a beam search is performed. The
E* algorithm [16] considers only bivariate functions. Also,
Sutton-Matheus’ algorithm [17] performs a regression, and thecorrelation between the squared error and the square of thevariable’s value [15] is used to combine variables.
These existing methods suﬀer from the following problems:ﬁrst, because combining two variables into a new one must
be done in order, a combinatorial explosion may occur whencomplex laws are sought for data consisting of a large number of variables, or the desired laws will be missed when someheuristic search parameters are inappropriate. Second, whensome powers appearing in a law are not restricted to integers,
the law may remain unknown unless some appropriate prototype functions 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 [8, 16].
We believe a connectionist approach has great potential tosolve the above problems. In order to directly learn a generalized polynomial term whose power values are not restricted tointegers, a computational unit called a product unit has been

proposed [2]; instead of calculating a weighted sum of 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 [13] to train networks containing these units [9].
Although some heuristic strategies such as multiple learningalgorithms have been proposed [9], their improvements overBP have been 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 approach called RF5for discovering numeric laws. Section 2 explains how neuralnetworks are used to discover a class of numeric 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, realdata, and time series data.

2. Law Discovery using Neural Nets
This section explains a connectionist problem formalization for numeric law discovery, ﬁrst proposed in [2]. Let{(x1 , y1 ), · · · , (xm , ym )} be a set of examples, where x t is ann-dimensional input vector and yt is a target value corresponding to xt . In this paper, a class of numeric laws expressed as
yt

=

c0 +

h


win
i1
ci xw
t1 · · · x tn

(1)

i=1

is considered, where each parameter ci or wij is an unknown
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 ﬁnite number of terms. Moreover, since thepower values are not restricted to integer, we can expect thatthe approximated polynomial function has a fewer number ofterms. Hereafter, (c0 , · · · , ch )T and (wi1 , · · · , win )T are expressed as c 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 we can assumexti > 0; then, Eq. (1) is equivalent to
yt

=

c0 +

h

i=1

ci exp

 n

j=1



wij ln(xtj ).

(2)

Equation (2) can be regarded as the feedforward computationof a three-layer neural network where the activation functionof each hidden unit is exp(s) = es . Here h, wi , and c denotethe 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, theoutput value of hidden
 unit i is described as vi (xt ; wi ) =
exp

n

j=1

wij ln(xtj ) and then the output value of the out-h

put unit is described as z(xt ; Φ) = c0 + i=1 ci vi (xt ; wi ).
The hidden units deﬁned by v(x; w) are called product units[2]. The discovery of numeric laws subject to Eq. (1) can thus
be deﬁned as the following learning problem in neural networks. That is, the problem is to ﬁnd the Φ that minimizes
the following objective function:1
(yt − z(xt ; Φ))2 .
2
m

f (Φ)

=

(3)

t=1

3. BPQ Algorithm
In our early experiments and as reported in earlier studies[9], the problem of minimizing Eq. (3) turned out to be quite
tough. Thus, in order to eﬃciently and constantly obtain goodresults, this paper employs a new second-order learning algorithm called BPQ [14]; by adopting a quasi-Newton method[4, 10] as a basic framework, the descent direction, ∆Φ, iscalculated on 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 ﬁrst-orderlearning 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 [4, 10] which converge more quickly byusing both gradient and curvature information, it is diﬃcult
to suitably scale up for large problems, and a large amountof computation is required for calculating the optimal steplength. BPQ can be reasonably scaled up by introducing astorage space parameter, and the computational 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 basic procedure forcalculating 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 thisapproximation is given by λ = −ζ  (0)/ζ  (0). Here, themethod for coping with the other cases is exactly the sameas described in [14]. For three-layer neural networks deﬁned by Eq. (3), we can eﬃciently calculate both ζ  (0) andζ  (0) as follows. By diﬀerentiatingm ζ(λ) and substituting
0 for λ, we obtain ζ  (0) = − t=1 (yt − z(xt ; Φ))z  (xt ; Φ),m ζ  (0) = t=1 (z  (xt ; Φ))2 − (yt − z(xt ; Φ))z  (xt ; Φ) . Now
is
that
the
derivativeofz(xt ; Φ)ddeﬁned as z (xt ; Φ) = dλz(xt ; Φ + λ∆Φ)|λ=0 , we obtain zt =
h
h



) and
zt = i=1 (2∆ci vit+ ci vit),∆c0 + i=1 (∆ci vit + ci vitn
where vi (xt ; wi ) = vi (xt ; wi ) j=1 ∆wij ln(xtj ), vi (xt ; wi ) =nvi (xt ; wi ) j=1 ∆wij ln(xtj ), 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 [11] or OSS [1], but BPQ worked the mosteﬃciently 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 thedata 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 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 [12] for this purpose.
The MDL ﬁtness value is deﬁned by
MDL=

0.5m log(MSE) + 0.5N log(m),(4)

where MSE represents the value of the mean squared error
deﬁned by1 (yt − z(xt ; Φ̂))2 .
m
m

MSE

=

(5)

t=1

Here, Φ̂ is a set of weights which minimizes Eq. (3), N is thenumber of parameters in Φ, and m is the number of examples.
Hereafter, our discovery method employing 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
artiﬁcial problem proposed by Sutton and Matheus [17] andour modiﬁed version.
The original problem is to restore a law described asy=

2 + 3x1 x2 + 4x3 x4 x5 ,(6)while a law for the modiﬁed problem is described asy=

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 , · · · , x 5 is randomly generated in the range of [0, 1], andthe 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 x 6 , · · · , x 9 is also randomly
generated in the range of [0, 1], each value of y is corrupted
by adding noise generated according to a normal distribution
with a mean of 0 and a standard deviation of 0.1, and thenumber 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 astandard deviation of 1; the initial values for the weights between the hidden and output units were set to 0, but thebias value at the output unit was initially set to the average
output value of all training examples. The iteration was terminated when the gradient vector was suﬃciently small, i.e.,∇f (Φ)2 /N < 10−8 , or the total processing time exceeded
100 seconds.

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

best
0.160
0.009
0.008
2.326
0.010
0.009

MSE valueavg.
s.d.
0.160 0.000
0.009 0.000
0.008 0.000
2.326 0.000
0.031 0.207
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 ofthem. Table 1 shows the basic statistics of MSE values, MDL
values, iterations, and processing times (sec.). 1 The best MSE
values were minimized when h = 3, while the best MDL values
were minimized when h = 2: this indicates the correct numberof hidden units was found again for the both problems. The
original and modiﬁed laws discovered by RF5 were
y
y

=
=

x−0.01x70.01 x−0.01
1.97 + 3.03x11.00 x20.97 x−0.01
45
8−0.03 −0.01 1.03 1.00 1.05 −0.01 −0.02+3.88x1x2
x3 x4 x5 x6
x72.01 + 3.00x−1.00
x23.001
1.02 0.50 −0.33 −0.01 −0.01+3.98x3 x4 x5
x6x9.

where the weight values were rounded oﬀ to the second decimal place. When the weight values were rounded oﬀ to theﬁrst decimal place, these laws become
y
y

=
=

2.0 + 3.0x11.0 x21.0 + 3.9x31.0 x41.0 x51.0
2.0 + 3.0x−1.0
x23.0 + 4.0x31.0 x40.5 x−0.3.
15Although some weight values were slightly diﬀerent, laws almost equivalent to the true ones were found. This shows that
RF5 is robust and noise tolerant to some degree. Note thatwithout 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.

5.2. Real data

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

iteration
avg. s.d.
68
4
93
9
784
82
90
12
228
134
753
127

time
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

The law discovered by RF5 isR

1.00 − 2.08ν 0.57 σ −0.57 .

=Although the data has some outstanding outlier examples, thediscovered law is very similar to the reference relation.
Kepler’s third law is a relation between the distance r to thesun and the revolution period T of ﬁve planets. The originallaw is described asT=

0.41r1.5 .

The law discovered by RF5 isT

0.41r1.50 + 0.19.

=

An undesired constant term appeared here, but the law isvery 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 Vof a quantity of enclosed air. The original law is described asV=

29.30/p.

The law discovered by RF5 isV=

29.05p−1.08 − 0.61.

5.3. Time series data

For real data, we used three data sets supporting HagenRubens’ law, Kepler’s third law, and Boyle’s law. 2 In thisexperiment, since the number of examples for each data set
was small, the number of hidden units was ﬁxed at 1. Notethat since we consider a constant term c 0 , 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 oﬀ to the second decimal
place. The other experimental conditions were exactly thesame as before.
Hagen-Rubens’ law is a relation among the frequency ν ofincident light, the electrical conductivity σ of metal, and thereﬂectance R. The original law is described as
R

=

1−2

  12νσ.

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

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

If time series data {x 1 , · · · , x t , · · ·} are generated according toa 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 expectedthat the law of chaotic time series can be identiﬁed by RF5.
Hereafter, a nonlinear autoregressive function deﬁned byxt=

z(xt−1 , · · · , x t−τ ; Φ),(8)is considered, where z(·) denotes the right hand-side of Eq. (2)and τ is the length of delayed inputs. For real time series data,
we used a data set recorded from a Far-Infrared-Laser in achaotic state [18] 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 of them. The
3 This data set was used in the Santa Fe Institute Time Series
Prediction and Analysis Competition.

50

MSE

MDL

40

2000

true value

predictied value

200

1600

30

1200
MDL

20

800

10

400

100
MSE

0

0
h=1 h=2 h=3 h=4 h=5 h=6 h=7
Figure 1: learning results

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 [18],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 forthe unknown data. Figure 2 shows a part of the predictions
by the discovered 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. Conclusion
To discover an underlying law from a set of numeric data, wehave 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 amountof noise. In the future, we plan to do further experiments toevaluate the proposed method using a wider variety of problems.

References[1] R. Battit. First- and second-order methods for learning
between steepest descent and newton’s method. NeuralComputation, 4(2):141–166, 1992.
[2] R. Durbin and D. Rumelhart. Product units: a computationally powerful and biologically plausible extension.
Neural Computation, 1(1):133–142, 1989.
[3] B.C. Falkenhainer and R.S. Michalski. Integrating quantitative and qualitative discovery in the abacus system.
In Y. Kodratoﬀ and R.S. Michalski, editors, Machinelearning: an artificial intelligence approach, Volume III,
pages 153–190. Morgan Kaufmann, San Mateo, CA, 1990.
[4] P.E. Gill, W. Murray, and M.H. Wright. Practical optimization. Academic Press, London, 1981.
[5] P. Langley. Bacon.1: a general discovery system. In
Proeeedings of the Seeond National Conferenee of theCanadian Soeiety for Computational Studies of Intelligenee, pages 173–180, 1978.

0
1200

1300

1400

Figure 2: prediction by discovered law[6] P. Langley, H.A. Simon, G. Bradshaw, and J. Zytkow.
Scientific discovery: computational explorations of thecreative process. MIT Press, Cambridge, MA, 1987.
[7] P. Langley and J. Zytkow. Data-driven approaches to
empirical discovery. Artificial Intelligence, 40:283–312,
1989.
[8] P. Langley and J. Zytkow. A robust approach to numericdiscovery. In Proc. seventh International Machine Learning Conference, pages 411–418, Austin, Texas, 1990.
[9] L.R. Leerink, C.L. Giles, B.G. Horne, and M.A. 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.
[10] D.G. Luenberger. Linear and nonlinear programming.
Addison-Wesley, Reading, MA, 1984.
[11] M.F. Møller. Supervised learning on large redundant
training sets. International Journal of Neural Systems,
4(1):15–25, 1993.
[12] J. Rissanen. Stochastic complexity in statistical inquiry.
World Scientiﬁc, 1989.
[13] D.E. Rumelhart, G.E. Hinton, and R.J. Williams. Learning internal representations by error propagation. In D.E.
Rumelhart and J.L. McClelland, editors, Parallel Distributed Processing, pages 318–362. MIT Press, 1986.
[14] K. Saito and R. Nakano. Partial BFGS update and efﬁcient step-length calculation for three-layer neural networks. Neural Computation, 9(1):239–257, 1997.
[15] T.D. 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.
[16] C. Schaﬀer. Bivariate scientiﬁc function ﬁnding ina sampled, real-data testbed.
Machine Learning,
12(1/2/3):167–183, 1993.
[17] R.S. Sutton and C.J. Matheus. Learning polynomial
functions by feature construction. In Proceedings of theEighth International Machine Learning Workshop, pages
208–212, Evanston, IL, 1991.
[18] A.S. Weigend and N.A. Gershenfeld. Time series prediction: forcasting future and understanding the past.
Addison-Wesley, 1994.