Universal Approximation Using Shuffled Linear Models

arXiv:1308.6498v2 [math.DS] 19 Nov 2013

Laurens Bliek
Faculty of Electrical Engineering, Mathematics, and Computer Science, Delft University
of Technology, Mekelweg 4, 2628 CD Delft, The Netherlands

Abstract
This paper proposes a specific type of Local Linear Model, the Shuffled Linear
Model (SLM), that can be used as a universal approximator. Local operating points are chosen randomly and linear models are used to approximate a
function or system around these points. The model can also be interpreted
as an extension to Extreme Learning Machines with Radial Basis Function
nodes, or as a specific way of using Takagi-Sugeno fuzzy models. Using the
available theory of Extreme Learning Machines, universal approximation of
the SLM and an upper bound on the number of models are proved mathematically, and an efficient algorithm is proposed.
Keywords: neural networks, system identification, nonlinear systems
1. Introduction
The approximation of nonlinear functions or systems is a problem that
appears in many disciplines [1] and that can be tackled by various approaches.
One of these approaches is to decompose the problem into smaller problems
and use a combination of simple solutions that only work for a small part
of the problem. The combination of these local models should then approximate the global function or system. The Local Linear Model (LLM) is a
key example of this strategy, where a dynamical system is partitioned into
multiple operating regimes, and for each regime a simple linear model is used
to model the system.
If the goal is to approximate a function, not to model a system, the
use of artificial neural networks has been a popular approach. One of the
most popular algorithms for approximating a function with neural networks
is the backpropagation algorithm. The Extreme Learning Machine (ELM)
algorithm as proposed by Huang et al. [2], is an alternative algorithm which
Preprint submitted to Neurocomputing

February 6, 2018

outperforms the backpropagation algorithm in many aspects. Instead of
tuning the parameters of the hidden layer and output layer, the hidden layer
parameters are initialised randomly and kept fixed while the output weights
are trained by a linear least-squares method. Since very efficient algorithms
are available for solving linear least-squares problems, the ELM algorithm is
an efficient method for training a feedforward neural network.
Although ELMs can use nodes with many functions that do not at all
resemble the activation functions observed in biological neurons [3], which
makes ELM a very general approach, to the author’s knowledge they have
not yet been extended in such a way that they include LLMs. The activation
functions can be very general, but the hidden network nodes typically take
one of two forms: either additive nodes, or Radial Basis Function (RBF)
nodes are used, while the output is a linear combination of these. Unfortunately, with this architecture, this type of network can not be made equivalent
to a LLM as defined in this paper. To prove that the same ELM algorithm
does work for LLMs, an extension has to be made.
The ELM algorithm is extended in this paper to include the use of LLMs,
and the approximation capabilities are analysed and proved mathematically.
Since ELMs use randomly initialised hidden layer parameters, the linear models used in the LLM in this paper will be accurate around a randomly chosen
local point. This can be seen as shuffling a deck of cards and handing them
to the players in a cardgame, where the cards consist of parameters that
indicate the location and impact factors of the local points, while the players
are the linear models. The proposed model in this paper is therefore called
Shuffled Linear Model (SLM). Note that not the linear models themselves,
but the localities where they are most valid are chosen randomly.
Extending the ELM algorithm to include LLMs gives several advantages,
like a better physical interpretation and the potential to identify nonlinear
systems, making the algorithm more fit for system identification. Besides
this, the extension could be used in the area of fuzzy modelling. Fuzzy
models are a popular approach to model nonlinear systems, and can be used
to combine both rigid mathematical concepts and vaguer linguistic concepts.
One of the most popular fuzzy models is the Takagi-Sugeno model [4] (TSM).
TSMs use linear consequences in their fuzzy rules, making them equivalent to
LLMs under some circumstances. This is also shortly analysed in this paper,
but the main point of this paper is the extension of the ELM algorithm.
The paper has the following structure: in Section 2, a summary of the
ELM approach will be given. Section 3 contains the proposed model with
2

