Return to GeoComputation 99 Index

Parameter estimation in neural spatial interaction modelling by a derivative free global optimization method

Manfred M. Fischer and Martin Reismann
Department of Economic and Social Geography, Wirtschaftsuniversitaet Wien, Augasse 2-6, A-1090 Vienna, Austria
E-mail: Manfred.M.Fischer@wu-wien.ac.at
E-mail: reismann@wigeo1.wu-wien.ac.at

Katerina Hlavácková-Schindler
Institute of Computer Science, Academy of Sciences of the Czech Republic,
Pod Vodárenskou vezí 2, 18207 Praha 8, Czech Republic
E-mail: katka@uivt.cas.cz

Abstract

Parameter estimation is one of the central issues in neural spatial interaction modelling. Current practice is dominated by gradient-based local minimization techniques. They find local minima efficiently and work best in unimodal minimization problems, but can get trapped in multimodal problems. Global search procedures provide an alternative optimization scheme that allows it to escape from local minima. Differential Evolution has been recently introduced as an efficient direct search method for optimizing real-valued multimodal objective functions (Storn and Price 1996). The method is conceptually simple and attractive, but little is known about its behaviour in real world applications. This paper explores this method as an alternative to current practice for solving the parameter estimation task, and attempts to assess its robustness, measured in terms of in-sample and out-of-sample performance. A benchmark comparison with back propagation of conjugate gradients illustrates the superiority of Differential Evolution. 

1. Introduction

The development of spatial interaction models is one of the major intellectual achievements and, at the same time, perhaps the most useful contribution of spatial analysis to social science literature. Since the pioneering work of Wilson (1970) on entropy maximization, there have been surprisingly few innovations in the design of spatial interaction models. The competing destinations version of Fotheringham (1983), the application of feed-forward neural networks (Openshaw 1993, Fischer and Gopal 1994), the use of genetic algorithms and genetic programming to breed new forms of spatial interaction models (Openshaw 1988, Diplock 1996, Turton, Openshaw and Diplock 1997) and neural network topology optimization (Fischer and Leung 1998) are the principal exceptions.

Neural spatial interaction models are termed neural in the sense that they are based on neural computational models, inspired by neuroscience. They are closely related to spatial interaction models of the gravity type, and can be understood as a special class of general feed-forward neural network models with a single hidden layer and sigmoidal transfer functions (Fischer 1998). This class of networks can provide approximation within an arbitrary precision (i.e., it has universal approximation property), as proven by Hornik et al. (1989).

Learning from examples, the problem for which feed-forward neural networks were designed in order to solve, is one of the most important research topics in artificial intelligence. A possible way to formalize learning from examples is to assume the existence of a function representing the set of examples, thus, enabling it to generalize. This can be called a function reconstruction from sparse data (or in mathematical terms, depending on the required precision, approximation or interpolation problem, respectively). Within this general framework, the main issues of interest are the representational capabilities of a given network model and the procedures for obtaining the optimal network parameters. In this contribution, the second issue, network training (i.e., parameter estimation), will be addressed.

Many training (learning) methods find their roots in function minimization algorithms, which can be classified as local or global minimization algorithms. Local minimization algorithms, such as the gradient descent and the conjugate gradient methods (for more details, see e.g., Fischer and Staufer 1999), are fast, but usually converge to local minima. In contrast, global minimization algorithms give better results, because in principle they can escape from local error minima. The computational cost for a single run of a global optimization routine is usually high (Fischer, Hlavackova-Schindler and Reismann 1999).

The next section provides a summarized description of single hidden layer neural spatial interaction models (Section 2). The parameter estimation problem is defined in Section 3 as a problem of minimizing the sum-of-squares error function along with a brief characterization of local search procedures that are generally used to solve this problem. Differential Evolution is outlined in some detail in Section 4. The testbed used for the evaluation of this parameter estimation procedure is Austrian interegional tele-communication traffic data (see Fischer and Gopal 1994), because this data set is known to pose a difficult problem to neural networks using gradient-based learning due to multiple local minima. Section 5 reports on a set of experimental tests carried out to identify an optimal parameter setting and assesses the efficacy of the DEM approach compared to multistart of back propagation of conjugate gradients. Section 6 summarizes the results achieved. 

2. The class of neural spatial interaction models under consideration

