A Constructive Learning Algorithm for an HME
Kazumi Saito and Ryohei Nakano
NTT Communication Science Laboratories
2 Hikaridai, Seika-cho, Soraku-gun, Kyoto 619-02 Japan
Telephone No. +81 774 95 1822, Fax No. +81 774 95 1839e-mail {saito,nakano}@cslab.kecl.ntt.jp
ABSTRACTA Hierarchical Mixtures of Experts (HME) model has been applied to several classes ofproblems, and its usefulness has been shown. However, deﬁning an adequate structure in advanceis required and the resulting performance depends on the structure. To overcome this problem,
a constructive learning algorithm for an HME is proposed; it includes an initialization method,a training method and an extension method. In our experiments, which used parity problems
and a function approximation problem, the proposed algorithm worked much better than theconventional method.

1. Introduction
A Hierarchical Mixtures of Experts (HME) model has been applied to several classes of problems, and its
usefulness has been shown [5, 6, 10]. However, deﬁning an adequate structure in advance is required andthe resulting performance depends on the structure.
To overcome this problem, we can consider two approaches: using pruning algorithms [8] or constructive
algorithms [7]. Pruning algorithms will play an important role in improving generalization performance of
networks, particularly for problems that involve many irrelevant input variables. However, a large amountof computation time must be spent to train larger than desirable initial networks; moreover, we generally
do not know what size is suitable for an initial network. Compared to pruning algorithms, we believe that
constructive algorithms are more suitable for the present problem. Although several algorithms have been
proposed for feed-forward networks or classiﬁcation trees (e.g.[3, 1]), little work has been done in an HME
context.
This paper proposes a constructive learning algorithm for HME. Section 2 explains the proposed algorithm, which includes an initialization method, a training method and an extension method. In Section 3 theproposed method is evaluated based on experiments that used parity problems and a function approximation
problem.

2. Constructive Algorithm
2.1. Basic deﬁnition
Let {(x1 , y1 ), · · · , (xm , ym )} be a set of examples, where xt denotes an n-dimensional input vector and yt atarget value corresponding to xt . An HME consists of several Expert Networks (EN s) and Gating Networks(GN s) expressed as a tree where leaves are EN s and the other nodes are GN s. In this paper, we only
consider binary HMEs, because an arbitrary tree can be reduced to a functionally equivalent binary tree.
Hereafter, the weight vector
by wi = (wi0 , · · · , win )T , and its output value is deﬁned
i is expressed 1 of ENT2
T
by ui (xt , yt ; wi ) = exp − 2 (yt − wi xt ) , where a means the transposed vector of a. The weight vector of
GNi is expressed by vi = (vi0 , · · · , vin )T , and since the soft-max function reduces to the sigmoidal function−1
in a binary HME, the output value of GNi can be deﬁned by gi (xt ; vi ) = 1 + exp(−viT xt ). Note thatwi0 and vi0 mean bias terms and xt0 is always set to 1. The simplest HME consisting of EN1 , EN2 , andGN1 is expressed as (GN1 , EN1 , EN2 ), and its output value is deﬁned by g1 u1 + (1 − g1 )u2 . Then, an
arbitrary HME can be expressed by a corresponding list structure such as (GN1 , EN1 , (GN2 , EN2 , EN3 )),and its output value can be recursively deﬁned by g1 u1 + (1 − g1 )(g2 u2 + (1 − g2 )u3 ).
2.2. Relation matrix and objective function
In order to describe the structure of variously extended HMEs, we introduce a matrix R = (rij ); here, rijrepresents a relation between ENi and GNj . An element only takes one of three values {1, −1, 0}: rij = 1
indicates that the product term gj ui appears in the HME; rij = −1 indicates that the term (1 − gj )ui
appears; rij = 0 indicates that ENi and GNj have nothing to do with each other. For example, the relation
matrix for (GN1 , EN1 , EN2 ) is R(2) = (1, −1)T , and the relation matrix for (GN1 , EN1 , (GN2 , EN2 ,EN3 )) is R(3) = ((1, −1, −1)T , (0, 1, −1)T ).
Hereafter, for the HME(c) whose number of EN s is c, a vector consisting of all weight vectors is simply
T
expressed as Φ(c) = (w1T , · · · , wcT , v1T , · · · , vc−1
)T , where N (c) = (2c − 1)(n + 1) is the dimension of Φ(c) .
(c)Then, by referring to the relation matrix R , the product term with respect to ENi can be expressed as
gj (xt ; vj ) ×(1 − gj (xt ; vj )).
(1)
hi (xt , yt ; Φ(c) ) = ui (xt , yt ; wi ) ×(c){j|rij =1}(c){j|rij =−1}