the definition of the SLM. Section 4 provides proof of the approximation
capabilities of the SLM and contains the proposed algorithm. Section 5
contains several ways to interpret the SLM. Conclusions are given in Section
7.
2. ELM summary
The ELM architecture is similar to a feedforward neural network with one
hidden layer. The output neurons have a linear activation function, while the
hidden neurons have a nonlinear activation function. The main difference
with traditional neural networks is in the tuning of the parameters: the
hidden layer parameters are initialised randomly and remain fixed, while only
the output weights are tuned. This leads to a linear least-squares problem.
Although the original ELM used neural alike hidden nodes, ELMs have
been generalised to include many activation functions that are not neuron
alike [5]. The most common example is the Radial Basis Function (RBF)
2
g(a, b, x) = e−b||x−a|| , with parameters a ∈ Rn , b ∈ R>0 , and input x ∈ Rn .
In this paper, only ELMs with RBF nodes are considered. For a network
with h hidden nodes, the output y ∈ Rm of the ELM is:
y=

h
X

βi g(ai , bi , x) =

i=1

h
X

βi gi (x),

(1)

i=1

with βi ∈ Rm the output weight vector connecting the i-th hidden node with
the m output nodes. For N inputs xj , this can be written in matrix notation
as
Y = GB,
(2)
where

and




g1 (x1 ) . . . gh (x1 )


..
..
G=

.
.
g1 (xN ) . . . gh (xN ) N ×h



y1T


Y =  ... 
,
T
yN N ×m

(3)




β1T


B =  ... 
.
T
βh h×m

3

(4)

If a similar N × m-matrix Y ∗ consisting of target output values y ∗ is
constructed, the goal is to minimise the sum of square errors:
E=

N
X

||yj − yj∗ ||2 = ||Y − Y ∗ ||2F = ||GB − Y ∗ ||2F ,

(5)

j=1

where ||·||F denotes the Frobenius norm. The output weights B that minimise
E, can be found by using the Moore-Penrose pseudo-inverse [6], here denoted
as †. This also makes sure that the norm of the output weights are minimised.
Then the optimal solution is:
B̂ = G† Y ∗ .

(6)

3. Proposed model
The model proposed in this paper is an extension of the ELM model. To
define the model, two other definitions are given first:
Definition 1. A Linear Model (LM) is a function of the form LMi (x) =
αi x + βi , with αi an m × n-matrix and βi ∈ Rm .
Definition 2. A Local Linear Model (LLM) is a weighted sum
Ph of LMs, where
the model weights pi can depend on the input: LLM (x) = i=1 pi (x)LMi (x)
P
= hi=1 pi (x)(αi x + βi ).
Now the proposed model will be defined:
Definition 3. A Shuffled Linear Model (SLM) is a LLM with RBF model
2
weights pi (x) = gi (x) = e−bi ||x−ai || , where the parameters ai ∈ Rn and
bi ∈ R>0 are drawn from aP
continuous probability distribution (independent
2
of the input): SLM (x) = hi=1 e−bi ||x−ai || (αi x + βi ).
Note that the SLM is equivalent to the ELM if all αi are zero. It is
therefore expected that, for nonzero αi , the approximation capabilities of
the SLM are better than those of the ELM if the number of local models is
equal to the number of hidden nodes h. This is what will be proved in this
paper.

4

4. Approximation capabilities of SLM and proposed algorithm
In this section it will be proved that a SLM, as an extension of an ELM,
can act as a universal approximator, and an upper bound on the number
of models will be given. First, it will be proved that the SLM can act as a
universal approximator. After that, several lemmas will be proved that lead
to the main theorem of this paper, where Lemma 2 is used in the proof of
Lemma 3, and Lemmas 1 and 3 are used in the proof of the main theorem.
This section concludes with a learning algorithm for the SLM.
4.1. Universal approximation with a SLM
The result of the first theorem follows directly from the approximation
capabilities of the ELM:
Theorem 1. SLMs are universal approximators. That is, given a continuous
target function f : Rm → Rn and a SLM with free parameters αi and βi ,
i = 1 . . . h, then with probability
one the parameters αi and βi can be chosen
R
in such a way that limh→∞ Rn ||SLM (x) − f (x)||2 dx = 0.
Proof. Choosing αi = 0, i = 1 . . . h, the SLM is equivalent to an ELM. For
an ELM, the parameters βi can be determined using linear regression to
guarantee universal approximation [5, Thm. II.1].
The interesting part is, of course, whether the extension from ELM to
SLM gives any improvements. It turns out that a SLM can also be trained
using linear regression with all free parameters, just like an ELM. But for the
same approximation accuracy, it needs less local models than the number of
hidden nodes in the ELM. This is simply due to the fact that a SLM has more
free parameters. The remainder of this section is devoted to the proof of these
results, after which an algorithm for training the SLM will be presented.
4.2. Necessary lemmas
Before presenting the next theorem, some lemmas will be proved. Most
proofs are similar to those of the ELM theorems [2, 5].
Lemma 1. Given a SLM with inputs xj and desired outputs yj∗ , j = 1 . . . N ,
P
∗ 2
the error E = N
j=1 ||SLM (xj ) − yj || can be minimised using linear regression.