Suppose we are interested in approximating an N-dimensional spatial interaction function , where as N-dimensional Euclidean real space is the input space and as 1-dimensional Euclidean real space is the output space. This function should estimate spatial interaction flows from regions of origin to regions of destination. (In practice, only bounded subsets of the spaces are considered.) The function F is not explicitly known, but given by a finite set of samples S = {(x^k, y^k), k = 1,..., K}, so that F(x^k) = y^k, k = 1,..., K. The set S is the set of pairs of input and output vectors. The task is to find a continuous function which approximates (or interpolates) set S. In real world applications, K is a small number and the samples contain noise.

To approximate F, we consider the class of neural spatial interaction models with one hidden layer, N input units, J hidden units and one output unit. consists of a composition of transfer functions so that the (single) output y of is

.

Vector x = (x₀,..., x_N) is the input vector augmented with a bias signal x₀ which can be thought as being generated by a 'dummy unit' whose output is clamped at 1. The w_jn's represent input to hidden connection weights and the w_j's hidden to output weights (including the biases). The symbol w is a convenient shorthand notation of the d = (J · (N + 1) + J + 1)-dimensional vector of all the w_jn and w_j network weights and biases (i.e. model parameters). and are differentiable non-linear transfer functions of the hidden units j = 1,..., J and the output unit, respectively. Following Fischer and Gopal (1994), we will consider only the case N = 3, i.e., the input space will be a closed interval of the three dimensional Euclidean space . The three input units correspond to the independent variables of the classical unconstrained spatial interaction model of the gravity type. They represent measures of origin propulsiveness, destination attractiveness, and spatial separation. The output unit corresponds to the dependent variable of the classical model and represents the spatial interaction flows from origin to destination.

One of the major issues in neural spatial interaction modelling includes the problem of selecting an appropriate member of model class in view of a particular real world application. This model specification problem includes both the choice of appropriate transfer functions and and the determination of an adequate network topology of (i.e., the number of hidden units J). Clearly, the model choice problem and the parameter estimation problem, i.e., the determination of an optimal set of model parameters, are intertwined in the sense that if a good model specification can be found, the success of which depends on the particular problem, then the step of parameter estimation (also termed network training), may become easier to perform.

In this contribution we focus only on the parameter estimation problem. Without loss of generality, we assume the transfer functions for all j = 1,..., J, and equal to the logistic function, thus, consider the special class of functions :

.

The approximation of then merely depends on the learning samples S, and the learning (training) algorithm that determines the parameter vector w from S and the number J of hidden units. 

3. The parameter estimation problem and local minimization procedures

Without loss of generality, we assume the neural spatial interaction model to have fixed topology, i.e., J is predetermined. Then, the goal of learning is to find suitable values w^* for the network weights of the model such that the underlying mapping represented by the training set , is approximated or learned, where k is the index of the training instance. y^k represents the desired network output (i.e. the spatial interaction flows) corresponding to x^k, k = 1,..., K. The process of determining optimal parameter values is called training or learning and can be formulated in terms of minimization of an appropriate error function (or cost function) E to measure the degree of approximation with respect to the actual setting of network weights. The most commonly used error function is the squared-error function of the patterns over the finite set of training data, so that the parameter estimation problem may be defined as the following minimization problem

,

where the minimization parameter is the weight vector w that defines the search space. In this way, the problem of network training has been formulated in terms of minimization of the error function E. This error function is a function of the adaptive model parameters, i.e., network weights and biases. The derivatives of this function with respect to the model parameters can be obtained in a computationally efficient way using the back-propagation technique (see, e.g., Gopal and Fischer 1996; and Fischer and Staufer 1999, for more details). The minimization of continuous differentiable functions of many variables is a problem that has been widely studied, and many of the non-linear minimization algorithms available are directly applicable to the training of neural spatial interaction models as described in Section 2. The general scheme of these algorithms can be formulated as follows:

(i) Choose an initial vector w in parameter space and set t = 1;

(ii) Determine a search direction d(t) and a step size so that

, t = 1, 2,...;

(iii) Update the parameter vector

, t = 1, 2,...;

(iv) If then set t = t + 1 and go to (ii), else return w(t + 1) as the desired minimum.

In this study we refer to the Polak-Ribiere variant of the conjugate gradient (CG) procedure (Press et al. 1992), a mature local optimization method, which is used as a benchmark in Section 5. This algorithm computes the sequence of search directions as

, t = 1, 2,...;

with

,

where is a scalar parameter in the form

for t = 1, 2,...,

w(t - 1)^T is the transpose of w(t - 1). The definition (8) of guarantees that for the sequence of vectors d(t)      holds. The parameter is chosen to minimize

