Predicting Station-level Hourly Demands in a Large-scale Bike-sharing Network: A Graph Convolutional Neural Network Approach

In plain words

Bike-sharing systems run out of bikes at some stations and fill up at others, so the operator needs to know how many people will want a bike at each station in the next hour. The usual approach predicts each station separately, which ignores that nearby or similarly used stations tend to move together.

This paper's model treats the network as a map: each station is a dot, and the model learns which dots should influence each other directly from the demand history, instead of relying on a link map a person drew up beforehand (for example, by distance alone). It adds a memory component, so the model can also use the recent sequence of demand, not only the last hour.

The authors report that this beats several simpler and more standard models on New York City data. The learned links also turn out to resemble, and go beyond, the hand-made links they tried.

TL;DR

• Station-level hourly bike demand is forecast with a graph convolutional network whose adjacency matrix (the graph filter) is a trainable parameter, not a predefined input §3.3.1.
• Two variants: a regular one (GCNN_reg-DDGF) and one with an LSTM block over time (GCNN_rec-DDGF) §3.3.2.
• Author-reported: on 272 Citi Bike stations, GCNN_rec-DDGF has the lowest test RMSE (2.12) and the highest R² (0.75) of 13 models tab. 3.
• Use it when stations are linked by a graph you cannot define well. Results are single-dataset and author-reported, and the public code does not reproduce the best (recurrent) variant.

Problem

Dock-based and free-floating bike systems have stations that run empty or full because demand fluctuates in space and time §1. Accurate station-level forecasts support rebalancing §1. Earlier station-level models were fitted one station at a time, and they did not use correlations between stations §2.1.3. Grid-based CNNs need a fixed grid, which is too coarse for station-level prediction or too sparse if the grid is fine §2.1.3 §1. Graph CNNs can handle the irregular station layout, but they need a predefined adjacency matrix, and choosing that matrix is hard §1 §3.2.

Assumptions: the demand history of each station is available hour by hour, and the network size is fixed (N = 272 stations in the study) §4. The paper does not study new stations or changing networks.

Method

Data. Each station has one demand value per hour, counted from bike check-outs at that station, over 2013-07 to 2016-06, which is 26,304 hours §4 tab. 2. Stations that were not present in all three years, or whose three-year demand is below 26,304 (under one bike per hour), are removed, leaving 272 stations §4. The first 22,304 hourly values are used for training, the next 2,000 for validation and the last 2,000 for testing §4.

Graph convolution. For an undirected graph with adjacency matrix A, the normalised Laplacian is L = I_N − D^(−1/2) A D^(−1/2) eq. 1. A first-order approximation of a spectral filter replaces the filter with a single term applied to à = A + I_N, with D̃ its degree matrix §3.1.1 eq. 6. A layer then maps the signal matrix X to D̃^(−1/2) à D̃^(−1/2) X Θ eq. 7. Stacking layers with ReLU gives the update in eq. 9: each hidden layer multiplies the previous activations by the normalised filter, then by a weight matrix, then applies ReLU §3.1.2 eq. 8 eq. 9. The output layer has the same form without ReLU eq. 10.

Four predefined graphs (baselines). Each one sets the entries of à to 1 or 0 by a threshold:

• Spatial distance (SD): connect stations within κ_SD eq. 11 §3.2.1.
• Demand (DE): aggregate origin–destination demand in both directions, connect pairs with total demand at least κ_DE eq. 12 §3.2.2.
• Average trip duration (ATD): divide total trip duration by DE, connect pairs whose average is at most κ_ATD eq. 13 §3.2.3.
• Demand correlation (DC): Pearson correlation of the two stations' hourly demand series (computed from training data only), connect pairs at or above κ_DC eq. 14 §3.2.4 §5.

Data-driven graph filter (the contribution). Treat the normalised filter  = D̃^(−1/2) à D̃^(−1/2) as an unknown symmetric N×N matrix of trainable parameters, and learn it together with the weights §3.3.1 eq. 15. The paper calls this the data-driven graph filter (DDGF). It needs no predefined edge set.

GCNN_reg-DDGF. A stack of convolution blocks (the DDGF applied to the signal) and feedforward blocks (a weight matrix plus ReLU). The paper's toy example uses three stations and shows the two steps per layer §3.3.2 fig. 1a.

GCNN_rec-DDGF. Per time step e, the station vector x^e (length N) is multiplied by Â, and the result feeds an LSTM with T cells and d hidden units per cell. The last hidden state h^t goes to a fully connected block that predicts the demand of all stations at t + 1 §3.3.2 fig. 1b. The LSTM gate equations are not written out; the paper refers to earlier work for them §3.3.2.

Targets and inputs. The model input is the last C0 hourly vectors of all stations, and the target is the next hour's vector for all stations §5. The paper uses C0 = 24 for the regular model tab. 2. Each training record starts at hour C0 §5.

Training and selection. Hyperparameters are chosen by grid search, judged on the validation set §5.1. Training uses a gradient method that the paper calls "classic Stochastic Gradient Descent" §5.1. Early stopping ends training when validation performance has not improved for s iterations, and the best model on validation is kept for testing §5.1.

