A divide and conquer method for symbolic regression ✩
Changtong Luoa,∗, Chen Chena,b , Zonglin Jianga,b

arXiv:1705.08061v2 [] 27 Jun 2017

a

State Key Laboratory of High Temperature Gas Dynamics, Institute of Mechanics,
Chinese Academy of Sciences, Beijing 100190, China
b
School of Engineering Sciences, University of Chinese Academy of Sciences,
Beijing, 100049, China

Abstract
Symbolic regression aims to find a function that best explains the relationship between independent variables and the objective value based on a given
set of sample data. Genetic programming (GP) is usually considered as an
appropriate method for the problem since it can optimize functional structure and coefficients simultaneously. However, the convergence speed of GP
might be too slow for large scale problems that involve a large number of
variables. Fortunately, in many applications, the target function is separable
or partially separable. This feature motivated us to develop a new method,
divide and conquer (D&C), for symbolic regression, in which the target function is divided into a number of sub-functions and the sub-functions are then
determined by any of a GP algorithm. The separability is probed by a new
proposed technique, Bi-Correlation test (BiCT). D&C powered GP has been
tested on some real-world applications, and the study shows that D&C can
help GP to get the target function much more rapidly.
Keywords: Mathematical modeling, Genetic programming, Symbolic
regression, Artificial intelligence, Divide and conquer

This work has been supported by the National Natural Science Foundation of China
(Grant No. 11532014).
∗
Corresponding author
Email addresses: luo@imech.ac.cn (Changtong Luo), chenchen@imech.ac.cn
(Chen Chen), zljiang@imech.ac.cn (Zonglin Jiang)
✩

Preprint submitted to Expert Systems with Applications

June 28, 2017

1. Introduction
Symbolic regression (SR) is a data-driven modeling method which aims to
find a function that best explains the relationship between independent variables and the objective value based on a given set of sample data (Schmidt and Lipson,
2009). Genetic programming (GP) is usually considered as a good candidate for SR since it does not impose a priori assumptions and can optimize
function structure and coefficients simultaneously. However, the convergence
speed of GP might be too slow for large scale problems that involve a large
number of variables.
Many efforts have been devoted trying to improve the original GP (Koza,
2008) in several ways. Some works suggest replacing its tree-based method,
with an integer string (Grammar Evolution) (O’Neill and Ryan, 2001), or
a parse matrix (Parse-Matrix Evolution) (Luo and Zhang, 2012). These
techniques can simplify the coding and decoding process but help little on
improving the convergence speed. Some other works suggest confining its
search space to generalized linear space, for example, Fast Function eXtraction (McConaghy, 2011), and Elite Bases Regression (Chen et al., 2017).
These techniques can accelerate the convergence speed of GP, even by orders
of magnitude. However, the speed is gained at the sacrifice of losing the
generality, that is, the result might be only a linear approximation of the
target function.
Fortunately, in many applications, the target function is separable or
partially separable (see section 2 for definitions). For example, in gas dynamics (Anderson, 2006), the heat flux coefficient St of a flat plate could be
formulated as
p
p
St = 2.274 sin(θ) cos(θ)/ Rex ,
(1)
and the heat flux qs at the stagnation point of a sphere as
p
qs = 1.83 × 10−4 v 3 ρ/R(1 − hw /hs ).

(2)

In equation (1), the two independent variables, θ and Rex , are both separable.
In equation (2), the first three variables v, ρ, and R are all separable, and
the last two variables, hw and hs , are not separable, but their combination
(hw , hs ) is separable. The function in equation (2) is considered partially
separable in this paper.
The feature of separability will be used in this paper to accelerate the
optimization process of symbolic regression. Some basic concepts on function separability are defined in Section 2. Section 3 describes the overall
2