for t = 1, 2,...

in the t-th iteration. This gives the automatic procedure for setting the step length, once the search direction d(t) has been determined.

This and other local optimization methods tend to have difficulties, when the surface is flat (i.e., gradient is close to zero), when gradients are in a large range, and when the surface is very rugged. When gradients vary greatly, the search may progress too slowly, when the gradient is small and may overshoot where the gradient is large. When the error surface is rugged, a local search from a random starting point generally converges to a local minimum close to the initial point with a worse solution than the global minimum. 

4. A new global search approach: the differential evolution method

Global search algorithms employ heuristics to allow escape from local minima. These algorithms can be classified as probabilistic or deterministic. Of the few deterministic global minimization methods developed, most apply deterministic heuristics to bring search out of a local minimum. Other methods, such as covering methods, recursively partition the search space into subspaces before searching. None of these methods operate well or provide adequate coverage when the search space is large as is usually the case in neural spatial interaction modelling.

Probabilistic global minimization methods rely on probability to make decisions. The simplest probabilistic algorithm uses restarts to bring a search out of a local minimum when little improvement can be made locally. More advanced methods rely on probability to indicate whether a search should ascend from a local minimum: simulated annealing, for example, when it accepts uphill movements. Other probabilistic algorithms rely on probability to decide which intermediate points to interpolate as new trial parameter vectors: random recombinations and mutations in evolutionary algorithms (see for example, Fischer and Leung 1998).

Central to global search procedures is a strategy that generates variations of the parameter vectors. Once a variation is generated, a decision has to be made whether or not to accept the newly derived trial parameter. Standard direct search methods (with few exceptions such as simulated annealing) use the greedy criterion to make the decision. Under this criterion, a new parameter vector is accepted if and only if it reduces the value of the error function. Although this decision process converges relatively fast, it has the risk of entrapment in a local minimum. Some stochastic search algorithms, such as genetic algorithms, and evolution strategies, employ a multipoint search strategy in order to escape from local minima.

The Differential Evolution Method (DEM), originally developed by Storn and Price (1996, 1997) is a global optimization algorithm that employs a structured, yet randomized parallel, multipoint search strategy, which is biased towards reinforcing search points at which the error function E(w) being minimized has relatively low values. DEM is similar to simulated annealing in that it employs a random (probabilistic) strategy. One of the apparent distinguishing features of DEM is its effective implementation of the parallel multipoint search. DEM maintains a collection of samples from the search space rather than a single point. This collection of samples is called population of trial solutions.

To start the stochastic multipoint search, an initial population P of, say M, d-dimensional parameter vectors P(0) = {w₀(0),..., w_M-1(0)} is created. For the parameter estimation problem at hand d = (5 · J) + 1. Usually this initial population is created randomly because it is not known a priori, where the globally optimal parameter is likely to be found in the parameter space. If such information is given it may be used to bias the initial population towards the most promising regions of the search space by adding normally distributed random deviations to this a priori given solution candidate. From this initial population, subsequent populations P(1), P(2),..., P(t),... will be computed by a scheme that generates new parameter vectors by adding the weighted difference of two vectors to a third. If the resulting vector yields a lower error function value than a predetermined population member, the newly generated vector will replace the vector that it was compared to, otherwise the old vector is retained. Similarly to evolution strategies, the greedy criterion is used in the iteration process and the probability distribution functions determining vector distributions are not a priori given. The scheme for generating P(t + 1) from P(t) may be summarized by two major stages: the construction of v(t + 1)-vectors from vector w(t)-vectors (stage 1), and the decision criterion whether or not the v(t + 1) - vectors should become members of the population (stage 2) representing possible solutions of the parameter estimation problem under consideration at step t + 1 of the iteration process.

Stage 1: For each population member w_m(t), m = 0, 1,...,.M - 1, a perturbed vector v_m(t + 1) is generated according to

v_m(t + 1) = w_best(t) + Q · (w_r1(t) - w_r2(t))

with r₁, r₂ integers chosen randomly from {0,...,.M -1} and mutually different. These integers also are different from the running index m. Q (0, 2] is a real constant factor that controls the amplification of the differential variation (w_r1(t) - w_r2(t)). The parameter vector w_best(t) which is perturbed to yield v_m(t + 1) is the best parameter vector of population P(t).

Stage 2:The decision whether or not v_m(t + 1) should become a member of P(t + 1), is based on the greedy criterion. If

