Network-Scale Traffic Modeling and Forecasting with Graphical Lasso and Neural Networks
Shiliang Sun Rongqing Huang Ya Gao
cs.LG stat.ML






[cor1]Corresponding author. Tel.: +86-21-54345186; fax:
+86-21-54345119. 
E-mail address: rqhuang09@gmail.com (R. Huang).

Department of Computer Science and Technology, East China Normal University, 500 Dongchuan Road, Shanghai 200241, China



Traffic flow forecasting, especially the short-term case, is an important topic in intelligent transportation systems (ITS).
This paper does a lot of research on network-scale modeling and forecasting of short-term traffic flows. Firstly, we propose
the concepts of single-link and multi-link models of traffic flow forecasting. Secondly, we construct four prediction models
by combining the two models with single-task learning and multi-task learning. The combination of the multi-link model and
multi-task learning not only improves the experimental efficiency but also the prediction accuracy. Moreover, a new multi-link
single-task approach that combines graphical lasso (GL) with neural network (NN) is proposed. GL provides a general methodology
for solving problems involving lots of variables. Using L1 regularization, GL builds a sparse graphical model making use of the
sparse inverse covariance matrix. In addition, Gaussian process regression (GPR) is a classic regression algorithm in Bayesian
machine learning. Although there is wide research on GPR, there are few applications of GPR in traffic flow forecasting. In this paper,
we apply GPR to traffic flow forecasting and show its potential. Through sufficient experiments, we compare all of the proposed approaches
and make an overall assessment at last.


Traffic flow forecasting, graphical lasso (GL), neural network (NN), Gaussian process regression (GPR).






§ INTRODUCTION

With the accelerated pace of modern life, more and more cars come into use. The increased use of cars brings convenience to the public on
the one hand, but on the other hand produces many social problems such as traffic congestion, traffic accidents and environmental
pollution. To make traffic management and traveler information services more efficient, intelligent transportation systems (ITS) emerges as
the times required, and the manual traffic management is no longer viable. The benefits of ITS can not be realized without the ability to
anticipate short-term traffic conditions. As an important aspect of traffic conditions, traffic flows can give insights into traffic conditions<cit.>.
Therefore, short-term traffic flow forecasting becomes one of the most important and fundamental problems in ITS. Short-term traffic flow forecasting is to
determine the traffic flows in the next time interval, usually in the range of five minutes to half an hour, using historical data<cit.>.
A good short-term traffic flow forecasting model can tell the right traffic condition in the near future and make traffic management more effective in turn.
In recent years, ITS, especially short-term traffic flow forecasting has already attracted great interest of researchers. Our work also focuses on the topic
of short-term traffic flow forecasting.


The items detected by ITS generally include traffic flow, volume and occupancy etc. Among all these items, traffic flow is considered to
be the typical metric of traffic condition on a certain link <cit.>. Traffic flow measures the number of vehicles passed
through in a defined time interval and lower traffic flow means heavier traffic congestion. Traditional traffic flow forecasting predicts
a future flow of a certain link only using the historical data on the same link, which is also called single-link traffic flow forecasting.
Obviously, single-link forecasting approaches ignore the relations between the measured link and its adjacent links. In fact, each link is
closely related to other links in the whole transportation system, especially their adjacent links. In this paper, we put forward the multi-link
forecasting models which take the relations between adjacent links into account. Sufficient experiments on real world data show that multi-link
approaches are superior to single-link approaches.

In the past decade, series of traffic flow forecasting approaches have been proposed,
such as time series based approaches <cit.>,
nonparametric methods <cit.>, local regression models <cit.>,
neural network approaches <cit.>, Kalman filtering <cit.>, Markov chain model <cit.> and so on.
Among all these approaches, neural-network-based forecasting approaches are considered as relatively effective methods due to their well-established models.
Typical neural-network-based forecasting methods mainly include back propagation (BP) neural network <cit.>, radial basis
function (RBF) neural network <cit.>, recurrent neural network <cit.>, time delayed neural network <cit.>, resource allocated networks <cit.>, etc. In this paper, we select BP neural networks to serve as the corresponding neural-network-based experiments.
The competitive results further verify the superiority of the proposed neural-network-based approaches.

Gaussian process regression (GPR) is a classic regression
algorithm basing on Bayesian theory. A Gaussian process is a
generalized Gaussian probability distribution and each process is
specified by its mean function and covariance function <cit.>. Because of
the characteristics of easy implementation, few parameters and
strong interpretability, GPR is studied widely in machine learning.
Furthermore, theoretical and practical developments over the last
decade have shown that Gaussian process is a serious competitor for
supervised learning applications <cit.>. However, there are few
applications of GPR in traffic flow forecasting. In this paper, we give a brief analysis of GPR and apply it to
traffic flow forecasting. Through sufficient tests of GPR in real-world data sets, we
point out the potential of GPR for traffic flow forecasting.

Graphical model is not rare in both statistics and computer science.
It is considered as an intersection of the two fields. In
statistics applications, there are often large-scale models with
thousands or even millions of variables involved. Similarly,
there are the same problems in machine learning applications, such
as biological information retrieval, language processing and so on.
Graphical lasso (GL) provides a general methodology for solving such
problems <cit.>. By using L1 regularization, GL
builds a sparse graphical model making use of the sparse inverse
covariance matrix. In this paper, we provide a detailed discussion of the GL algorithm in theory and apply it to multi-link traffic flow
forecasting. With the further information extracted by GL, and combining with the BP neural networks,
we construct a new multi-link single-task traffic flow prediction model, which we refer to as GL_NN.