5

Proof. Let γi = [αi βi ] be the m × (n + 1)-matrix of the free parameters
for model i, and let Γ = [γ1 . . . γh ]T be the (n + 1)h × m-matrix of all free
xj
parameters. Let also zj =
, then the output of the SLM can be written
1
as:
h
X
SLM (x) =
gi (x)γi z.
(7)
i=1

Now the error E can be rewritten as follows:
E =

N
X

||SLM (xj ) − yj∗ ||2

j=1
N
h
X
X
||
gi (xj )γi zj − yj∗ ||2
=
j=1

i=1




g1 (xj )zj


..
∗ 2
=
||ΓT 
 − yj ||
.
j=1
gh (xj )zj
N
X

(8)

= ||ΓT K T − Y ∗ T ||2F ,
= ||KΓ − Y ∗ ||2F ,
with

and




y1∗ T


,
Y ∗ =  ... 
∗ T
yN
N ×m

(9)


g1 (x1 )z1T . . . gh (x1 )z1T


..
..
.
K=

.
.
T
T
g1 (xN )zN . . . gh (xN )zN N ×(n+1)h

(10)



This sum of squared errors can be minimised using linear regression, for
example by using the minimum norm least-square solution
Γ̂ = K † Y ∗ .
This completes the proof.

6

(11)

Lemma 2. Given N distinct samples (xj , yj∗ ) ∈ Rn ×Rm , if a ∈ Rn is chosen
randomly from a continuous probability distribution, then with probability one
we have ||xj − a|| =
6 ||xj 0 − a||, for j 6= j 0 , j, j 0 = 1 . . . N .
Proof. The proof is an adaptation to one of the proofs of Huang and Babri [7].
Consider the set V (xj , xj 0 ) = {a ∈ Rn : ||xj − a|| = ||xj 0 − a||}. This set is a
hyperplane in Rn . Since this is an (n−1)-dimensional
surface in Rn , and there
S
0
are only finite j, j , the union U = j6=j 0 V (xj , xj 0 ) is a finite union of (n − 1)dimensional surfaces for j 6= j 0 , j, jR0 = 1 . . . N . So, for any probability density
function f , we have P (a ∈ U ) = U f (a)da = 0. Therefore, the probability
that ||xj − a|| = ||xj 0 − a|| is 0 for randomly chosen a, and the result of the
lemma follows.
Lemma 3. For a SLM , matrix K from equation (10) has full rank with
probability one if the matrix Z = [z1 . . . zN ]T has full rank.
Proof. For i = 1 . . . h, and k = 1 . . . n + 1, the (i × k)-th column of K is given
2
by [gi (x1 )z1k . . . gi (xN )zNk ]T , where gi (x) = g(ai , bi , x) = e−bi ||x−ai || . Let
I = (u, v) be any interval from R, with u < v. Let c : I → RN be the curve
defined as
c(bi ) = [gi (x1 ) . . . gi (xN )]T ,
(12)
seen as a function of bi , for bi ∈ I. Using a similar proof as that of Tamura
and Tateishi [8], it can be proved by contradiction that c(bi ) does not belong
to a subspace with dimension less than N .
Suppose it does, then there exists a vector w 6= 0 orthogonal to this
subspace:
T

w (c(bi ) − c(u)) =

N
X

wj g(ai , bi , xj ) − wT c(u) = 0.

(13)

j=1

Without loss of generality, assume wN 6= 0, then
N −1

X wj
1 T
w c(u) −
g(ai , bi , xj ).
g(ai , bi , xN ) =
wN
wN
j=1