# GCNN-DDGF forward pass  [§3.3.2] [eq. 9] [eq. 10] [eq. 15]
# Input:  signals (B, N, C0), N = 272 stations, C0 = 24 hours
# Output: next-hour demand (B, N)

def gcnn_ddgf_forward(signals, filters, weights, biases):
    # Start from the raw input window  [§5]
    hidden = signals                                       # (B, N, C0)
    for layer in hidden_layers:
        # Mix each station with its neighbours through the learned filter
        # [eq. 15] [§3.3.1]
        mixed = graph_filter(filters[layer], hidden)       # (B, N, C_prev)
        # Feature transform and ReLU  [eq. 9] [code: gcn.py#L35-L36]
        hidden = relu(mixed @ weights[layer] + biases[layer])
    # Output layer: the paper applies the filter here too [eq. 10]; the
    # code does not (see Paper vs code)  [code: gcn.py#L43]
    return hidden @ weights_out + bias_out                 # (B, N)

# Graph filter used at each layer  [eq. 15] [§3.3.1]
# Symmetrisation and absolute values are code-only, not in the paper
# [code: gcn.py#L31-L32] [code: utils.py#L7]

def graph_filter(raw, signal):
    # Make the learned matrix symmetric  [code: gcn.py#L31]
    symmetric = 0.5 * (raw + transpose(raw))              # (N, N)
    # Use magnitudes, so negative entries count as positive  [code: utils.py#L7]
    magnitude = abs(symmetric)
    # Symmetric degree normalisation, D^(-1/2) A D^(-1/2)  [eq. 15]
    degree = row_sums(magnitude)
    scale = diag(degree ** -0.5)
    normalised = scale @ magnitude @ scale
    return normalised @ signal

# Training with validation-based early stopping  [§5.1]
# Settings used by the code are listed under Paper vs code.

def train(train_set, val_set, test_set, patience):
    # Track the best validation RMSE seen so far  [eq. 16] [§5.1]
    best_val = infinity
    steps_without_gain = 0
    test_at_best = None
    for epoch in range(max_epochs):
        for batch in shuffled(train_set, batch_size):
            # Gradient step on the training RMSE  [eq. 16] [§5.1]
            update_parameters(batch)
        val_rmse = rmse(model(val_set.inputs), val_set.targets)
        if val_rmse < best_val:
            best_val = val_rmse
            # Report test error at the best validation epoch  [§5.1]
            test_at_best = rmse(model(test_set.inputs), test_set.targets)
            steps_without_gain = 0
        else:
            steps_without_gain += 1
        if steps_without_gain == patience:
            break
    return test_at_best

Metrics. Test RMSE is the root of the mean squared error over all hours and stations eq. 16. MAE is the mean absolute error eq. 17. R² compares the squared error with the variance of the test demand around its mean ȳ eq. 18. The peak-hour RMSE in tab. 3 uses only 7:00 AM to 9:00 PM, since demand is near zero at other times §5.2.

Hyperparameters

Table 2 gives the selected values; Table 1 gives the search ranges. The notation {start:step:end} in Table 1 includes the end value tab. 1 §5.1.

Hardware and software: a GTX 1080 GPU, 64 GB RAM, Ubuntu 16.04 and TensorFlow with Python §5.1. Training time is not reported §5.

Results

Author-reported, on the 2,000-hour test set after selection on validation. Setup: 272 Citi Bike stations in New York City, 2013–2016 §4. Baselines are the four predefined-graph GCNNs, a per-station MLP, LSTM, XGBoost, LASSO, linear and RBF SVR, and the historical average (HA) §5.2.

Validation RMSE for the six GCNN models tab. 2:

Test set, all 13 models tab. 3:

The paper reports no confidence intervals and no variance across runs §5, so small RMSE gaps between neighbouring rows cannot be judged from the paper alone.

Learned-graph analysis §5.3, author-reported and descriptive: the learned graph has 1,565 edges with weight at least 0.15 after normalising to [0, 1] §5.3; eight communities cover 251 of 272 stations §5.3; the weighted degree correlates with neighbour count with R² = 0.96 fig. 6 §5.3.1. The authors read the learned edges as partly matching the SD, DE and DC matrices and partly not §5.3.2 fig. 8 fig. 9.

Paper vs code