Parts of our work have been presented recently at international conferences <cit.>.
In this paper, we combine and extend them to give a more systematical analysis.
The remainder of this paper is organized as follows. First, we introduce the four prediction models
basing on NNs. Next, we give the introduction of GPR and GL, respectively, which are closely
related to our work. Then, all the corresponding experiments and discussions on GPR are
presented in the section of experiments.  Finally, conclusions are given in the last section.



§ PREDICTION MODELS WITH NEURAL NETWORKS

Due to the excellent ability in handling complex problems, and the
characteristics of self-learning, self-organizing and self-adaptation,
neural networks (NNs) usually perform well in machine learning problems. On the other hand,
multi-task learning (MTL) is widely applied in computational intelligence and has shown competitive performance.
The main difference between MTL and single-task learning (STL) can be that, with the same inputs, MTL has multiple outputs but STL
has only one output at a time. As to the multiple tasks in MTL, there is only one main task and the others are extra tasks assisting
the learning of the main task. More details about MTL and STL can be found in <cit.>.
In this paper, we further combine the single-link and multi-link models with single-task learning and multi-task learning to
construct four prediction models. The four models are single-link single-task learning (SSTL), single-link
multi-task learning (SMTL), multi-link single-task learning (MSTL)
and multi-link multi-task learning (MMTL), respectively.



 §.§ Single-Link Model

Traditional traffic flow prediction models are single-link models, which predict the future flow of one certain road link using only the historical data of the same link. Combining the single-link model with single-task learning and multi-task learning, we construct two models which are SSTL and SMTL,  respectively. The main difference of the two models lies in the different numbers of outputs, which is also the difference of STL and MTL approaches in a very narrow sense. Following the settings in <cit.>, we set the number of outputs as 3(one main task and
two extra tasks) for our MTL approaches. For one link, we use the first 5 historical traffic flows to predict the next one. That is, the number of inputs is 5.




Take link Ba in Fig. 1 as an example. Record the traffic flow of road link Ba at time interval n as t_a(n), and then
the corresponding five historical traffic flows are respectively t_a(n-5),⋯, t_a(n-1). In the single-link model, we predict t_a(n) using
t_a(n-5), ⋯, t_a(n-1). Basing on NNs, the five historical flows are five inputs. In SSTL, t_a(n) is the one and only output.
While in SMTL, there are three outputs t_a(n-1), t_a(n) and t_a(n+1). Among the three outputs, t_a(n) is the main task,
t_a(n-1) and t_a(n+1) are extra tasks to assist the prediction of the main task. Note that the selection of extra tasks are not
specified. Here we follow the settings in previous experiments <cit.>. Diagrams of the two single-link models SSTL and SMTL are
shown in Fig. 2a and Fig. 2b.







 §.§ Multi-Link Model

Obviously, the single-link prediction model is inefficient because it just predicts one future flow of only one link at a time.
Worse still, it does not make use of the relevant information between adjacent links to improve the prediction.
In fact, since vehicles always come from one link and go to other links in the whole transportation system, traffic
flows of all links, especially adjacent links in the whole traffic network are relevant. Therefore, taking the relevance between adjacent links into account,
we combine the multi-link model with single-task learning and multi-task learning to construct the multi-link single-task learning model (MSTL)
and the multi-link multi-task learning model (MMTL). The multi-link model can simultaneously predict multiple traffic flows of multiple links using historical flows of all links
in the same junction at a time. Certainly, using the historical data from multiple links in different junctions to predict the flow on one link of them is the special case of the multi-link model. For example, GL_NN is a special case of the multi-link model.

Again take the map in Fig. 1 as an example. We can see that the junction B connects three links Ba, Bb and Bc. In multi-link
models, we simultaneously predict future flows of the three links, using all the historical data of the three links. Therefore,
similar to the analysis in single-link model, there are 3×5=15 inputs in multi-link models, which are respectively
t_a(n-5), ⋯, t_a(n-1), t_b(n-5), ⋯, t_b(n-1), and t_c(n-5), ⋯, t_c(n-1). Combing with single-task learning and multi-task learning, there
are three outputs t_a(n), t_b(n) and t_c(n) in MSTL, while 3×3=9 outputs t_a(n-1), t_a(n), t_a(n+1), t_b(n-1), t_b(n), t_b(n+1), t_c(n-1), t_c(n) and t_c(n+1) in
MMTL. Diagrams of the two multi-link models are shown in Fig. 3a and Fig. 3b. For the sake of clarity, in Fig. 3b, we draw the three outputs corresponding to the three links Ba,
Bb and Bc each in a box. Through the four diagrams shown as Fig. 2a, Fig. 2b, Fig. 3a and Fig. 3b, we can get a better understanding of single-link and multi-link models.







§ GAUSSIAN PROCESS REGRESSION

Gaussian process regression (GPR) is an important Bayesian machine
learning approach which places a prior distribution over the function space. All inferences of
GPR are taken place in the function space. In supervised learning
applications, it aims at the conditional distribution of the targets
given the inputs but not the distribution of the inputs. Moreover,
Gaussian process is a generalized Gaussian probability distribution.
The random variables in Gaussian probability distribution are
scalars or vectors (multivariate case). In Gaussian process
case, the random processes are represented as functions. There are several ways to
interpret GPR. Here we give a brief inference on the theoretical basis of GPR from the function-space
view. More details can be seen in <cit.>.