(14)

On both sides of this equation are functions of bi . Taking the s-th partial
derivative to bi on both sides gives the following equation:
N −1

X wj ∂ s
∂s
g(a
,
b
,
x
)
=
−
g(ai , bi , xj ),
i
i
N
∂bsi
wN ∂bsi
j=1
7

s = 1, 2, . . . .

(15)

Now, g is infinitely differentiable w.r.t. bi , with
||2

∂s
g(ai , bi , xj )
∂bsi

= (−||xj −

ai ||2 )s e−bi ||xj −ai . Using this fact and the result of Lemma 2, it follows that
with probability one, Equation (15) actually contains an infinite number of
equations that are linear in parameters wj , for s = 1, 2, . . .. The number of
free parameters wj , however, is N . This gives a contradiction. Hence, c(bi )
does not belong to a subspace with dimension less than N .
Since the above contradiction holds for bi from any interval I ⊆ R, it
is possible to choose h values b1 , b2 , . . . , bh from any continuous probability
distribution over R, such that the following matrix G has full rank:
G = [c(b1 ) c(b2 ) . . . c(bh )].

(16)

This is the same N × h matrix as the hidden layer output matrix from
equation 3. The matrix Z = [z1 . . . zN ]T is an N × (n + 1)-matrix. Let
M = G ⊗ Z be the Kronecker product of these matrices. Then M is an
N 2 × (n + 1)h-matrix. Since both G and Z have full rank, and for matrices
the general rule rank(A ⊗ B) = rank(A)rank(B) holds, the matrix M has
rank r = min(h, N ) min(n + 1, N ).
Now, matrix K appears actually inside matrix M . For matrix M , we
have:


g1 (x1 )z1T . . . gh (x1 )z1T
..
..




.
.

T
T 
 g1 (x1 )zN . . . gh (x1 )zN 


..
..
.
(17)
M =G⊗Z =
.
.


 g (x )z T . . . g (x )z T 
 1 N 1
h N 1 


..
..


.
.
T
T
g1 (xN )zN
. . . gh (xN )zN
The N × (n + 1)h-matrix K from equation (10) follows from this matrix
by removing rows until only the rows with indices (j − 1)(N + 1) + 1 are
left, for j = 1, . . . , N . It follows that the rank of matrix K is equal to
min(N, r) = min(N, (n + 1)h), so matrix K has full rank.
4.3. Improvement of the SLM as an extension of an ELM
Using the Lemmas that appeared in this section, the main theorem of
this paper, which is similar to that of Huang [2, Thm. 2.2] can be presented.

8

From the proofs of the Lemmas, the following definitions are needed:


xj
zj =
,
1
Z = [z1 . . . zN ]T .

(18)
(19)

Theorem 2. Given any  > 0, and given N distinct samples (xj , yj∗ ) ∈
Rn × Rm , if matrix Z has full rank, then there exists h ≤ N/(n + 1) such
that aPSLM with h local models can be trained using linear regression to get
∗ 2
<  with probability one.
E= N
j=1 ||SLM (xj ) − yj ||
Proof. The error E can actually become zero by choosing h = (n + 1)N .
Since Z has full rank, from Lemma 3 it follows that matrix K also has full
rank, with probability one. Then K is also invertible, since the size of K is
N ×(n+1)h = N ×N . The solution showed at the end of the proof of Lemma
1 is now actually equal to Γ̂ = K −1 Y ∗ , and the error E = ||KΓ − Y ∗ ||2F is
zero.
Since zero error can be achieved by choosing h = (n + 1)N , there always
exists a h ≤ (n + 1)N to let the error be as small as desired.
Note that the condition that Z has full rank can be easily satisfied, for
example by adding noise to the training samples. Theorem 2 shows that the
SLM needs less local models than the number of hidden nodes in an ELM.
This does not necessarily imply a decrease in computation time, since the
matrix of which a pseudo-inverse needs to be computed, has size N ×(n+1)h,
not size N ×h. An ELM with h hidden nodes and a SLM with h/(n+1) local
models are therefore comparable in both approximation capabilities and in
computational efficiency. The construction of the matrix for which a pseudoinverse needs to be computed, however, can be done more efficiently due to
the block structure of K, so there is some gain in computational efficiency
for a large number of models or hidden neurons. There might also exist a
more efficient way to calculate the pseudo-inverse of this block matrix.
It is questionable whether it is desired to use less nodes or models in an
approach that uses random nodes and models, since the model performance
should not depend highly on parameters that are chosen randomly.
Still, there are several advantages of the SLM, compared to an ELM.
As mentioned above, the matrix K can be constructed efficiently because it
contains a block structure, and maybe in the future an efficient method to
calculate the pseudo-inverse of this matrix could be found. Besides this, the
9