work flow of the proposed method, divide and conquer. Section 4 presents
a special technique, bi-correlation test (BiCT), to determine the separability
of a function. Numerical results are given in Section 5, and the concluding
remarks are drawn in Section 6.
2. Basic concepts
The proposed method in this paper is based on a new concept referred to
as partial separability. It has something in common with existing separability
definitions such as reference (Berenguel et al., 2013) and (d’Avezac et al.,
2011), but is not exactly the same. To make it clear and easy to understand,
we begin with some illustrative examples. The functions as follows could all
be regarded as partially separable:
z = 0.8 + 0.6 ∗ ( u2 + cos(u) ) + sin(v + w) ∗ (v − w) ;

(3)

z = 0.8 + 0.6 ∗ ( u2 + cos(u) ) − sin(v + w) ∗ (v − w) ;

(4)

z = 0.8 + 0.6 ∗ ( u2 + cos(u) ) ∗ sin(v + w) ∗ (v − w) ;

(5)

where the boxed frames are used to indicate sub-functions, u is separable with
respect to z, while v and w themselves are not separable, but their combination (v, w) is considered separable. A simple example of non-separable
function is f (x) = sin(x1 + x2 + x3 ).
More precisely, the separability could be defined as follows.
Definition 1. A scalar function with n continuous variables f (x) (f : Rn 7→
R, x ∈ Rn ) is said to be partially separable if and only if it can be rewritten
as
f (x) = c0 ⊗1 ϕ1 (I1 x) ⊗2 ϕ2 (I2 x) ⊗3 · · · ⊗m ϕm (Im x)
(6)
where the binary operator ⊗i could be plus (+), minus (−), times(×). Ii is a
sub-matrix of the identity matrix, and Ii ∈ Rni ×n . The set { I1 , I2 , · · · , Im }
m
P
is a partition of the identity matrix I ∈ Rn×n ,
ni = n. The sub-function
i=1

ϕi is a scalar function such that ϕi : Rni 7→ R. Otherwise the function is
said to be non-separable.

3

In this definition, the binary operator, division (/), is not included in ⊗
for simplicity. However, this does not affect much of its generality, since the
sub-functions are not preset, and can be transformed as ϕ̃i (·)= 1/ϕi (·) if only
ϕi (·) 6= 0.
A special case is that all variables are separable, which could be defined
as follows.
Definition 2. A scalar function with n continuous variables f (x) (f : Rn 7→
R, x ∈ Rn ) is said to be completely separable if and only if it can be rewritten
as equation (6) and ni = 1 for all i = 1, 2, · · · , m.
3. Divide and conquer method
As above mentioned, many practical problems have the feature of separability. To make use of this feature to accelerate the optimization process of
genetic programming, a new method, divide and conquer (D&C), is proposed.
It works as follows.
First, a separability detection process is carried out to find out whether
the concerned problem is separable (at least partial separable) or not. The
variables are identified one by one, and then their combinations. Once it (or
the variable combination) is identified as separable, a sub-function ϕi (xi ) (or
ϕi (Ii x) for variable combinations) will be assigned. In this way, the structure
of target function f (x) could be divided into a set of sub-functions based on
the separability information: ϕi (Ii x), i = 1, 2, · · · , m.
Then, the sub-functions ϕi (Ii x) (i = 1, 2, · · · , m) are optimized and
determined one by one, using any of genetic programming algorithms, including classical GP, Grammatical Evolution (O’Neill and Ryan, 2001) and
Parse-Matrix Evolution (Luo and Zhang, 2012). When optimizing one subfunction, variables not involved ((I − Ii )x) are fixed as constants. That is,
only a small number of variables (Ii x, which is only a subset of {x1 , x2 , · · · , xn })
need be considered. This means the sub-function determination should be
much easier than evolving the target function f (x) directly.
For example, in Equation (3), Equation (4), or Equation (5), the
sub-function ϕ1 (u) = u2 + cos(u) (or ϕ2 (v, w) = sin(v + w) ∗ (v − w)) has
less number of variables and complexity than the original function.
Thus, optimizing them one by one is much easier for GP.
Finally, these sub-functions are properly combined to form the target
function, which is referred to as a function recover process. This process
4

Separability Detection
Function Division
Sub−Function Determination
Function Recover

Figure 1: Work flow of divide and conquer for symbolic regression