Suppose that we have a training set W={(𝐱_i, y_i )| i= 1,
…, n}, where n is the number of observations,
𝐱_i denotes the i-th D-dimensional input variable and
y_i is the corresponding target which is always a real value in the
regression case. Aggregating the inputs as column vectors for all
n case, we can get a D× n design matrix X.
Similarly, collect the targets in a vector as y. Then the
training set W can be written as W= (X, y). In the
same way, we represent a test set M={(𝐱_i^*, y_i^* )| i=
1, …, n^* } as M= (X^*, y^*), where X^* is a
D× n^* matrix.

A Gaussian process is specified by its mean function and
covariance function <cit.>. Define the mean function m(𝐱) and
covariance function k(𝐱,𝐱') of a Gaussian
process f(𝐱) as


    [                               m(x) = E[f(x)],; k(x,x') = E[(f(x) - m(x))(f(x') - m(x'))]   . ]

Then, the Gaussian process can be written as

    f(x) ∼GP(m(x),k(x,x')).


Generally, for the simplicity of notation and computation, the mean
function is set to be zero. In addition, Gaussian process has a
definition that a Gaussian process is a collection of random
variables and any finite number of them have a joint Gaussian
distribution. According to this definition, we can get a consistency
principle. That is, if there is a distribution (y_1, y_2 ) ∼
N(μ, Σ ), then there exists y_1∼ N(μ _1 ,
Σ _11 ), where Σ _11 is the relevant submatrix of
Σ. This consistency principle plays an important role in the
inference of the GPR algorithm, which is also known as the marginalization
property.

Since prediction with noise-free observations is the special case
of prediction with noise observations, we take the noise case to do the inference of GPR. In the noise case, the relation
between the observed target value y and the function value
f(𝐱) is

    y =  f(𝐱) + ε ,

where ε is the additive independent identically
distributed Gaussian noise with variance σ _n^2. In the GPR
inference, the covariance function should be predefined. We take the
squared exponential (SE) covariance function as an example.


    [                    cov (f(x_p ),f(x_q )) = k(x_p ,x_q ); = σ _f exp[  - (x_p  - x_q )^T P^ - 1 (x_p  - x_q )/2], ]

where P is a diagonal matrix with l_1, l_2, …, l_D being the diagonal elements, D is the dimension of the input
space, and σ _f^2 is the signal variance which controls the
global variation. From formula (<ref>), we note that the
covariance between the outputs can be written as a function of the
inputs. Following the independence assumption and with the specified
covariance function, we can easily get the prior of the training set X. That is

    cov (𝐲) = K(X,X) + σ
    _n^2 I.

It is easy to find out that formula (<ref>) can also be seen as the
prior in noise-free case by removing the noise term. According to the
independence assumption, the noise term is a diagonal matrix.

With the prior, we can further get the prior joint distribution of
the observed target values 𝐲 and the function values
𝐟^* of the test samples:

    [ [   𝐲; 𝐟^*;     ]] ∼ N( 0,  [ [ K(X,X) + σ _n^2 I         K(X,X^* );         K(X^* ,X)      K(X^* ,X^* );                   ]]).

To get the posterior distribution over functions, we need to reject
those functions that disagree with observations of the prior. In
probabilistic terms, this operation can be done easily. According
to properties of marginal distribution and conditional distribution <cit.>, we can get

    [                                  f^* | X,y,X.^*  ∼ N(f̅^* ,cov (f^* )),;       f̅^*  ≜ E[f^* | X,y,X.^* ] = K(X^* ,X)[K(X,X) + σ _n^2 I]^ - 1 y,; cov (f^* ) = K(X^* ,X^* )- K(X^* ,X)[K(X,X) + σ _n^2 I]^ - 1 K(X,X^* ).;                                                                         ]

Formula (<ref>) gives the distribution of the function values
𝐟^*. In practical applications, the mean
function value 𝐟̅^* is evaluated as the output of
GPR. In fact, the GPR algorithm simultaneously outputs the variance of
prediction values, which is also considered as the potential
capability of GPR compared to other regression algorithms. A faster and more stable algorithm using Cholesky decomposition to compute
the inverse covariance matrix in formula (<ref>) also can be found in <cit.>.



§ GRAPHICAL LASSO

Graphical lasso (GL) is an algorithm to construct a sparse graphical
model by applying the lasso penalty to the inverse covariance matrix. There is
a basic model assumption that the observations have a multivariate
Gaussian distribution with mean μ and covariance matrix Σ <cit.>. The key to build a sparse graphical model is to make
the inverse covariance matrix as sparse as possible. If the ij-th
component of Σ ^ - 1 is zero, then there is no link between
the two variables i and j in the sparse graphical model.
Otherwise, there exists a link between the two variables. In recent years, series of
approaches have been proposed to solve this problem. All the
approaches can be classified into two types: the approximate approaches and the exact
approaches. The approximate approaches estimate the sparse graphical
model by fitting a lasso model to each variable and using the others
as predictors <cit.>. The exact approaches solve the maximization
of the L1-penalized log-likelihood problem. There are also several ways
to solve the exact problem. For example, interior point
optimization methods <cit.> and blockwise coordinate descent (BCD)
algorithm <cit.>, etc. Thereinto, the BCD algorithm based on GL is appreciated
as a relatively efficient method <cit.>. Below, for completeness, we introduce the GL algorithm. More detailed information
can be found in related references.



 §.§ Problem Setup

Assume that we are given N observations independently drawn from a
p-variate normal Gaussian distribution, with mean μ and
covariance Σ. Let S denote the empirical covariance
matrix. Thus we have

    S = 1/N∑_k = 1^N (x_k  - μ )
    (x_k  - μ )^T ,

where x_k denotes the k-th observation. According to the
independence assumption, we can easily get the likelihood L on the
given data set.

    [                                                       L = ∏_k = 1^N p(x_k ;Σ ); 1pt  = ∏_k = 1^N {1/(2π )^p/2| Σ|^1/2e^ - (x_k  - μ )^T Σ ^ - 1 (x_k  - μ )/2};      = 1/(2π )^Np/2| Σ|^N/2e^ - ∑_k = 1^N (x_k  - μ )^T Σ ^ - 1 (x_k  - μ )/2.;                                                                                ]

Then, the log-likelihood can be written as

    [ log L =  - log [(2π )^Np/2| Σ|^N/2 ] - ∑_k = 1^N (x_k  - μ )^T Σ ^ - 1 (x_k  - μ )/2;     =  - Np/2log 2π  - N/2log| Σ| - 1/2∑_k = 1^N (x_k  - μ )^T Σ ^ - 1 (x_k  - μ ) .;                                                                                      ]

Because the GL algorithm has to solve the problem of maximizing the L1-penalized
log-likelihood, we make a few transformations on formula
(<ref>). Removing the constant term in formula (<ref>) and
combining it with formula (<ref>), we get

    [     log L ∝  - log| Σ| - 1/N∑_k = 1^N (x_k  - μ )^T Σ ^ - 1 (x_k  - μ ); ∝log| Σ ^ - 1| - trace(1/N∑_k = 1^N (x_k  - μ )^T Σ ^ - 1 (x_k  - μ ) ); ∝log| Σ ^ - 1| - trace(Σ ^ - 1·1/N∑_k = 1^N (x_k  - μ ) (x_k  - μ )^T );                                      ∝log| Σ ^ - 1| - trace(Σ ^ - 1· S);                                    ∝log| Σ ^ - 1| - trace(S ·Σ ^ - 1 ).;                                                                         ]

Therefore, the exact problem that the GL algorithm solves can be written as

    Σ ^ - 1  = max_X ≻
    0log| X | - trace(S · X) - ρ X
    _1,

where X _1 is the L1-norm of matrix X, which is
the sum of the absolute values of the elements of X, and ρ is
the penalty parameter which controls the extent of
penalization <cit.>.



 §.§ Formula Transformations

Focusing on formula (<ref>), a series of
transformations are carried out to get an equivalent form which can be easily solved. Firstly,
we write the problem as

    max_X ≻ 0min_ U _∞⩽ρlog| X
    | + trace(X,S + U),

where U _∞ denotes the maximum absolute value element
of the symmetric matrix U <cit.>. Exchanging the max and min, formula (<ref>) is transformed as follows.

    maxlog| X | - trace(X,S + U).

Computing the derivative of formula (<ref>) over X, we
obtain X=(S+U)^-1. Replace X in formula (<ref>) with (S+U)^-1 and
it follows that the dual problem of formula (<ref>) becomes

    min_ U _∞⩽ρ
    - log| (S + U) | - p,

where p is the dimension of matrix X, and the relation between the primal and the dual variables
is X=(S+U)^-1. To write neatly, set M=S+U. Then, the dual
problem of the primal maximum L1-penalized log-likelihood problem is

    Σ  = max{log| M |:M - S_∞⩽ρ}.

According to the series of transformations above, we find that we finally
estimate Σ in the dual problem (<ref>) while the inverse
covariance matrix Σ ^ - 1 in the primal problem (<ref>).
Moreover, we also observe that the diagonal elements of Σ and
S have a relation as

    Σ _ii  = S_ii  + ρ,

which holds for all i.



 §.§ Block Coordinate Descent (BCD) Algorithm

Let W be the estimation of Σ. The block coordinate descent
(BCD) algorithm solves problem (<ref>) by optimizing
cyclically over each row and column of W, until achieving the given
convergence condition. More details about BCD algorithm can be seen
in <cit.>. As to the GL approach we discussed, BCD algorithm plays as a
launching point.

Divide W and S into blocks as

    W = ( [   W_11   w_12; w_12^T   w_22;        ]),S = ( [   S_11   s_12; s_12^T   s_22;        ]).

The BCD algorithm updates w_12 through solving the quadratic  programming


    w_12  = min _y { y^T W_11^ - 1 y:y
    - s_12_∞⩽ρ},

which is solved by an interior point procedure. Permuting the rows
and columns to make the target column always be the last one, BCD
solves the problem like formula (<ref>) for each column
and updates its estimation for W when all columns have been processed. This process is repeated until convergence.

The dual problem of formula (<ref>) shown as the
following is also deduced in <cit.>.

    min _β{1/2W_11^1/2β
    - b^2  + ρβ_1 },

where b = W_11^ - 1/2 s_12. If β solves formula (<ref>), then the solution of formula (<ref>) is
w_12=W_11β. It can be easily found that formula (<ref>)
resembles a lasso regression problem, which is the launching
point of GL approach. There is a verification on the
equivalence between the solutions of formula (<ref>) and (<ref>) given in <cit.>.



 §.§ Algorithm Description and Realization

According to the lasso problem (<ref>) achieved by the BCD algorithm,
the GL approach solves and updates this problem recursively. Details of
the GL algorithm can be described as

1. Set W= S+ ρ I, where I is the identity matrix. Then the
diagonal of W remains unchanged in all the following steps.

2. For each row and column of W, solve the lasso problem (<ref>) and obtain the solution β.

3. Compute w_12 by w_12  = W_11β, and replace the corresponding row and column of W with w_12.

4. Repeat the above steps 2 and 3 until convergence.

5. Compute the inverse matrix of W, which is also the required inverse covariance matrix Σ ^ - 1.

In <cit.>, there also gives a relatively cheap method to execute the
last step of the GL algorithm above. From the achieved
sparse matrix Σ ^ - 1, GL builds the desired sparse
undirected graphical model. Each row or column of matrix Σ ^
- 1 represents a node in the graphical model. The row or column is
corresponding to a variable of the multi-variable data. Therefore,
p-dimensional data has p nodes in the graphical model. Whether
there is a link between two nodes is determined by the
corresponding component of matrix Σ ^ - 1 being zero or not.
If the component is zero, then there is no link in the graphical
model. That is, the two variables are conditionally independent given other
variables. In the next section, combing with the multi-link traffic flow
prediction model, we give an instance of building the sparse graphical
model by GL.




§ EXPERIMENTS



 §.§ Data Description

The data sets used in this paper are vehicle flow rates recorded
every 15 minutes, which were gathered along many road links by the
UTC/SCOOT system of Beijing Traffic Management <cit.>. The unit
of the data is normalized as vehicles per hour (vehs/h). For the
short-term traffic flow forecasting, we carry a one-step prediction
and take 15 minutes as the prediction horizon. That is, we predict
the traffic flow rates of the next 15-minute interval every time.

From the urban traffic map, we select a portion including 31 road
links shown as Fig. 1. Each circle node in the figure represents a
road junction which combines several road links. The arrows show the
directions of traffic flows from the upstream junctions to the
corresponding downstream junctions. Paths without arrows denote no
traffic flow records. Raw data are taken from March 1 to March 31,
2002, totally 31 days. Considering malfunctions of traffic flow
detectors, we wiped away the days with empty data. Finally, the
remaining data we used are of 25 days and have totally 2400 sample points.
We divide the data into two parts, the first 2112 samples as
training data and the rest as test data.




 §.§ Model Building with GL

In multi-link traffic flow prediction case, we can use certain
historical traffic flows of all the links in the whole map. Through
building a sparse graphical model by GL, we extract the informative historical
traffic flows provided by all the links. Basing on the data set
described above, we take 6 continuous traffic flows of each link to
build the sparse graphical model. Because there are 31 links, we
will get a 186×186 inverse variance matrix (6×31=186).
For traffic flow forecasting, the first 5 historical traffic flows of 31 links are all used to
predict the 6-th traffic flow of one link. As to a predicted link, when building the graphical model,
we need the 6 traffic flows on the predicted link and the other 30 links' first 5 historical traffic flows.
Therefore, to a certain link, there are at most 186-30=156 nodes in the graphical model.

Still take link Ba as an example. In the single-link prediction model, we
predict the traffic flow Ba(n) using the continuous historical
traffic flows Ba(n-5), …, Ba(n-1) on link Ba. While in the
multi-link prediction model, we consider the historical traffic
flows of all adjacent links or all the links in the whole traffic
map. The latter one seems to be more comprehensive, but it also brings
too much computation. Fortunately, this problem can be easily solved by
the GL approach. With the sparse graphical model built by GL,
we can extract the most relevant historical traffic flows to predict the predicted traffic flow.
In modeling of link Ba, we just consider the components of the corresponding column
or row in the inverse covariance matrix. If the component is zero, it means there
is no relevance or very little relevance between the two variables. Then there is no link
between the two corresponding variables in the graphical model. For example, let variable i
represent the predicted traffic flow and variable j represent some historical traffic flow.
If there is no link between variables i and j, it means historical traffic flow j contributes nothing to the prediction of traffic flow i.
Fig. 4 gives the sparse graphical model of link Ba built by GL.





By comparing Fig. 1 and Fig. 4 with link Ba, we can see that the prediction of link Ba is not only relevant to the three traffic
flows Ba(t-3), Ba(t-2), Ba(t-1) on link Ba itself but also the other five traffic flows Eb(t-1), Fe(t-2), Fe(t-1) and Hl(t-1), Ib(t-1) of
link Eb, Fe, Hl and Ib respectively. There are only 8 variables considered relevant to the prediction of Ba(t), which is much fewer
than all the 155 variables considered in the general multi-link model. Therefore, GL further extracts the relevant information based on
our previous multi-link prediction model.



 §.§ Experimental Settings

In the design of NNs, a three-layer BP neural network is selected.
On the one hand, a three-layer NN can approximate arbitrary bounded
and continuous functions <cit.>, and on the other hand, more
layers will make the network more complex.
Besides, BP NNs are well known for their good self-learning capability.
The number of input and output units is determined by the dimension of the experimental data.
For example, in single-link models based on NNs, as we use 5 historical traffic flows to predict the traffic flow of the next time interval,
the number of input units is 5, while the number of output units is 1 in the SSTL model and 3 in the SMTL model. The
case of the multi-link model also can be inferred from the representation of the multi-link model in the section of multi-link model. In GL with NN
case, the number of input units depends on the extracted dimension
of the GL algorithm and the number of output units is the same as the SSTL
model. Obviously, different links will have different numbers of input
units in GL with the NN case. For all approaches based on NNs, the number
of hidden units is computed by the empirical formula shown below.

    n = √(n_i  + n_0 )  + a,

where n, n_i and n_0 respectively denote the number of hidden,
input and output-layer units, a is a constant that can be chosen
between 1 and 10<cit.>. To obtain a relatively optimal construction of NN,
we try different values of a from 1 to 10 with an interval 1 and finally
choose the one with best performance. As to the transfer
functions of NN, sigmoid function is selected between the input layer
and the hidden layer, purelin function is selected between the hidden
layer and the output layer. The trainlm function is selected as the
training function, because it is based on Levenberg-Marquardt algorithm
and can converge rapidly with a high prediction accuracy.

In the realization of GPR, we need to specify the covariance
function and find a method to optimize the parameters. In our
experiments, we choose the following squared exponential (SE)
covariance function.


    k(x_p ,x_q ) = σ _f^2 exp[  -
    (𝐱_p  - 𝐱_q )^TP^ - 1
    (𝐱_p  - 𝐱_q )/2],

where P is a diagonal matrix with diagonal elements being l_1, l_2,
…, l_D, D is the dimension of the input space, σ _f^2 is
the signal variance which controls the global variation. Therefore,
together with the noise variance σ _n^2 involved in GPR
algorithm, there are totally D+2 parameters. Following the
suggestion in <cit.>, we initially set l_1, l_2, …, l_D, σ
_f^2 all as 1 and σ _n^2 as 0.1. Then, we optimize these
parameters by maximizing the marginal likelihood

    log p(𝐲|X) =  - 1/2𝐲^𝐓(K + σ _n^2 I)^ - 1𝐲- 1/2log| K + σ _n^2 I| - n/2log 2π ,

which can also be easily achieved by y ∼ N(0,K + σ _n^2
I). We use the gradient descent algorithm to minimize the
negative marginal likelihood to get the optimal parameters.

Similarly, there are two parameters to be specified in GL. One is the penalty parameter ρ, and the other is the
lower limit value which determines the relevance between two
variables. As to the selection of the penalty
parameter, we follow the suggestion in Section 2.3 of <cit.>.
There is also a statement saying that, if we follow the suggestion
to choose the penalty parameter, the error rate of estimating the
graphical model can be controlled. Since the GL algorithm builds the
sparse graphical model according to the inverse covariance matrix,
we need a lower limit value to screen out the nonzero
components of the inverse covariance matrix that represent effective information. In
our experiments with GL, we set the component of the inverse
covariance matrix as zero if it is less than 5e-4. That is, we think that
there is little relevance between the two variables when the
corresponding component in the inverse covariance matrix is so
small.





 §.§ Results

To examine the rationality of the multi-link model, we first compute the correlation coefficients
of adjacent links. There are 10 junctions in the traffic map of Fig. 1, and there are
totally 36 correlation coefficients. We list the 36 correlation coefficients in Table <ref>.

From Table <ref>, we can see that the minimum, the maximum and the mean of the 36 correlation coefficients are
0.5988, 0.9637 and 0.8790, respectively. Thereinto, almost all of the 36 correlation coefficients are larger than 0.8,
which means high correlation between variables. Therefore, as to the real-world data sets we used in the experiments,
it is reasonable and meaningful to test the proposed multi-link approaches on them.

In this paper, the mentioned approaches for traffic flow forecasting include SSTL, SMTL, MSTL, MMTL, GPR and GL_NN, respectively.
We test them on the 31 real-world traffic flow data sets described above, which are collected from the 31 road links of Fig. 1.
To get a complete evaluation of all the proposed approaches, the historical average, which we mark as Hist_Avg, is adopted as a base line for comparison.
We adopt root mean square error (RMSE) and mean absolute relative error (MARE) to evaluate the
prediction performance of different approaches. RMSE and MARE are formulated as follows.

    RMSE=√(1/N∑_i = 0^N (t(i) -t'(i))^2 ) ,

and

    MARE=1/N∑_i = 0^N |t(i) -t'(i)|/t(i)  ,

where t'(i) is the prediction of t(i), and N is the number of test samples.

In Table <ref> and Table <ref>, we present the experimental results of
MARE and RMSE on the 31 road links corresponding to all the compared approaches.




In order to get a comprehensive comparison, we compare the experimental results from two views: the global and the local.
Table <ref> and Table <ref> locally give the MAREs and the RMSEs of the 31 links corresponding to all
the compared approaches. Fig. 5 globally shows the sum of the RMSEs of the 31 links corresponding to all the compared approaches.




Besides these, we further compare all the proposed methods using t-test. The t-test results of all the compared approaches are listed on Table <ref>.


Firstly, we compare all the proposed approaches with the base line Hist_Avg. According to the experimental results shown
in Table <ref>, Table <ref> and Fig. 5, of the 31 data sets, the number of data sets on which SSTL, SMTL, MSTL, MMTL, GPR and GL_NN performs better than Hist_Avg is
19, 21, 22, 22, 24 and 29, respectively, which means that all the proposed traffic flow forecasting approaches are superior to Hist_Avg.
Secondly, we compare the performance of all the proposed approaches with each other such as the comparison of single-link approaches with multi-link approaches and
single-task approaches with multi-task approaches. GPR is a single-link single-task prediction approach in nature, when evaluating its performance, we compare it
with SSTL. As to the comparison criterion MARE, GPR outperforms SSTL on 28 data sets according to Table <ref>. And according to t-test results shown in Table <ref>,
we can find that GPR is significantly better than SSTL. In the GPR column of Table <ref>, we marked the components in bold which show that GPR is better than SSTL. We can see that, in the total
31 links, there are 14 links showing that GPR outperforms SSTL. However, there are another 5 links (italics in the GPR column of Table <ref>) with error rate difference of GPR and SSTL less than 1. In Fig. 5, GPR is globally better than SSTL. Therefore, we still can conclude that GPR outperforms SSTL in traffic flow forecasting. According to the experimental results based on RMSE in Table <ref>, in the columns corresponding to SSTL, SMTL, MSTL, MMTL, we marked the two best components of the four approaches in bold. The numbers of boldfaces
are respectively 11, 13, 17 and 21. Therefore, we can get the conclusion that, multi-link approaches perform better than single-link approaches and multi-task learning approaches are better than single-link learning approaches. This is also why MMTL performs best in the four approaches. GL_NN constitutionally belongs to the multi-link single-task prediction approaches.
The difference between GL_NN and MMTL is that GL_NN extracts the relevant information. From Table <ref>, GL_NN performs better than MMTL on 22 data sets. In the GL_NN column of
Table <ref>, we marked the components in bold which show that GL_NN is better than MMTL. It is easy to find out that there are 21 links of all 31
links showing that GL_NN is better than MMTL. According to the t-test result of GL_NN and MMTL in Table <ref>, it shows that GL_NN outperforms MMTL to some extent.
The results fully verify the superiority of GL used in extracting information through building the sparse graphical model. GL_NN performs best on traffic flow forecasting compared to all
the other proposed approaches.



 §.§ Discussions on GPR

As a single-link single-task approach, GPR performs better than SSTL on traffic flow forecasting.
In this section, we would give special illustrations of GPR on the potential capability.
Actually, the GPR algorithm outputs two terms that are the mean and the variance to be predicted.
Exactly speaking, GPR gives the distribution of the targets rather than the exact values. When computing the prediction errors, we use the mean as the prediction value.
Take link Kd as an example, Fig. 6 shows the practical prediction results. The star curve represents the actual values while the dot curve represents
the prediction values. The shaded area is the fluctuating range of targets predicted by GPR. As we can see, with the shaded part,
more actual targets can be contained in the prediction scope.





Reducing the noise variance in Fig. 6, we can get a new prediction
figure of link Kd shown as Fig. 7. From Fig. 7, we can see that the
fluctuating range gets larger when the noise variance is reduced.
The shaded area in Fig. 7 can even contain all the actual targets.
This is the potential capability of GPR. GPR can get more precise predictions by adjusting more appropriate parameters.
Therefore, in the area that it needs to predict only the output scope rather than precise values, GPR is the approach well worth considering.







§ CONCLUSIONS

Due to the disadvantage of the traditional single-link traffic flow forecasting model,
we propose the multi-link model that predicts traffic flows using historical data from all the adjacent links.
By combining the single-link model and multi-link model with single-task learning and multi-task learning, we propose
four basic traffic flow forecasting approaches, SSTL, SMTL, MSTL and MMTL. Graphical Lasso(GL) is an effective approach in extracting
the relevant information of variables of complex problems by building the sparse graphical model.
We make use of GL to extract the most informative historical flows from all the links in the whole transportation system,
and then construct a BP neural network with the extracted data to predict traffic flows. We refer to the approach combining GL with NN as
GL_NN. GL_NN is aslo a multi-link traffic flow forecasting approach, but it is more efficient than MMTL. The test of GL_NN on real-world traffic
flow forecasting shows competitive results. In addition, we apply GPR to traffic flow forecasting and discuss its potential.
Competitive experimental results reveal the superiority of GL_NN to other proposed approaches.
Moreover, the results further verify that multi-link approaches outperform
single-link approaches and multi-task learning approaches outperform single-task learning approaches in traffic flow forecasting.

In the future, three interesting aspects can be considered. Firstly,
the potential of the multi-link model in traffic flow forecasting
should be further studied. Secondly, GL can be combined with other
approaches not just with NNs. Thirdly, considering that in practical
applications a much lower traffic flow estimation of one road can probably attract vehicles from
adjacent roads and thus causes subsequent traffic difficulties, the
traffic flow prediction methods discussed in this paper can be further
enhanced by investigating microcosmic prediction errors and taking
some precautionary actions for lower estimations.




§ ACKNOWLEDGEMENTS


This work is supported in part by the National Natural Science
Foundation of China under Project 61075005, and the Fundamental
Research Funds for the Central Universities.


99

[Abdulhai et al.(1999)]AbdulhaiNN99
Abdulhai, B., Porwal, H., and Recher, W., 1999. Short term freeway flow
prediction using genetically-optimized time-delay based neural
networks, in: Proceedings of the 78th Annual Meeting of the Transportation Research Board, Washington D.C., USA.

[Abdulhai et al.(2002)]AbdulhaiNN02
Abdulhai, B., Porwal, H., and Recher, W., 2002. Short-term traffic flow prediction using neuro-genetic algorithms,
Journal of Intelligent Transportation Systems, vol. 7, 3–41.

[Banerjee et al.(2008)]Banerjee08
Banerjee, O., Ghaoui, L., and Aspremont, A., 2008. Model selection
through sparse maximum likelihood estimation, Journal of Machine
Learning Research, vol. 9, 485–516.

[Chen and Grant(2001)]ChenSequential01
Chen, H. and Grant-Muller, S., 2001. Use of sequential learning for
short-term traffic flow forecasting, Journal of Transportation
Research Record, Part C: Emerging Technologies, vol. 9(5), 319–336.

[Chen and Chen(2007)]ChenEnsemble07
Chen, L. and Chen, C., 2007. Ensemble learning approach for
freeway short-term traffic flow prediction, in: Proceedings of IEEE
International Conference on System of Systems Engineering, 1–6. 

[Caruana(1997)]Caruana97
Caruana, R., 1997. Multitask learning, in: Proceedings of International
Joint Conference on Machine Learning, vol. 28(1), 41–75.

[Davis and Nihan(1991)]DavisRegression91
Davis, G. and Nihan, N., 1991. Non-parametric regression and
short-term freeway traffic forecasting, Journal of Transportation
Engineering, vol. 177(2), 178–188.

[Davis(1990)]Davis90
Davis, G., 1990. Adaptive forecasting of freeway traffic
congestion, Journal of Transportation Research Record, vol. 1287, 29–33.

[Dahl et al.(2008)]Dahl08
Dahl, J., Vandenberghe, L., and Roychowdhury, V., 2008. Covariancem
selection for non-chordal graphs via chordal embedding,
Optimization Methods and Software, vol. 23(4), 501–520. 

[Duda et al.(2001)]Duda01
Duda, R., Hart, P., and Stork, D., 2001. Pattern
Classification, John Wiley and Sons, New York.

[Friedman et al.(2008)]Friedman08
Friedman, J., Hastie, T., and Tibshirani, R., 2008. Sparse inverse
covariance estimation with the graphical lasso, Biostatistics, vol. 9(3), 432–441.

[Gao and Sun(2010)]Gao2010
Gao, Y. and Sun, S., 2010. Multi-link traffic flow forecasting using
neural networks, in: Proceedings of the Sixth International Conference on Natural
Computation (ICNC), 398–401. 

[Gao et al.(2011)]Gao2011
Gao, Y., Sun, S., and Shi, D., 2011. Network-scale traffic modeling and
forecasting with graphical lasso, in: Proceedings of the Eighth International
Symposium on Neural Networks (ISNN), 151–158. 

[Hall and Mars(1998)]HallArtificial98
Hall, J. and Mars, P., 1998. The limitations of artificial neural
networks for traffic prediction, in: Proceedings of the Third IEEE Symposium on
Computers and Communications, 8–12.

[Jordan(2004)]JordanGraphical04
Jordan, M., 2004. Graphical models, Statistical Science, vol. 19(1), 140–155.

[Jin and Sun(2008)]Jin08
Jin, F. and Sun, S., 2008. Neural network multitask learning for
traffic flow forecasting, in: Proceedings of International Joint
Conference on Neural Networks (IJCNN), 1898–1902.

[Lee and Fambro(1999)]LeeSubsets99
Lee, S. and Fambro, D., 1999. Application of subsets autoregressive
integrated moving average model for short-term freeway traffic
volume forecasting, Journal of Transportation Research Record, vol. 1678, 179–188.

[Moorthy and Ratcliffe(1988)]MoorthyShort88
Moorthy, C. and Ratcliffe, B., 1988. Short term traffic
forecasting using time series methods, Journal of Transportation
Planning and Technology, vol. 12(1), 45–56.

[Meinshausen and Bühlmann(2006)]MeinshausenHighDimensional06
Meinshausen, N. and Bühlmann, P., 2006. High dimensional graphs
and variable selection with the lasso, The Annals of Statistics,
vol. 34, 1436–1462. 

[Okutani and Stephanedes(1984)]OkutaniDynamic84
Okutani, I. and Stephanedes, Y., 1984. Dynamic prediction of
traffic volume through kalman filter theory, Journal of
Transportation Research Record, Part B, vol. 18(1), 1–11.

[Park et al.(1998)]ParkRadialBasis98
Park, B., Messer, C., and Urbanik II, T., 1998. Short-term freeway
traffic volume forecasting using radial basis function neural
network, Journal of Transportation Research Record, vol. 1651, 39–47.

[Rasmussen and Williams(2006)]RasmussenGaussian06
Rasmussen, C. and Williams, K., 2006. Gaussian processes for machine
learning, the MIT press.

[Smith and Demetsky(1997)]SmithTraffic97
Smith, B. and Demetsky, M., 1997. Traffic flow forecasting:
comparison of modeling approaches, Journal of Transportation
Engineering, vol. 123(4), 261–266.

[Smith and Demetsky(1994)]SmithTrafficApproaches94
Smith, B. and Demetsky, M., 1994. Short-term traffic flow prediction
approaches. Prediction: neural network approach, Journal of
Transportation Research Record, vol. 1453, 98–104.

[Sun and Zhang(2007)]Sun07
Sun, S. and Zhang, C., 2007. The selective random subspace predictor
for traffic flow forecasting, IEEE Transactions on Intelligent
Transportation Systems, vol. 8(2), 367–373.

[Ulbricht(1994)]UlbrichtMulti94
Ulbricht, C., 1994. Multi-recurrent networks for traffic
forecasting, in: Proceedings of the Twelfth National Conference on
Artificial Intelligence, vol. 1, 883–888.

[William and Hoel(2003)]WilliamSeasonal03
William, B. and Hoel, L., 2003. Modeling and forecasting
vehicular traffic flow as a seasonal ARIMA process, Journal of
Transportation Engineering, vol. 129(6), 664–672.

[Wang and Xiao(2003)]WangRadialBasis03
Wang, X. and Xiao, J., 2003. A radial basis function neural
network approach to traffic flow forecasting, in: Proceedings of IEEE
Intelligent Transportation Systems, vol. 1, 614–617.

[Yu et al.(2003)]YuMarkov03
Yu, G., Hu, J., Zhang, C., Zhuang, L., and Song, J., 2003. Short-term
traffic flow forecasting based on markov chain model, in: Proceedings
of IEEE Intelligent Vehicles Symposium, 208–212.


[Zhang and Sun(2010)]ZhangMultiple10
Zhang, Q. and Sun, S., 2010. Multiple-view multiple-learner active learning, Pattern Recognition, vol. 43, 3113–3119.