SLM is not as black-box as an ELM because the output matrix Γ actually
shows the direct (linear) relations between input and output near several
operating points of the input space, whereas the output weights of an ELM
only show relations between the less interpretable hidden layer and the outputs. Finally, the SLM allows for a clear physical interpretation due to its
similarities with local linear models from the area of system identification,
and its similarities with fuzzy inference models; see Sections 5.3 and 5.4.
4.4. Proposed algorithm
The proposed SLM learning algorithm is shown in Figure 1. After running
the output of the SLM for input x is given by SLM (x) =
 the algorithm,

g1 (x)z
 
x

.
T 
..
Γ 
.
 , with z =
1
gh̃ (x)z

SLM Algorithm
1. Given N input-output pairs (xj , yj∗ ) ∈ Rn × Rm , randomly generate parameters ai ∈ Rn and bi ∈ R>0 from a continuous
probability distribution independent of the input.
2
2. Choose the number of models h and let gi (x) = e−bi ||x−ai || for
i = 1, . . . 
, h. 
xj
3. For zj =
, calculate the N × (n + 1)h-matrix
1


g1 (x1 )z1T . . . gh (x1 )z1T


..
..
K=
.
.
.
T
T
g1 (xN )zN . . . gh (xN )zN
∗ T
4. Put the desired outputs yj∗ in matrix Y ∗ = [y1∗ . . . yN
] and calculate Γ = K † Y ∗ by using the Moore-Penrose pseudo-inverse.

Figure 1: Proposed SLM algorithm.

10

5. Interpretation of the model
Several interpretations of the SLM are possible. These interpretations
will be described in the following subsections.
5.1. ELM interpretation
The main interpretation of the SLM in this paper is the extension of an
ELM with RBF nodes. The output weights βi of an ELM can be seen as
constant functions, since they do not depend directly on the input. After
training an ELM, the output weights are the same for every possible input
of the network. In the SLM extension, the output weights are not constants,
but functions of the input. Instead of multiplying the outputs of the hidden
layer gi (x) with constant output weights βi , they are multiplied with a linear
function of the input: αi x+βi . This also changes the interpretation of output
weight to linear model, and from hidden layer output to model weight.
Since only the parameters in the hidden layer are fixed, a SLM has more
free parameters than an ELM. The vector βi is of size m and the matrix αi is
of size m×n. If h1 is the number of local models in a SLM and h2 the number
of hidden neurons in an ELM, this implies that the SLM has h1 · m · (n + 1)
free parameters, while an ELM has only h2 · m free parameters. A logical
conclusion is that an ELM needs n + 1 times as many neurons as the number
of local models in a SLM, for similar approximation capabilities. This is true
under the conditions shown in Theorem 2.
5.2. RBF interpretation
A SLM can be interpreted as a RBF network, but there are some important differences. In the first place, it is a generalisation of the RBF network.
Several generalisations of RBF networks exist [9, 10], but the one used in this
paper is similar to the one proposed by Hunt, Haas and Murray-Smith [11].
Using the notation
from this paper, a standard RBF network can be deP
noted as y = hi=1 βi gi (x), while a generalised RBF network is denoted as
P
y = hi=1 θi (x)gi (x), letting the output weights depend on the input 1 .
Using θi (x) = αi x + βi , where the output weights are a linear function of the
input, this is equivalent to a SLM.
1

The generalised RBF as defined by Hunt, Haas and Murray-Smith, differs from a
standard RBF network in three aspects. Only one of these aspects is considered here,
namely that the output weights can depend on the input.

11