is rather simple, and all traditional regression algorithms are qualified to
accomplish this mission.
The work flow of D&C could be described in Figure 1.
4. Bi-correlation test
4.1. Description
The main idea of the proposed divide and conquer (D&C) method is to
make use of the separability feature to simplify the search process. Therefore,
the most important and fundamental step (the separability detection process
in Fig. 1) is to determine whether the concerned problem is separable (at
least partial separable) or not. To fulfill this task, a special technique, Bicorrelation test (BiCT), is provided in this section.
Consider independent variables x1 , x2 , · · · , xn and the dependent f as
n+1 random variables, and the known data as sample points. Recall that
the linear relation and correlation coefficient of two random variables has the
following relations.
Lemma 1. The two random variables ξ and η are linearly related with correlation coefficient 1 (i.e., ρξη = 1) if and only if there exists two constants
a, b (b 6= 0), such that P {η = a + bξ} = 1.
The correlation coefficient ρ could be estimated by the sample correlation

5

coefficient r, which is defined as follows.
N
¯ (ηi − η̄)
1 X (ξi − ξ)
r=
·
n − 1 i=1
σξ
ση

P
where N is the number of observations in the sample set,
is the summation
¯
symbol, ξi is the ξ value for observation i, ξ is the sample mean of ξ, ηi is
the η value for observation i, η̄ is the sample mean of η, σξ is the sample
standard deviation of ξ, and ση is the sample standard deviation of η.
Only continuous model functions are considered in this paper. As a result,
the conclusion of Lemma 1 could be simplified as follows.
• The two random variables f A and f B are linearly related (f B = a +
bf A ) if and only if the sample correlation coefficient r = 1 for any given
sample set.
Studies shows that the functional separability defined in the above section (See equation 6) could be observed with random sampling and linear
correlation techniques.
Without the loss of generality, a simple function with three variables
(f (x) = f (x1 , x2 , x3 ), xi ∈ [ai , bi ], i = 1, 2, 3) is considered to illustrate the
implementation of the bi-correlation test. To find out whether the first variable x1 is separable, two correlation tests are needed.
The first correlation test is carried out as follows. A set of random sample points in [a1 , b1 ] are generated, then these points are extended to three
dimensional space with the rest variables (x2 and x3 ) fixed to a point A. We
get a vector of function values f (A) = f (x1 , A) = (f1A , f2A , · · · , fNA ), where
N is the number of sample points. Then these points are extended to three
dimensional space with fixed x2 and x3 to another point B, We get another
vector f (B) = f (x1 , B) = (f1B , f2B , · · · , fNB ). It is obviously that the two
vectors f (A) and f (B) will be linearly related if x1 is separable. However, it
is easy to show that this linear relation could NOT ensure its separability.
Then it comes to the second correlation test. Another set of random
sample points in [a2 , b2 ]×[a3 , b3 ] are generated, then these points are extended
to three dimensional space with the rest variable(s) (x1 in this case) fixed to
a point C, and get a vector f (C) = f (C, x2 , x3 ). Similarly, another vector
f (D) = f (D, x2 , x3 ) is obtained. Again, the two vectors f (C) and f (D) needs
to be linearly related to ensure the separability of x1 .
6

4.2. Proposition
Without the loss of generality, suppose we have a scalar function f (x)
with n continuous variables (f : Rn 7→ R, x ∈ Ω ⊂ Rn , and Ω = [a1 , b1 ] ×
[a2 , b2 ] × · · · × [an , bn ]), and need to find out whether the first m variable
combination (x1 , x2 , · · · , xm ) are separable. Let the matrix X1 be a set of N
random sample points from the subset [a1 , b1 ] × [a2 , b2 ] × · · · × [am , bm ] ⊂ Rm ,
and


(1)
(1)
(1)
x1
x2
· · · xm
 (2)
(2)
(2) 
 x
x2
· · · xm 
X1 =  1
.
 ··· ··· ··· ··· 
(N )
(N )
(N )
x1
x2
· · · xm