E(v_m(t + 1)) < E(w_m(t))

then v_m(t + 1) is assigned to w_m(t + 1) otherwise the old value w_m(t) is retained as w_m(t + 1).

The algorithm is stopped after a fixed number of iterations.

5. Performance test results

This section presents results of experiments and describes the performance of the Differential Evolution Method to solve the parameter estimation problem (3) for neural spatial interaction models (2) with three inputs, a single hidden layer with a fixed number of hidden units and a single output unit. The output unit represents the intensity of telecommunication flows from one origin region to a destination region, and the input units of the three independent variables of the classical gravity model: the potential pool of telecommunication activities in the origin region, the potential draw of telecommunication activities in the destination region, and a factor representing the inhibiting effect of geographic separation from the origin to the destination region. The ultimate goal of such neural spatial interaction models is to exhibit good generalization performance, i.e., to make good predictions for new inputs. The spatial interaction modelling prediction accuracy (generalization measured in terms of out-of-sample performance) is generally more important than fast learning.

The model performance is measured in this study by the average relative variance ARV(S) of a set S of patterns given by (Fischer and Gopal 1994)

where y^k denotes the target value and the average over the K desired values in S. ARV(S) provides a normalized mean squared error metric for assessing the in-sample and out-of-sample performance of trained neural spatial interaction models. ARV(S) = 1 if the estimate is equivalent to the mean of the data (i.e, (x^k, w) = ). In the following experiments, the ARV set {ARV1, ARV2} will refer to the average relative error corresponding to the {training set, test set}. 

5.1 The data

The experiments were conducted using Austrian telecommunication flow data (see Fischer and Gopal 1994 for more details). The data set was constructed from three data sources: a (32, 32)-interregional telecommunication flow matrix, a (32, 32)-distance matrix, and gross regional products for the 32 telecommunication regions. It contains 992 4-tuples (x₁, x₂, x₃, y), where the first three components represent the input vector x = (x₁, x₂, x₃) and the last component the target output of the neural spatial interaction model, i.e., the telecommunication intensity from one region of origin to another region of destination. (Input and target output data were preprocessed to logarithmically transformed data scaled into [0, 1].) The telecommunication data stem from network measurements of carried telecommunication traffic in Austria in 1991, in terms of erlang, which is defined as the number of phone calls (including facsimile transmission) multiplied by the average length of the call (transfer) divided by the duration of measurement. This data set was randomly and disjointedly divided into two separate subsets: about two thirds of the data were used for parameter estimation only, and one third as a test set for assessing the generalization performance. In comparison to Fischer and Gopal (1994), the data used were updated to take some measurement errors into account. Consequently, the results of the experimental work described below cannot be directly compared to the previous results. 

5.2 Experimental work and results

To get an optimal solution, it is necessary to use reliable settings for the control parameters Q and M, and an optimal number of hidden units J. There is no other way to a priori define useful combinations of control parameter values than by extensive computational simulation. In order to identify good settings, simulation runs with different combinations of parameter values have been performed. Since every function evaluation has the same complexity, iterations to converge to the optimal ARV2 value were used as a measure of learning time. Simulation runs, which did not converge within 7,000 iterations, were stopped. (Note that population size M multiplied by the number of iterations gives the number of function evaluations.) Each experiment (i.e., a fixed combination of DEM control parameter values) was repeated six (ten) times, the network model being initialized with a different set of random weights from [-0.3, 0.3] before each trial. Simulation runs were realized on PC Pentium 400 Mhz and partly on DEC Alpha 375 Mhz using the egcs-1.0.3 C-compiler.

Table 1 presents the results of the first series of experiments illustrating the effect of the Q-parameter [0.7, 0.8, 0.9, 1.0, 1.1] that controls the amplification of the differential mutation. The cardinality M of the population was set to M = 1000. A fixed topology of 10 hidden units is assumed. In-sample (out-of-sample) performance is measured in terms of ARV1 (ARV2).

There is strong evidence of the robustness of the algorithm (measured in terms of standard deviation) both with respect to the choice of the DEM-parameter Q and to the choice of the initial random state of the simulation runs. Q = 0.9 leads to the best result in terms of average out-of-sample performance (ARV2 = 0.2280). Results for Q = 1.0 are inferior (ARV2 = 0.2440). Even though there is no evidence that Q = 0.9 outperforms other settings such as Q = 0.7, 0.8 and 1.1 in a statistically significant way. Q = 0.9 has been deliberately chosen for the second series of experiments.