In the second place, a SLM differs from a RBF network in the way the
hidden layer parameters are chosen. It is customary in RBF networks to
let the location of the centers of the Gaussians ai be the same as some of
the training samples xj , or else to determine them by some sort of clustering algorithm. In a SLM, both the centers ai as well as the widths bi of
the Gaussians are taken from a continuous probability distribution, and are
therefore independent of the training samples 2 . This has the advantage that
no prior information about the training samples is necessary, and that it
takes practically no time to determine the parameters.
5.3. TSM interpretation
The generalised RBF network has been shown to be equivalent to the
Takagi-Sugeno fuzzy model [4] under some circumstances [11]. The TakagiSugeno model (TSM) consists of a number of fuzzy if-then rules with a fuzzy
premise and a linear consequence, e.g.
Ri : if x1 is A1 ∧ x2 is A2 ∧ · · · ∧ xm is Am
then θi (x) = αi x + βi

(20)

The circumstances under which the TSM is equivalent to the generalised
RBF network, and therefore also to the SLM, are: 1) the number of fuzzy
rules is equal to the number of RBF units, or local models for the SLM,
2) the membership functions within each rule are all Gaussians, and 3) the
operator used for the fuzzy ∧ is multiplication. Using multiplication for the
if-then part as well, and letting gik (x) denote the fuzzy membership function
of the fuzzy set Ak of fuzzy rule Ri , this result is straightforward:
T SM (x) =

h Y
n
X

gik (xk )(αi x + βi )

i=1 k=1

=

h Y
n
X

2

e−bi (xk −aik ) (αi x + βi )

i=1 k=1

2

The probability distribution of the parameters ai and bi can be any continuous probability distribution, for example the normal or the uniform distribution. However, the
distribution of bi should be such that P (bi ≤ 0) = 0, since bi ∈ R>0 . In theory this is the
only constraint, but in practice there might be more constraints due to rounding errors in
the numerical evaluations of the Gaussians.

12

=

=

h
X
i=1
h
X

e−bi

Pn

k=1

(xk −aik )2

(αi x + βi )

(21)

2

e−bi ||x−ai || (αi x + βi )

i=1

= SLM (x).
A SLM can therefore be seen as a TSM with randomly generated Gaussian
membership functions for the premises of the fuzzy rules. This equivalence
relation gives the following result:
Corollary 1. Theorem 2 also holds for the Takagi-Sugeno fuzzy model with
h fuzzy rules, if the membership functions for each rule are Gaussians with
parameters chosen randomly from continuous probability distributions, and
the operator used for the fuzzy ∧ is multiplication.
5.4. LLM interpretation
LLMs like in definition 2 usually appear in the context of nonlinear dynamical systems, where a complex (nonlinear) system is approximated by
taking a weighted sum of simpler (linear) models. These linear models can
be seen as linearisations of the nonlinear system around an operating point,
and the weights are used to describe the transitions between the operating
points. The problem of using a LLM to approximate a nonlinear dynamical
system, belongs to the area of system identification. Several overviews of
identification techniques for local models exist [12, 13].
The input of such a dynamical LLM consists of past inputs and outputs
or states of the system, while the output of the dynamical LLM should give
the next output or state of the system, e.g.
st+1 = f (st , ut ),

(22)

where s is the output or state, and u the control input of the dynamical
system. If the function f is differentiable, it can be linearised around an
operating point (s(i) , u(i) ), to get:
f (st , ut ) ≈ f (s(i) , u(i) ) +

∂f (i) (i)
(s , u )(st − s(i) )
∂st

∂f (i) (i)
(s , u )(ut − u(i) )
∂ut
= Ai st + Bi ut + Oi ,
+

13

(23)

where Ai and Bi are the partial derivatives of f w.r.t. s and u respectively,
and Oi = f (s(i) , u(i) ) − Ai s(i) − Bi u(i) . The local linear model is then a
weighted sum of these linearisations:
LLM (st , ut ) =

h
X

gi (st , ut )(Ai st + Bi ut + Oi ),

(24)

i=1