The rest variables xm+1 , xm+2 , · · · , xn are fixed to two given points A and
B in the subset [am+1 , bm+1 ] × [am+2 , bm+2 ] × · · · × [an , bn ] ⊂ Rn−m , i.e.,
(A)
(B)
(B)
(B)
xA = (xm+1 (A) , xm+2 (A) , 
··· ,x
n ), xB = (xm+1 , xm+2 , · · · , xn ).
1
xm+1 (A) xm+2 (A) · · · xn (A)
1 
 xm+1 (A) xm+2 (A) · · · xn (A) 
 
, and
Let the matrix X2 (A) =   xA = 
 ···
···
··· ··· 
 · · ·
xm+1 (A) xm+2 (A) · · · xn (A)
1
 
1
1 
 
X2 (B) =   xB .
 · · ·
1




Let the extended matrix XA = X1 X2 (A) , and XB = X1 X2 (B) .
Let f A be the vector of which the i-th element is the function value of
the i-th row of matrix XA , i.e., f A = f (XA ), and f B is similarly defined,
f B = f (XB ).

Lemma 2. The two vectors f A and f B are linearly related if the function
f (x) is separable with respect to the first m variable combination (x1 , x2 , · · · , xm ).
Proof. Since the first m variable combination (x1 , x2 , · · · , xm ) are separable,
from definition 1, we have f (x) = ϕ1 (x1 , x2 , · · · , xm )⊗ϕ2 (xm+1 , xm+2 , · · · , xn ).
Accordingly, the vector f A = f (XA ) = ϕ1 (X1 ) ⊗ ϕ2 (X2A ) = ϕ1 (X1 ) ⊗ kA ,
where ⊗ is a component-wise binary operation, and kA = ϕ2 (xA ) is a scalar.

7

Similarly, the vector f B = ϕ1 (X1 ) ⊗ kB . As a result,

if ⊗ is times
 kA /kB · f B
fA =
kA − kB + f B if ⊗ is plus

kB − kA + f B if ⊗ is minus

which means the two vectors f A and f B are linearly related.

On the other hand, if the first m variables are fixed to two given points
C and D, and the rest of n − m variables are randomly sampled. A similar
proposition
 could
 be concluded as follows. Let

1
x1 (C) x2 (C) · · · xm (C)
1 
  x1 (C) x2 (C) · · · xm (C) 
 
,
X1 (C) =   x1 (C) x2 (C) · · · xm (C) = 
 ···
··· ···
··· 
 · · ·
x1 (C) x2 (C) · · · xm (C)
1
 
1
1 

 
X1 (D) =   x1 (D) x2 (D) · · · xm (D) ,
 · · ·

 1
(1)
(1)
(1)
xm+1 xm+2 · · · xn
 (2)
(2) 


 xm+1 x(2)
· · · xn 
m+2
X2 = 
 , the N × n matrix XC = X1 (C) X2 ,
··· ··· ··· 
 ···
(N )
(N )
(N )
xm+1 xm+2 · · · xn


and XD = X1 (D) X2 . Let f C be the vector of which the i-th element is
the function value of the i-th row of matrix XC , i.e., f C = f (XC ), and f D
is similarly defined, f D = f (XD ).
Lemma 3. The two vectors f C and f D are linearly related if the function
f (x) is separable with respect to the first m variable combination (x1 , x2 , · · · , xm ).
Proof. Since the first m variable combination (x1 , x2 , · · · , xm ) are separable,
from definition 1, we have f (x) = ϕ1 (x1 , x2 , · · · , xm )⊗ϕ2 (xm+1 , xm+2 , · · · , xn ).
Accordingly, the vector f C = f (XC ) = ϕ1 (X1C ) ⊗ ϕ2 (X2 ) = kC ⊗ ϕ2 (X2 ),
where ⊗ is a component-wise binary operation, and the scalar kC = ϕ1 (xC ).
Similarly, the vector f D = kD ⊗ ϕ2 (X2 ). As a result,

kC /kD · f D
if ⊗ is times
fC =
kC − kD + f D if ⊗ is plus or minus