Now, the objective function of the HME(c) can be deﬁned by the logarithmic likelihood
 c

m


(c)
(c)
L(Φ ) =log
hi (xt , yt ; Φ ) .
t=1

(2)

i=1

Thus, the maximum likelihood estimation problem is ﬁnding the Φ(c) that maximizes L(Φ(c) ).
2.3. Main algorithm
The basic idea of the proposed method is that after training the HME(c) , the ENb with the largest error
is selected among the ENi , i = 1, · · · , c; then, by replacing ENb with (GNc , ENb , ENc+1 ), an extended
HME(c+1) is constructed. This is based on the “divide-and-conquer” approach employed by existing methods
such as CART [1]. The main algorithm is described as follows:
step 1: Initialize Φ(2) (w1 , w2 and v1 ), and set R(2) = (1, −1)T , c = 2;
step 2: Train all weights included in the HME(c) ;step 3: Terminate the iteration if a stopping criterion is satisﬁed;
step 4: Select an ENb to be extended, calculate R(c+1) , and initialize wc+1 and vc ;step 5: Set c = c + 1, and return to Step 2;
2.4. Initialization method (Step 1)
In our preliminary experiments using parity problems, when the initial values for Φ(2) were set to random(2)values near 0, most trials convergedmto Φ̂ (wT1 = 2w2 = ŵ, v1 = 0), where ŵ is the solution of a regression
problem: ﬁnd ŵ that minimizes t=1 (yt − w xt ) . Moreover, even though a small vector ∆Φ was addedto Φ̂(2) , it usually held that L(Φ̂ + ∆Φ) < L(Φ̂), and in most trials Φ(2) came back to the same point Φ̂(2) .
(2)Here, we analyze how this occurs by investigating the Hessian matrix ∇2 L(Φ0 ). In general, we canmm
assume that t=1 xt xTt is positive deﬁnite. Also note that t=1 (yt − ŵT xt )xt = 0. When Φ̂(2) is set as
w1 = w2 = ŵ and v1 = 0, we can see that ∇L(Φ̂(2) ) = 0, because∇wi L(Φ̂(2) ) =mmhi (xt , yt ; Φ̂(2) )(yt − wiT xt )1
xt =(yt − ŵT xt )xt = 0,
2
(2) )
2
h
(x
,
y
;
Φ̂
k
t
t
k=1
t=1
t=1

∇v1 L(Φ̂(2) ) =m
h1 (xt , yt ; Φ̂(2) )(1 − g1 (xt ; v1 )) − h2 (xt , yt ; Φ̂(2) )g1 (xt ; v1 )xt = 0,2(2) )k=1 hk (xt , yt ; Φ̂
t=1

where h1 (xt , yt ; Φ̂(2) ) = h2 (xt , yt ; Φ̂(2) ) and g1 (xt ; v1 ) = 0.5 for Φ̂(2) . On the other hand, by calculating
second-order derivatives, we can see that ∇2 L(Φ(2) ) usually becomes semi-negative deﬁnite, when the target
values are normalized in the range of [0, 1]. Thus, Φ(2) is a weak local maximum.
Fortunately, we can easily escape from this weak local maximum by using ∆Φ such that ∆w1 = ∆w2 = 0
v1 ) around0.5, eachand ∆v1 is set to a random vector. Here, in order to randomly scatter the value of g1 (xt ; nmvalue of {v11 , · · · , v1n } is set to a random value in the range of [−1, 1], and v10 is set to − i=1 v1i ( t=1 xti ).
That is, v1 becomes a random hyper-plane that includes the gravity of input vectors; training examples
will be divided into two classes eachm of which includes much the same number of examples. Then, L(Φ̂) =L(Φ̂ + ∆Φ), ∇w1 L(Φ̂ + ∆Φ) = t=1 g1 (xt ; v1 )(y − ŵT xt )xt = 0, and the value of the objective function
can be steadily increased in the next iteration.
2.5. Training method (Step 2)
2.5.1. Second-order learning algorithm
The conventional learning algorithm for an HME [6, 10] is based on the Expectation Maximization (EM)algorithm [2]. However, since this algorithm requires Iterative Reweighted Least Squares (IRLS) in the
M-step, numerical instabilities are likely [10]. Since our preliminary experiments also suﬀered from thisproblem, we employed a second-order learning algorithm [9] based on a quasi-Newton method [4].
In this algorithm, the following are repeated until convergence: after calculating the gradient vector, thesearch direction (∆Φ) is calculated on the basis of the partial BFGS (Broydon-Fletcher-Goldfarb-Shanno)
update; the optimal step-length λ that maximizes L(Φ + λ∆Φ) is calculated as the maximal point of asecond-order approximation. Below, we show that the optimal step-length can be eﬃciently calculated for
the HME(c) networks deﬁned in Section 2.2.
Since λ is the only variable in L(·), we can express L(Φ + λ∆Φ) simply as ζ(λ); then the maximal pointof a second-order approximation is given when λ = −ζ (0)/ζ  (0). By diﬀerentiating ζ(λ) and substituting
0 for λ, we obtain

c
c
m c
m
2


hi (xt , yt ; Φ)hi (xt , yt ; Φ)hi (xt , yt ; Φ)


i=1
i=1
c
+ i=1
ζ (0) =,
ζ (0) =− c.
c
i=1 hi (xt , yt ; Φ)i=1 hi (xt , yt ; Φ)i=1 hi (xt , yt ; Φ)
t=1
t=1
d
hi (xt , yt ; Φ + λ∆Φ)|λ=0 , we obtain
Now that the derivative of hi (xt , yt ; Φ) is deﬁned as dλhi (xt , yt ; Φ) = hi (xt , yt ; Φ)αi ,hi (xt , yt ; Φ) = hi (xt , yt ; Φ)αi + hi (xt , yt ; Φ)βi .

where
αi
βi

= (y − wiT xt )∆wiT x +=∆vjT x −{j|rij =1}



gj (xt ; vj )∆vjT x{j|rij ∈{1,−1}}−(∆wiT x)2 −gi (xt ; vj )(1 − gj (xt ; vj ))(∆vjT x)2 .
{j|rij ∈{1,−1}}Here, recall that ∆wi and ∆vi mean the search directions with respect to wi and vi , respectively.
2.5.2. Weight reduction
As construction of the HME proceeds, the magnitude of weight vectors in gating networks generally becomesvery large; since these output values approach 0 or 1, their derivatives almost become 0. Since these weights
can hardly be modiﬁed, the magnitude of such weights should be reduced; this reduction is equivalent
to decreasing the gain (or slope) of the sigmoid nonlinearity. In our preliminary experiments, when themagnitude of all gating weights was simultaneously reduced, a small reduction did not have much eﬀect,
but a large reduction was likely to destroy what the HME had learned so far. Since it is diﬃcult to know anadequate reduction in advance, we employed a method where the weight reduction was performed separately.
(new)
(old)
After training the HME(c) , the HME(c) is trained again by vk= γvk . This process is iterated fromk = 1 to k = c − 1. In our experiments, γ was set to 0.1.

Table: 1: Constructive HME for parity problems

parity bits
converged trials (out of 100)
average number of EN s
standard deviation number of EN s
minimal number of EN s

4
100
3.07
0.26
3

5
100
3.85
0.43
3

6
100
4.01
0.10
4

7
100
4.20
0.47
4

8
99
5.23
0.75
5

2.6. Extension method (Step 4)
In order to extend the HME(c) , the ENi that maximizes the following formula is selected:
Err(ENi ) =m(yt − wiT xt )2 ×gj (xt ; vj ) ×(c){j|rij =1}

t=1

(1 − gj (xt ; vj )).

(3)

(c)
{j|rij =−1}

In this formula, the ﬁrst squared value represents the error of each example, and the remaining productterm represents the probability that the example is assigned to ENi ; thus, the total summation can be
regarded as the expected error of ENi . The ENb with the largest error is selected.
The values of R(c+1) are determined as follows. First, since ENc+1 is located under the former locationof ENb , for 1 ≤ j ≤ c − 1, each element rc+1,j is set equal to rbj . Then, since GNc takes eﬀect only on ENb
and ENc+1 , rbc = 1, rc+1,c = −1, and for i = b, each element ric is set equal to 0.
The values of wc+1 and vc are initialized using the weight initialization method described above; namely,
wc+1 = wb , each value of {vc1 , · · · , vcn } is set to a random value in the range of [−1, 1], and vc0 is set to


n
m






vc0 = −vci 
xti ×
gj (xt ; vj ) ×(1 − gj (xt ; vj )) .
i=1

t=1

(c)

{j|rbj =1}(c){j|rbj =−1}

Here, vc becomes a random hyper-plane that includes the gravity of weighted input vectors with respect toENb .

3. Evaluation by experiments
3.1. Parity problems
By using 4- to 8-bit parity problems, the proposed method was evaluated. In this experiment, all possible
input patterns were used as training examples, the target values were set to 0 or 1, and the maximum(c)(c)< 10−8 ,number of EN s was set to 8. The iteration of Step 2 was terminated when ∇L(Φc ) /N
while the iteration of the main loop was terminated when each example satisﬁed i=1 Err(ENi ) < 0.1.
Table 1 shows the learning results; trials were performed 100 times for each number of bits. The minimum
number of EN s to solve an n-bit parity problem is given by n/2 + 1. Table 1 shows that most trials of the
proposed method converged by using the minimum number of EN s. The conventional method [10] required
many more EN s to solve these parity problems. For example, in the case of the 8-bit parity problem, it
required 64 EN s in order to achieve a success rate of 67%.
3.2. Function approximation problem
Consider a piecewise linear function that consists of several lines whose slope is −4 or 4; these lines arealternately connected to each other from (x, y) = (0, 2) to (4, 2), x is an input value and y is a target output
value. In this experiment, the value of x was randomly generated in the range of [0, 4], and the corresponding
value of y was calculated from x, where each value of y was corrupted by adding noise generated according

2Y

GN1 GN2 GN7 GN6 GN3

GN4 GN5

2

YGN5 GN2 GN4 GN1 GN7GN3 GN6
EN6

0

X
0 EN EN 1 EN
2
3 EN EN 431
7 EN8 EN5 EN4
6 2
EN1

0
2
EN1 EN3 EN41EN7 EN EN
8
2

GN1
GN2

GN5
EN2

EN6

EN5

EN4

X

GN4
GN7

(a) Result of constructive algorithm

GN3

GN2

GN6
EN7

4

GN1
EN3

GN3
GN4

3
EN5

EN8

EN1

EN2

GN5
EN3

EN4

GN6
EN5

EN6

GN7
EN7

EN8

(b) Results of 3-level balanced binary tree

Fig. 1: HMEs for function approximation

to a normal distribution with a mean of 0 and a variance of 0.1. The total number of training examples was
set to 100, and the same stopping criteria as in the previous experiment was employed. In all 10 trials, ourconstructive algorithm found desirable results. An example of the learning results is shown in Figure 1(a).
Its structure is expressed by(GN1 , EN1 , (GN2 , (GN3 , (GN4 , (GN5 , EN2 , EN6 ), EN5 ), (GN6 , EN4 , (GN7 , EN7 , EN8 )), EN3 ))).
On the other hand, when a 3-level balanced binary tree expressed by(GN1 , (GN2 , (GN4 , EN1 , EN2 ), (GN5 , EN3 , EN4 )), (GN3 , (GN6 , EN5 , EN6 ), (GN7 , EN7 , EN8 ))),was given in advance, we were not able to obtain a reasonable result during 10 trials. Figure 1(b) shows
an example of the learning results. In order to obtain a desirable result using the 3-level balanced binary
tree, the decision boundary of GN1 must be x = 2. Actually, in the case of Figure 1(b), the decisionboundary of GN1 was x ≈ 1.5; thus, redundant EN s appeared when x < 1.5. Conversely, more EN s
were required when x > 1.5. When the structure of the HME is ﬁxed in advance, learning will be morediﬃcult because the decision boundaries of several GN s are predetermined. Figure 2 shows the learning
process of our constructive algorithm; although the decision boundary of GN1 was x ≈ 0.5, our constructive
algorithm was able to obtain a desirable result because new GN s with adequate boundaries were successively
generated.

4. Conclusion
We have proposed a constructive learning algorithm for HMEs that includes an initialization method, a
training method and an extension method. In experiments using parity problems and a function approximation problem, the proposed algorithm worked better than the conventional algorithm in respect that
our algorithm obtained the desirable results with the minimum size of HME. In the future, we plan to do
further comparisons using a wider variety of problems.

[#EN=2]2[#EN=5]
2

GN1

GN1

GN2

GN3

GN4

EN2
0

0
0

1
EN1
GN1 GN2

2

0 EN EN 1
3
1[#EN=4]
GN1 GN2
2

2[#EN=3]
2

3

4

0
1
EN4 EN3 EN1[#EN=6]
GN1 GN222
3
EN5 EN2
GN3 GN4

4
GN5

EN2
0

0
0
1
EN4 EN3 EN1 EN2

3

4

3

4

GN3

0
0
1
EN4 EN3 EN1[#EN=7]
GN1 GN2
2

2
EN5

3

4
EN6 EN2
GN6 GN3 GN4 GN5

2

3

0
2

0

1
EN3 EN1 EN7

EN5 EN4

4
EN6 EN2

Fig. 2: Learning process of constructive algorithm

References[1] L. Breiman, J. Friedman, R. Olshen, and C. Stone, Classification and Regression Trees. Montrey, CA:
Wadsworth International Group, 1984.
[2] A. Dempster, N. Laird, and D. Rubin, “Maximum likelihood from incomplete data via EM algorithm,”
J. Royal Statist. Soc. Ser. B (methodology), vol. 39, pp. 1–38, 1977.
[3] S. Fahlman and C. Libiere, “The cascade-correlation learning architecture,” in Advances in Neural
Information Processing Systems 2 (D. Touretzky, ed.), (Los Altos, CA), pp. 524–532, Morgan Kaufmann,
1990.
[4] P. Gill, W. Murray, and M. Wright, Practical optimization. London: Academic Press, 1981.
[5] R. Jacobs, M. Jordan, S. Nowlan, and G. Hinton, “Adaptive mixtures of local experts,” Neural Computation, vol. 3, no. 1, pp. 79–87, 1991.
[6] M. Jordan and R. Jacobs, “Hierachical mixtures of experts and EM algorithm,” Neural Computation,
vol. 6, no. 2, pp. 181–214, 1994.
[7] T. Kwok and D. Teung, “Constructive feedforward neural networks for regression problems: a survey,”
tech. rep., HKUST-CS95-43, 1995.
[8] R. Reed, “Pruning algorithms - a survey,” IEEE Trans. Neural Networks, vol. 4, no. 5, pp. 740–747,
1993.
[9] K. Saito and R. Nakano, “A connectionist approach to numeric law discovery,” in Machine Intelligence,
vol. 15, 1995 (to appear).
[10] S. Waterhouse and A. Robinson, “Classiﬁcation using hierarchical mixtures of experts,” in Proc. of
NNSP, pp. 177–186, 1994.