where gi is a scalar-valued function that represents the validity of the local
model. The relationship between these kinds of dynamical LLMs and TSMs
has also been investigated [14].
Going out of the area of dynamical systems, if the input of the function
f is denoted as x rather than (st , ut ), the LLM of equation 24 is equal to the
one from Definition 2.
Now, for a SLM, the weights are Gaussian functions of the input, but
the parameters of the Gaussians (i.e. their centers and widths, or the operating points and regimes) are chosen randomly. This implies that the local
models are supposed to be good approximations of the global system around
randomly chosen operating points, while the traditional approach is to use
a gradient-descent type algorithm to determine the operating points and
regimes. Theorem 2 shows that even with randomly chosen operating points
and regimes, it is possible to approximate any function with any desired accuracy, as long as the number of local models h is high enough, and that this
can be done with a direct algorithm.
6. Experiments and results
As an example, in this section the SLM algorithm is compared with the
ELM algorithm for modelling the following dynamical system:
dx1
= x2 ,
dt
dx2
= λ(x21 − 1)x2 − x1 .
dt

(25)

These are the equations of a Van der Pol oscillator. It is a nonlinear system with a limit cycle, a phenomenon that does not occur in linear systems.
The experiment is implemented in Matlab, on a computer with a 2.80 GHz
processor.
14

6.1. Data
The model inputs are the x1 and x2 values of the Van der Pol oscillator
with λ = 1 for timesteps t = 1, . . . , 1000, while the desired model outputs are
the x1 and x2 values for timesteps t + 1. The values are computed with the
Euler forward method. The data is seperated in three phases: the learning
phase, the generalisation phase, and the simulation phase. In each phase, the
dynamical system is run for 10 different initial conditions, giving 104 data
samples.
6.2. Models used
The data is tested on a SLM with 100 local models. The parameters ai
are drawn from a zero-mean normal distribution with variance 2, while the
parameters bi are drawn from a uniform distribution over (0, 1). The ELM
uses the same distributions for the parameters, but with 300 hidden neurons
since h(n + 1) = 300 in this case (see Section 5.1 for a comparison between
the number of hidden neurons in an ELM and the number of local models in
a SLM).
The experiment consists of three phases: a learning phase, a generalisation phase and a simulation phase. In the learning phase, the SLM is trained
using the algorithm from Figure 1, and the ELM is trained using the ELM
algorithm proposed by Huang et al. [2]. In the generalisation phase, only the
inputs of the system are used, and the models have to produce the outputs
using the model parameters obtained during the training phase. In the simulation phase, only the initial conditions of the system are given, and the
models have to use the model output of the previous timestep as an input
for future timesteps, also using the model parameters obtained during the
training phase.
6.3. Results and discussion
The whole process has been repeated 100 times. The average results are
shown in Table 1, where MSE stands for Mean Squared Error. The output
of the training phase of the last run is shown in Figure 2, for the SLM case
only. The ELM case showed similar results, as well as the generalisation and
simulation phases for both algorithms.
From the results we see that the SLM algorithm takes less computation
time. As expected, this is not due to the calculation of the pseudo-inverse,
which has actually become slower in the SLM algorithm (the reason for this is
15

unclear). The construction of the matrix that is to be inverted has been done
more efficiently in the SLM approach by making use of its block structure.
The errors during all three phases are similar in both methods, while the
number of local models for the SLM was three times less than the number
of hidden neurons for the ELM. This is in line with the theory presented in
this paper.

Table 1: Comparing the errors and computation time of the SLM and ELM algorithm for
the modelling of a Van der Pol oscillator.

Total
computation
time
Computation time for
pseudo-inverse
MSE during training
MSE during generalisation
MSE during simulation

Average over
100
SLM
tests
9.0952 sec.

Average over
100
ELM
tests
17.8737 sec.

1.8038 sec.

1.5514 sec.

2.3872 · 10−11
2.2645 · 10−10

1.3108 · 10−11
3.1529 · 10−9

7.6113 · 10−5

7.5266 · 10−5

7. Conclusions
An extension to Extreme Learning Machines (ELM) was proposed in this
paper, together with an algorithm and theorems concerning the approximation capabilities of the Random Local Linear Model (SLM). Under certain
conditions, the SLM can achieve similar approximation capabilities as the
ELM while using less local models than the number of hidden nodes of an
ELM. Although this is not necessarily an advantage, the SLM allows for a
clear physical interpretation and can be run more efficiently than an ELM
when the number of local models and hidden nodes is high.
The SLM is similar to Local Linear Models from the area of system identification and, under some circumstances, to the well-known Takagi-Sugeno
Model from the area of fuzzy modelling, for which an upper bound on the
number of randomly generated rules is given. This also makes them more fit
16