which means the two vectors f C and f D are linearly related.
8

The above lemmas show that two function-value vectors must be linearly
related if the target function has the separability feature, while the separable
variables (or their complement variables) are fixed. These are necessary
conditions for the separability identification of target function. The sufficient
and necessary conditions are given as follows.
Theorem 1. The function f (x) is separable with respect to the first m variable combination (x1 , x2 , · · · , xm ) if and only if both of the flowing statements
are true.
(1) Any two function-value vectors with fixed (x1 , x2 , · · · , xm ) are linearly related;
(2) Any two function-value vectors with fixed (xm+1 , xm+2 , · · · , xn ) are linearly related.
Proof. From Lemma 2, and Lemma 3, we can conclude that the necessary
conditions of the theorem hold. The sufficient conditions can be proved by
contradiction. Suppose the separable form f (x) = ϕ1 (x1 , x2 , · · · , xm ) ⊗
ϕ2 (xm+1 , xm+2 · · · , xn ) can not be derived from the above two conditions.
Thus, there is at least one non-separable variable presented in both subfunctions, ϕ1 and ϕ2 . Without loss of generality, we assume xm to be this
non-separable variable. That is,
f (x) = ϕ1 (x1 , x2 , · · · , xm ) ⊗ ϕ2 (xm , xm+1 , xm+2 · · · , xn ) .

(7)

Similarly, the process of sampling for the first correlation test can be given
as
 (1)
(1)
(1) 
x1
x2
· · · xm
 x(2) x(2) · · · x(2) 
m 

2
,
(8)
X1 =  1.
.
.. 
..
 ..
. 
(N )
(N )
(N )
x1
x2
· · · xm
 (1)
(A)
(A)
(A) 
xm xm+1 xm+2 · · · xn
 x(2) x(A) x(A) · · · x(A) 
m 
 m
(A)
m+1
m+2
X̃2 =  .
,
(9)
.
.
.. 
.
.
.
 .
.
.
. 
(N )

xm

(A)

(A)

(A)

xm+1 xm+2 · · · xm

9

and

(B)
(B)
(B) 
xm+1 xm+2 · · · xn
(B)
(B)
(B)

xm+1 xm+2 · · · xm 


(B)
(10)
X̃2 = 
.
.
.
.
.
.
.

.
.
. 
(N )
(B)
(B)
(B)
xm xm+1 xm+2 · · · xm
h
i
h
i
(A)
(B)
′
′
Let the extended matrix XA = X1 X̃2
, and XB = X1 X̃2
.
Thus,


(A)
f ′A = f (X ′ A ) = ϕ1 (X1 ) ⊗ ϕ2 X̃2
= ϕ1 (X1 ) ⊗ α,
(11)



(1)

xm
(2)
xm
..
.

and

f ′B

′

= f (X B ) = ϕ1 (X1 ) ⊗ ϕ2



(B)
X̃2



= ϕ1 (X1 ) ⊗ β
(12)


(A)
and
are defined, where α and β are function-value vectors of ϕ2 X̃2


(B)
ϕ2 X̃2
respectively. As a result,
f ′A


′
 γ·f B
α − β + f ′B
=

β − α + f ′B

if ⊗ is times
if ⊗ is plus
if ⊗ is minus

(13)

where γ = (α1 /β1 , α2 /β2 , · · · , αN /βN ). From the lemmas, we know that, two
vectors f ′A and f ′ B are linearly related if they are in the relation of f ′A =
k1 ·f ′ B +k2 , where k1 and k2 are constant scalars, k1 6= 0. But, from the above
discussion, the components of all the three vectors, γ, α − β and β − α, are

(1)
(2)
(N )
not constant, due to the randomness of sample points xm , xm , · · · , xm .
This contradicts the supposition that the two vectors f ′A and f ′ B are linearly
related, and so Equation (7) cannot hold.
4.3. Notes on BiCT
The proposed technique is called bi-correlation test (BiCT) since two
complementary correlation tests are simultaneously carried out to determine
whether a variable or a variable-combination is separable.
The above process is illustrated with two sub-functions, and it could be
extended to determine the separability of a function with more sub-functions.
However, if the binary operators ⊗1 , ⊗2 , · · · , ⊗m in equation (6) are mutually
different with mixed times and plus or minus, the extension might be a little
10