The official repository is leilin-research/GCGRNN at commit c44f9e255bcd65ebbcdc6e63beeb12d79161402a README.md. Its GCNN-DDGF_bike_sharing/ directory holds the regular variant only README.md. It has no LSTM or recurrent block for the NYC data (the word "recurrent" appears only in the name of the repository's other model, GCGRNN) README.md. The code therefore cannot be checked against GCNN_rec-DDGF.

The repository has no LICENSE file, so its code is not licensed for reuse. Nothing here is quoted from it; the citations are for checking facts only.

Discrepancies for the regular variant:

• Architecture. Paper: one hidden layer of width 40, C2 = 0 tab. 2. Code: three hidden layers of widths 10, 10 and 20 GCNN_bike_sharing.py. The width 10 lies outside the paper's search sets for C1 and C2 tab. 1.
• Output layer graph filter. Paper eq. 10 applies the filter in the output layer eq. 10. Code's output is a plain linear map of the last hidden layer, with no adjacency gcn.py.
• Filter per layer. Eq. 15 writes one Â, with no layer index eq. 15. Code learns a separate N×N matrix per hidden layer gcn.py.
• Symmetrisation and absolute values. Not in the paper. Code symmetrises the learned matrix and takes absolute values before normalising gcn.py utils.py.
• Biases. Eq. 9 has no bias. Code adds one per layer gcn.py.
• Optimiser. Paper: classic SGD §5.1. Code: Adam gcn.py. An RMSProp line with the learning-rate decay is commented out gcn.py. The decay of 0.9 is defined but not used GCNN_bike_sharing.py.
• Learning rate. Paper: 0.005 tab. 2. Code: 0.01 GCNN_bike_sharing.py.
• Early-stopping patience. Paper: 20 tab. 2. Code: 200, although its comment says 20 GCNN_bike_sharing.py. The code counts patience in epochs gcn.py.
• Training size. Paper: 22,304 hours §4. Code: 20,000 records for training, with 2,000 validation and 2,000 test GCNN_bike_sharing.py.
• Shuffling. Code shuffles the training records each run, unseeded GCNN_bike_sharing.py. Paper: not stated.
• Input scaling. Paper: min–max to [0, 1] §5. Code: standardisation with the mean and standard deviation of the whole training input GCNN_bike_sharing.py utils.py.
• Dropout. Code sets keep probability 1 GCNN_bike_sharing.py, and the forward function takes no dropout argument gcn.py. Paper: not stated.
• Epochs. Code: up to 500 GCNN_bike_sharing.py. Paper: not stated.
• Selection. Code records test RMSE at the epoch of best validation RMSE gcn.py. This matches the paper's description §5.1.
• Metrics. Code computes RMSE only. MAE and R² in tab. 3 are not computed in the repository gcn.py.
• Framework. Code uses TensorFlow 1.x graph and session APIs gcn.py gcn.py. The paper says TensorFlow without a version §5.1.
• Data. The code loads a pre-built 272 × 26,304 pickle GCNN_bike_sharing.py. The paper's aggregation steps are not in the repository, so they cannot be checked against code.

Pitfalls

• Only the regular variant is implemented publicly. The recurrent variant, which the paper reports as best, has no code. Implement it from the text and Table 2 only; the LSTM gate equations are not written in the paper §3.3.2.
• Parameter count. Each learned filter has N² = 272 × 272 entries, about 74 thousand values, before the weights §4 eq. 15. The paper does not say how symmetry is enforced; the code uses 0.5(A + Aᵀ) gcn.py.
• Zero-degree rows. The normalisation divides by the square root of the row sum. If a learned row is all zeros, its degree is zero and the scaling is undefined utils.py. The paper does not discuss this.
• Self-loops. Eq. 6 adds the identity to the adjacency before normalising eq. 6. For the learned filter, the diagonal is just another trainable value; the paper does not say whether the identity is added eq. 15.
• Unstated initialisation. The paper gives no initial values for the filter or weights. The code draws them from a normal distribution with standard deviation 0.5 gcn.py.
• Run-to-run variance. Neither the shuffle nor the initialisation is seeded in the code GCNN_bike_sharing.py gcn.py. The paper does not say whether it repeated runs §5, so the gap between the two DDGF variants (2.12 vs 2.35) has no reported uncertainty.
• Selection on the test set. The code computes the test error every epoch and keeps the one at the best validation epoch gcn.py. It does not use the test error for selection, but the test set is evaluated repeatedly, so it should not be used for further tuning.
• Recurrent model input. For the LSTM variant, C0 is "N/A" in Table 2; the sequence length T = 24 plays that role, and each time step feeds one hourly vector tab. 2 §3.3.2.
• Interpretation of learned edges. Weighted degree and edge weights come from a single trained model. The paper's claim that the learned graph finds links not visible in SD, DE or DC is based on rank comparisons for four stations §5.3.2 fig. 9.
• Demand is mostly near zero outside peak hours. The paper reports the 7:00 AM–9:00 PM RMSE separately for that reason §5.2. Busy stations also have much larger error, which the paper links to their weighted degree §5.3.1 fig. 7.
• Edge direction. The learned graph is treated as undirected, and the authors list directed graphs as future work §6.

Lineage

• Builds on: spectral graph convolution, as in Defferrard et al. (2016) and Kipf and Welling (2016), and the first-order filter used by the latter §2.2 §3.1.1. Neither has a card here.
• Builds on: the LSTM §3.3.2 (Hochreiter and Schmidhuber, 1997). No card here.
• Compared against: four GCNNs with predefined graphs (SD, DE, ATD, DC), an LSTM, XGBoost, MLP, LASSO, SVR and the historical average §5.2 tab. 3.
• Related code in the same repository: GCGRNN, which the README describes as "improved from" GCNN-DDGF, for traffic speed and volume data README.md. Not compared in the paper.
• Superseded by: not known from content in this wiki.