Figure 2: Modelling a Van der Pol oscillator with a Random Local Linear Model. The
black dots are the locations of the centers of the Gaussians (parameters ai ), which are
distributed randomly. Around each black dot, a certain linear model is valid. Adding
the linear models together in a nonlinear fashion results in a global model that accurately
approximates the nonlinear dynamics of the oscillator.

to be used for modelling dynamical systems. Compared to these techniques,
the training of a SLM can be done more efficiently, since the operating points
or fuzzy rules can be chosen randomly and independent from the input, which
implies that less prior knowledge is required. The training algorithm then
uses linear regression, which can be performed much faster than most other
parameter estimation methods.
The theorems and proofs of this paper could be extended to more general
cases. However, the link with system identification and fuzzy modelling
techniques already makes sure that the SLM could be an alternative to some
very popular techniques, and it allows much room for further investigation.

17

Acknowledgment
This research was performed at Almende B.V. in Rotterdam, the Netherlands. The author would like to thank Anne van Rossum and Dimitri Jeltsema for their supervision during this research, and Giovanni Pazienza for
providing valuable feedback on this paper.
References
[1] O. Nelles, Nonlinear system identification: from classical approaches to
neural networks and fuzzy models, Springer, 2001.
[2] G.-B. Huang, Q.-Y. Zhu, C.-K. Siew, Extreme learning machine: theory
and applications, Neurocomputing 70 (1) (2006) 489–501.
[3] G.-B. Huang, L. Chen, C.-K. Siew, Universal approximation using incremental constructive feedforward networks with random hidden nodes,
Neural Networks, IEEE Transactions on 17 (4) (2006) 879–892.
[4] T. Takagi, M. Sugeno, Fuzzy identification of systems and its applications to modeling and control, Systems, Man and Cybernetics, IEEE
Transactions on (1) (1985) 116–132.
[5] G.-B. Huang, L. Chen, Convex incremental extreme learning machine,
Neurocomputing 70 (16) (2007) 3056–3062.
[6] C. R. Rao, S. K. Mitra, Generalized inverse of a matrix and its applications, J. Wiley, New York.
[7] G.-B. Huang, H. A. Babri, Upper bounds on the number of hidden
neurons in feedforward networks with arbitrary bounded nonlinear activation functions, Neural Networks, IEEE Transactions on 9 (1) (1998)
224–229.
[8] S. Tamura, M. Tateishi, Capabilities of a four-layered feedforward neural
network: four layers versus three, Neural Networks, IEEE Transactions
on 8 (2) (1997) 251–255.
[9] S. A. Billings, H.-L. Wei, M. A. Balikhin, Generalized multiscale radial
basis function networks, Neural Networks 20 (10) (2007) 1081–1094.

18

[10] F. Fernández-Navarro, C. Hervás-Martı́nez, J. Sanchez-Monedero, P. A.
Gutiérrez, Melm-grbf: a modified version of the extreme learning machine for generalized radial basis function neural networks, Neurocomputing 74 (16) (2011) 2502–2510.
[11] K. J. Hunt, R. Haas, R. Murray-Smith, Extending the functional equivalence of radial basis function networks and fuzzy inference systems,
Neural Networks, IEEE Transactions on 7 (3) (1996) 776–781.
[12] R. Murray-Smith, T. A. Johansen, Multiple model approaches to nonlinear modeling and control (1997).
[13] V. Verdult, Non linear system identification: a state-space approach,
Twente University Press, 2002.
[14] T. A. Johansen, R. Shorten, R. Murray-Smith, On the interpretation
and identification of dynamic takagi-sugeno fuzzy models, Fuzzy Systems, IEEE Transactions on 8 (3) (2000) 297–313.
Laurens Bliek received the B.Sc. degree in applied mathematics at Delft
University of Technology in 2011, and is currently graduating for the M.Sc.
degree of the same study. His minor, B.Sc thesis, and many of his electives,
had to do with artificial intelligence and neural networks. His research interests include neural networks, dynamical systems theory and optimisation.

19