difficult. This issue will be left for the future work. Hereafter, we assume
that the binary operators in equation (6) are the same, i.e, ⊗i = times or
⊗i = plus or minus for all i = 1, 2, · · · , m, for simplicity. In this case, the
extension process is very easy and omitted here.
To enhance the stability and efficiency of the algorithm, the distribution
of sample points should be as uniform as possible. Therefore, controlled sampling methods such as Latin hypercube sampling (Beachkofski and Grandhi,
2002) and orthogonal sampling (Steinberg and Lin, 2006) are preferred for
sample generation.
For the correlation test, any of correlation methods could be used. That
is, Pearson’s r method, Spearman’s rank order correlation, and Kendall’s τ
correlation are all effective for BiCT.
Take the function f (x) = 0.8 + 0.6 ∗ (x21 + cos(x1 )) ∗ sin(x2 + x3 ) ∗ (x2 −
x3 ), x ∈ [−3, 3]3 , as an example, the first sample set consists of 13 uniformly distributed points in [−3, 3], and the second sample set consists of
169 uniformly distributed points in [−3, 3]2 . The correlation tests could be
illustrated as in Fig. 2. As can be seen that the function-value vectors f A
and f B are linearly related (Fig. 2(b)), in which the variable x2 and x3 are
fixed when considering the first variable x1 (Fig. 2(a)). Similarly, to find
out the separability of variable combination (x2 , x3 ), the first variable x1 is
fixed (Fig. 2(c)). The corresponding function-value vectors f C and f D are
linearly related (Fig. 2(d))
The D&C method with BiCT technique is described with functions of
explicit expressions. While in practical applications, no explicit expression is
available. In this case, some modifications need to adapt for D&C method.
In fact, for data-driven modeling problems, a surrogate model of black-box
type could be established as the underlying target function (Forrester et al.,
2008) in advance. Then the rest steps are the same as above discussions.
5. Numerical results
5.1. Analysis on Computing time
The computing time (t) of a genetic programming with the proposed
divide and conquer (D&C) method consists three parts:
t = t1 + t2 + t3

(14)

where t1 is for the separability detection, t2 for sub-function determination,
and t3 for function recover. Note that both the separability detection and
11

6
A

1

f

B

4

f

0
−1
f

B

0

−2

−2

−3

−4

−4
−2

−1

0
x1

1

2

−5

3

1

2

3

4

5

6

fA

(b) f A ↔ f B

(a) (x2 , x3 ) are fixed
150
100
50
D

−6
−3

f

f

2

0
−50
−100
−150
−100

−50

0
C

f

(d) f C ↔ f D

(c) x1 are fixed

Figure 2: Demo of separability detection process of BiCT

12

50

100

Table 1: A mapping table for parse-matrix evolution

a·1
T
a·2 , a·3
expr
a·4
f → fk

-5
√
·
-5
λ2
-1
skip

-4
ln
-4
λ1
0
f1

-3 -2
cos /
-3 -2
f f2
1
f2

-1
-1
f1

0
skip
0
1.0

1
+
1
x1

2
*
2
x2

3
sin
···
···

4
exp
d
xd

5
(·)2

function recover processes are double-precision operations and thus cost much
less time than the sub-function determination process. That is, t ≈ t2 .
It is obviously that the CPU time for determining all sub-functions (t2 )
is much less than that of determining the target function directly (td ). Next,
a typical genetic programming, parse matrix evolution (PME), is taken as
the optimization driver (other GP algorithms should also work) to show the
performance of D&C.
Suppose that the dimension of the target function is d, the height of
the parse-matrix h, and the mapping Table as in Table 1, then the parsematrix entries a·1 ∈ −5, −4, ..., 4, 5, a·j ∈ −5, −4, · · · , d, (j = 1, 2), and
a·4 ∈ −1, 0, 1. Thus the parse-matrix (aij )h×4 have as many as (11 ∗ (6 + d) ∗
(6 + d) ∗ 3)h possible combinations. Thus the CPU time of determining the
target function directly satisfies
td ∼ (11 ∗ (6 + d) ∗ (6 + d) ∗ 3)h .

(15)

This means that the searching time of determining a target function will
increase exponentially with model complexity. Using D&C, only the subfunctions are needed to determinate, and each sub-function has less dimension d and less complexity h. Therefore, it will cost much less CPU time.
In fact, by D&C powered GP, the CPU time will increase only linearly with
dimensions, provided that the target function is completely separable.
Take equation (2) in Section 1 as an example. Without D&C, to search
directly, the control parameters of PME should be set as follows: d = 5,
h ≥ 9. From equation (15), the order of required CPU time td = O(2.58·1032).
Using D&C powered PME, the required CPU time will be much less.
In fact, after the separation detection, the function is divided into for subfunctions as follows.
√
√
qs = 1.83 × 10−4 v 3 · ρ · 1/ R · (1 − hw /hs )
13

√
For the sub-function v 3 , 1/ R, and 1 − hw /hs , the control parameters of
√
PME should be set as d = 1, h ≥ 2. For the sub-function ρ, d = 1, h ≥ 1.
As a result, t2 ≈ 4∗O(2.61·106) = O(107) by equation (15), which means the
D&C method could reduce the computational effort by orders of magnitude.
5.2. Program timing
Next, two illustrative examples are presented to show how much time the
D&C technique could save in practical applications.
Again, equation (1) and equation (2) are set as the target function, respectively. For equation (1), the sample set consists 100 observations uniformly distributed in [1,10] degree and [1000, 10000] (i.e., θ=1:10; Rex =
1000:1000:10000). The angle θ is fixed to 5 degree while detecting the subfunction f1 (Rex ), and the Renold number Rex is fixed to 5000 while detecting
the sub-function f2 (θ).
For equation (2), the sample set consists 30000 observations uniformly
distributed in a box in R5 (i.e., v = 500 : 100 : 1000; ρ = 0.0001 : 0.0001 :
0.001; R = 0.01 : 0.01 : 0.1; hw = 10000 : 10000 : 50000; hs = 100000 :
100000 : 1000000). The velocity of free-stream v, density of air ρ, radius of
nose R, Wall enthalpy hw , and total enthalpy hs are fixed to 800 m/s2 , 0.0005
kg/m3 , 0.05 m, 20000 J/kg, and 200000 J/kg, respectively, while detecting
the sub-functions.
In both tests, the program stops when the current model is believed good
can be regarded as the
enough: 1 − R2 < 1.0 · 10−10 , where R2 = 1 − SSE
SST
fraction of the total sum of squares that has ‘explained by’ the model. To
suppress the affect of randomness of PME, 10 runs are carried out for each
target function, and the averaged CPU time on a PC using a single CPU
core (Intel(R) Core (TM) i7-4790 CPU @3.60GHz) is recorded to show its
performance. The test results (see Table 2, and Table 3) show that D&C
technique can save CPU time remarkably. For equation (1), PME needs
about 12 minutes and 26 seconds to get an alternative of the target function
without D&C technique, while the D&C powered PME needs only about 11.2
seconds, which is much faster than the original algorithm. Similar conclusion
could also seen from the test results of equation (2) (see Table 3). Note that
the total time of D&C powered PME includes t1 , t2 and t3 (See Equ. 14),
and t1 + t3 ≈ 0.2 for Equ. (1), 0.3 for Equ. (2).

14

Table 2: Performance of PME on detecting Equ. (1) (with and without D&C)

Target Function CPU time
For D&C powered PME
f1 (Rex )
3s
f2 (θ)
8s∗
Total time
11.2s
PME without D&C
f (Rex , θ)
12m26s∗
Total time
746s

Result Expression
√
St = 0.1978/ Rex
St = 0.03215 ∗ θ − 0.01319 ∗ θ3
√
St = 2.274 ∗ θ ∗ cos(0.9116 ∗ θ)/ Rex

PME failed to get the exact result, but always result in an alternative
function with fitting error of zero in double precision (i.e., 1 − R2 = 0.0).
∗

Table 3: Performance of PME on detecting target Equ. (2)(with and without D&C)

Target Function CPU time
For D&C powered PME
f1 (v)
3s
f2 (ρ)
2s
f3 (R)
4s
f4 (hw , hs )
9s
Total time
18.3s
PME without D&C
f (v, ρ, R, hw, hs) 85m43s
Total time
5143s

Result Expression
qs
qs
qs
qs

= 1.647 · 10−5 ∗ v 3
√
= 3.77 · 105 ∗ pρ
= 6.49 · 103 ∗ 0.08442/R
= 9370 − 9370 ∗ hw /hs

qs = 0.000183 ∗ v 3 ∗

15

p

ρ/R ∗ (1 − hw /hs )

6. Conclusion
The divide and conquer (D&C) method for symbolic regression has been
presented. The main idea is to make use of the separability feature of the
underling target function to simplify the search process. In D&C, the target
function is divided into a number of sub-functions based on the information
of separability detection, and the sub-functions are then determined by any
of a genetic programming (GP) algorithms.
The most important and fundamental step in D&C is to identify the
separability feather of the concerned system. To fulfill this task, a special
algorithm, bi-correlation test (BiCT), is also provided for separability detection in this paper.
The study shows that D&C can accelerate the convergence speed of GP
by orders of magnitude without losing the generality, provided that the target
function has the feature of separability, which is usually the case in practical
engineering applications.
References
References
Anderson, J., 2006. Hypersonic and High-Temperature Gas Dynamics (2nd
ed.). American Institute of Aeronautics and Astronautics, Inc., Virginia.
Beachkofski, B., Grandhi, R., APRL 2002. Improved distributed hypercube
sampling. In: 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural
Dynamics, and Materials Conference, Structures, Structural Dynamics,
and Materials and Co-located Conferences. AIAA paper no. 2002-1274.
Denver, Colorado.
URL https://doi.org/10.2514/6.2002-1274
Berenguel, L., Casado, L. G., Garca, I., Hendrix, E. M. T., Messine, F.,
2013. On interval branch-and-bound for additively separable functions with
common variables. Journal of Global Optimization 56 (3), 1101–1121.
Chen, C., Luo, C., Jiang, Z., 2017. Elite bases regression: A real-time algorithm for symbolic regression.
URL https://arxiv.org/abs/1704.07313

16

d’Avezac, M., Botts, R., Mohlenkamp, M. J., Zunger, A., 2011. Learning
to predict physical properties using sums of separable functions. SIAM
Journal on Scientific Computing 33 (6), 3381–3401.
Forrester, A., Sobester, A., Keane, A., 2008. Engineering design via surrogate
modelling: a practical guide. John Wiley & Sons, Hoboken, New Jersey.
Koza, J. R., 2008. Genetic programming: on the programming of computers
by means of natural selection. MIT Press, Cambridge, MA.
Luo, C., Zhang, S.-L., 2012. Parse-matrix evolution for symbolic regression.
Engineering Applications of Artificial Intelligence 25 (6), 1182–1193.
McConaghy, T., 2011. FFX: Fast, Scalable, Deterministic Symbolic Regression Technology. Springer New York, New York, NY, pp. 235–260.
O’Neill, M., Ryan, C., Aug. 2001. Grammatical evolution. IEEE Trans. Evol.
Comp 5 (4), 349–358.
URL http://dx.doi.org/10.1109/4235.942529
Schmidt, M., Lipson, H., 2009. Distilling free-form natural laws from experimental data. Science 324 (59236), 81–85.
Steinberg, D. M., Lin, D. K. J., 2006. A construction method for orthogonal
latin hypercube designs. Biometrika 93 (2), 279–288.

17

