2.1 Preliminary data analysis

When examining the bicycle-sharing data sets, we observe strong daily usage patterns, which motivate the development of our time-dependent SBMs. In Figure 1, we illustrate these usage patterns by plotting the number of trips by starting hour for each city. In each city, weekday activity spikes during morning and evening commuting hours (and Los Angeles also has a mid-day spike), whereas weekend trips peak in early afternoon. Similar patterns have been observed previously for bicycle-sharing systems in New York City and many other cities (Austwick et al. Reference Austwick, O’Brien, Strano and Viana2013; Etienne and Latifa Reference Etienne and Latifa2014; Taiyeb Reference Taiyeb2014; Xie and Wang Reference Xie and Wang2018; Zhu et al. Reference Zhu, Khan, Kats, Santosh Bamne and Sobolevsky2018).

For individual stations, the morning and evening peaks for arriving trips (i.e., in-edges) and departing trips (i.e., out-edges) are often unbalanced. One direction has a stronger morning peak, and the opposite direction has a stronger evening peak. In the bottom panel of Figure 1, we show examples of this imbalance for the six busiest stations in San Francisco. (In Appendix A.1, we further explore this issue using singular-value decompositions of the trip data.) These patterns suggest how stations are used in commuting, and they motivate our time-dependent identification of stations into “home” and “work” types.

Another characteristic of our data that we incorporate into the design of our models is the strong positive Pearson correlation coefficient between the total (summed over all hours) in-degree and out-degree of each station; it is $0.99$ in New York City, $0.98$ in San Francisco, and $0.91$ in Los Angeles. Accordingly, we apply degree-correction to the overall degree of each station, rather than to the in-degrees and out-degrees separately (see Section 3.1). This correlation is an intrinsic feature of docked bicycle-sharing systems; at some point, bicycles must be returned to a station before any new trips can originate from it. (The use of trucks to redistribute bicycles within a system can loosen this requirement.) More generally, the strong correlation also relates to a famous axiom of human mobility: for each current of travel, there is an associated countercurrent (Ravenstein Reference Ravenstein1885).

3.1 Time-dependent mixed-membership stochastic block model (TDMM-SBM)

We now describe the framework for our mixed-membership SBM. Let $i,j \in \mathcal{N}$ (with $|\mathcal{N}| = N$) be nodes, which represent bicycle-sharing stations; let $g,h \in \mathcal{K}$ (with $|\mathcal{K}| = K$) be blocks. We treat both the number of nodes and the number of blocks as fixed and given. To each bicycle-sharing network, we associate a three-dimensional array of size (N,N,T), where T is the number of time slices (i.e., time layers). Our model is also valid when layers have meanings other than time slices. However, because of our application to bicycle-sharing systems, we describe the layers as time layers for convenience. We consider hourly groupings of bicycle trips based on their starting times, and we do not include interlayer edges.^{Footnote 2} The array entry $A_{ijt}$ is the observed number of trips from station i to station j with a starting time that is at least t and less than $t+1$. Our networks are directed, so we count each trip that starts at node i in hour t and ends at node i (i.e., self-edges) exactly once in $A_{iit}$. We use $\tilde{A}_{ij}=\sum_{t=0}^{23}A_{ijt}$ to denote the edge weights of the associated time-aggregated matrix.

For each node i, there is a latent length-K vector of real numbers $C_{ig}\in [0,1]$. These vectors, which we will estimate, signify the mixed-membership block assignments of the nodes. The block-assignment parameter $C_{ig}$ indicates the “strength” of node i in block g. For each ordered pair g,h of blocks and each time $t\in \{0, 1, \ldots, 23\}$ (where $t = 0$ represents the hour that starts at midnight), there is a parameter $\omega_{ght}$, which we call the “block-connectivity” parameter (i.e., the “block-to-block” activity parameter), that represents the directed activity from block g to block h in layer t. Note that $\omega_{ght}$ need not be equal to $\omega_{hgt}$ for a directed network; in our application, this captures asymmetries in the numbers of trips with respect to reversing origins and destinations. We also define the notation $\tilde{\omega}_{gh}=\sum_{t=0}^{23} \omega_{ght}$ for our time-aggregated matrix.

For each pair of nodes, i and j, we assume that the number of trips that depart from i in hour t and arrive at j at any time is Poisson-distributed with mean $\mu_{ijt} = \sum_{g,h} C_{ig}\omega_{ght} C_{jh}$. Our use of the Poisson distribution follows Karrer and Newman (Reference Karrer and Newman2011) and Peixoto (Reference Peixoto2019), facilitates computation, and is standard for modeling count data (although overdispersion is a concern).

For identifiability (i.e., to allow unambiguous inference), we apply the constraint $\sum_i C_{ig}=1$ for all g. This does not constrain the set of possible models in terms of realizable mean edge activities $\mu_{ijt}$. Consider a model with unconstrained parameters $\omega_{ght}$ and $C_{ig}$. A model with parameters $\omega'_{ght}$ and $C'_{ig}$ such that $C'_{ig}=\frac{C_{ig}}{\sum_{j} C_{jg}}$ and $\omega'_{ght}=\omega_{ght}\left(\sum_{j} C_{jg}\right)\left(\sum_{j}C_{jh}\right)$ gives an equivalent parameterization, because the means of the distributions of edge weights are equal to the those of the model with unconstrained parameters. That is, $\mu'_{ijt}=\sum_{g,h}C'_{ig}\omega'_{ght}C'_{jh}=\sum_{g,h}C_{ig}\omega_{ght}C_{jh}=\mu_{ijt}$.

Given that $\sum_{i} C_{ig}=1$, we can think of $C_{ig}$ as the proportion of the total activity of block g from the activity of node i. We interpret $C_{ig}$ relative to $C_{ih}$ as how strongly block g is associated with node i relative to how strongly block h is associated with node i. The expected total number of trips at node i is $\sum_g C_{ig}\sum_{h,t} (\omega_{ght}+\omega_{hgt})$. By contrast, $\sum_g C_{ig}$ is a measure of the activity of node i in which we do not weight $C_{ig}$ by the total activity of its corresponding block. We use $\sum_g C_{ig}$ to determine the size of node i when plotting the block assignments that we infer using the TDMM-SBM (see the left panels of Figures 3 and 6) because it helps ensure that we do not overlook blocks with important usage patterns but relatively low activity. The parameters $C_{ig}$, with i fixed and $g \in \mathcal{K}$, extend the degree-correction parameter for SBMs that was introduced by Karrer and Newman (Reference Karrer and Newman2011) to mixed block membership. We elaborate on this connection in Section 3.2, where we introduce a model that specifies that nodes are assigned to only one block.

Figure 2. In the TDD-SBM networks that we construct, the red dots indicate the mean weights of the edges from block g to block h in time layer t, where $g,h \in \{1,2,3\}$ and $t \in \{0, 1, \ldots, 15\}$. The title of each panel indicates the blocks g and h (e.g., “1 to 2”).

Figure 3. Classification of bicycle-sharing stations in downtown Los Angeles using (left) a two-block TDMM-SBM and (right) a two-block TDD-SBM. The sizes of the nodes take continuous values. In the left panel, we scale the size of each node based on its value of $\sum_gC_{ig}$. In the right panel, we scale the size of each node based on its degree.

We now compute the likelihood function that we will optimize to obtain maximum-likelihood estimates (MLEs) of parameters. We assume conditional independence between the Poisson-distributed hourly numbers of trips along each edge, given model parameters, so the likelihood of the data is

(2) \begin{align} L(G;\,\mathbf{\omega},\mathbf{C})=\prod_{t=0}^{23}\prod_{i,j\in\mathcal{N}} \frac{(\mu_{ijt})^{A_{ijt}}}{A_{ijt}!}\exp\left(-\mu_{ijt}\right)\,,\end{align}

where $\mathbf{\omega}$ and $\mathbf{C}$ are the model parameters (i.e., $\mathbf{\omega}=\{\omega _{ght}\}$ and $\mathbf{C}=\{C_{ig}\}$). Note that $\mu_{ijt}=\sum_{g,h} C_{ig}\omega_{ght} C_{jh}$ is a function of these parameters, the set $\mathcal{N}$ of nodes of the network is fixed and given, and the number K of blocks is also fixed and given. The unnormalized log-likelihood is

(3) \begin{align} \ell(G;\,\mathbf{\omega},\mathbf{C})=\sum_{t=0}^{23}\sum_{i,j\in\mathcal{N}}\left[A_{ijt}\log\left({\mu_{ijt}}\right)- \mu_{ijt}\right]\,,\end{align}

where we omit the constant ${-}\!\sum_{i,j,t}\log\left(A_{ijt}!\right)$ because it does not affect the maximum of the function. In (3) and throughout this paper, we use “log” to denote the natural logarithm.

3.2 Time-dependent discrete stochastic block model (TDD-SBM)

We derive a single-membership SBM from our mixed-membership SBM by making the extra assumption, for each node $i\in \mathcal{N}$, that $C_{ig}>0$ for only one block $g \in \mathcal{K}$. For our single-membership SBM, we introduce some new notation to aid our description and promote consistency with the notation in other work (Karrer and Newman Reference Karrer and Newman2011; Zhu et al. Reference Zhu, Yan and Moore2013). For a given node i, the block g for which $C_{ig}>0$ is the block $g_i$ that is associated with node i. Therefore, we describe this model as a discrete-membership model. We introduce a scalar parameter $\theta_i=C_{ig_i}$ for each node i; it is a multilayer extension of the degree-correction parameter in Karrer and Newman (Reference Karrer and Newman2011). The expectation of the Poisson distribution of the value of an edge from node i to node j at time t is $\theta_i\theta_j\omega_{g_ig_jt}=C_{ig_i}\omega_{g_{i},g_{j}t} C_{jg_{j}}{=\mu_{ijt}}$. We retain the sum constraints of our mixed-membership model, so $\sum_{i\in g}\theta_i = 1$ for all g.

We compute optimal values of the parameters $\mathbf{\omega}$ and $\mathbf{\theta}=\{\theta_i\}_{i\in\mathcal{N}}$. As in our TDMM-SBM, we take $\mathcal{N}$ and K to be fixed and given. Again dropping the constant term $-\!\sum_{i,j,t}\log\left(A_{ijt}!\right)$, the log-likelihood of our discrete-membership SBM is

(4) \begin{align} \ell(G;\,\mathbf{\omega}, \mathbf{\theta}) = \sum_{t}\sum_{g,h}\sum_{i\in g, j\in h} \left[A_{ijt}\log\theta_i+A_{jit}\log \theta_j+ A_{ijt}\log \omega_{ght} -\theta_i\theta_j \omega_{ght}\right].\end{align}

This resembles the degree-corrected SBM in equation (14) of Karrer and Newman (Reference Karrer and Newman2011), but we have adapted it to directed networks and summed over time layers.

We find explicit formulas for the conditional MLEs of $\theta_i$ and $\omega_{ght}$, given the block memberships $\{g_i\}_{i \in \mathcal{N}}$. In the following calculations, the absence of t in the subscript of a parameter and the presence of a tilde designates a sum over all t. Specifically, we define $\tilde{A}_{ij}=\sum_{t=0}^{23} A_{ijt}$ and $\tilde{\omega}_{gh}=\sum_{t=0}^{23} \omega_{ght}$. We differentiate $\ell$ with respect to $\omega_{ght}$ to obtain

\begin{align*} \frac{\partial}{\partial \omega_{ght}} \ell=\frac{\sum_{i \in g, j \in h}A_{ijt}}{\omega_{ght}}-1\,,\end{align*}

where we have used the blockwise sum constraints on $\theta$ to simplify $\sum_{i\in g, j\in h}\theta_i\theta_j$ to 1. Therefore, the conditional MLE for $\omega_{ght}$ is

(5) \begin{align} \hat{\omega}_{ght}=m_{ght}\,,\end{align}

where $m_{ght}$ is the sum of the weights of the edges from nodes in block g to nodes in block h in layer t. That is, $m_{ght} = \sum_{i\in g,j\in h}A_{ijt}$.

We then differentiate $\ell$ with respect to $\theta_i$ to obtain

\begin{align*} \frac{\partial}{\partial\theta_i}\ell=\frac{\sum_j \tilde{A}_{ij}+\sum_{j} \tilde{A}_{ji}}{\theta_i}-\sum_{h}\tilde{\omega}_{g_i h}-\sum_{h}\tilde{\omega}_{h g_i}\,.\end{align*}

At $\hat{\omega}_{ght}$, the conditional MLE for $\theta_i$ is

(6) \begin{align} \hat{\theta_i}=\frac{\sum_j \big(\tilde{A}_{ij}+\tilde{A}_{ji}\big)}{\sum_{g}\big(\tilde{m}_{gg_i}+\tilde{m}_{g_ig}\big)} =\frac{k_i}{\kappa_{g_i}}\,,\end{align}

where $k_i=\sum_j \left(\tilde{A}_{ij}+\tilde{A}_{ji}\right)$ is the degree of node i (i.e., the sum of the in-degree and the out-degree^{Footnote 3} of i) over all time layers, $\tilde{m}_{gh}=\sum_{t=0}^{23} m_{ght}$, and $\kappa_{g}=\sum_{h} \left(\tilde{m}_{gh}+\tilde{m}_{hg}\right)$ is the sum of the in-degrees and out-degrees of all nodes in block g over all time layers. The term $2\tilde{m}_{gg}$ in the equation for $\kappa_g$ implies that we count each intra-block edge twice: once for departing from g and once for arriving at g. Similarly, $k_i$ includes the term $2\tilde{A}_{ii}$, so we count self-edges twice in this term.

The parameter MLE $\hat{\theta}_i$ in our TDD-SBM for directed, multilayer networks is analogous to the MLE of a degree-correction parameter in a degree-corrected SBM for undirected networks with one layer (Karrer and Newman Reference Karrer and Newman2011). Indeed, the MLE for the degree-correction parameter of a node in the latter model is the degree of that node divided by the sum of the degrees of all nodes in its block. Another similarity between the degree-corrected SBM of Karrer and Newman (Reference Karrer and Newman2011) and our TDD-SBM is that at the MLE of the TDD-SBM, the expectation of the (time-aggregated) degree of a node i is equal to the degree of node i from the observed data. That is, $\sum_{t}\sum_{j}\left(\mu_{ijt}+\mu_{jit}\right) = k_i$. (See Appendix A.2 for the proof.) For our mixed-membership SBM, we are not aware of such a precise relationship between the data and the expected value of any model statistics.

We now substitute the above conditional parameter estimates into the unnormalized log-likelihood of the TDD-SBM to obtain

\begin{align*} \sum_{t}\left[\sum_{i,j} \left(A_{ijt}\log \left(\frac{k_i}{\kappa_{g_i}}\right)+A_{jit}\log\left(\frac{k_j}{\kappa_{g_j}}\right)\right)+\sum_{g,h} m_{ght}\log m_{ght} - \sum_{g,h} m_{ght}\right] \\[3pt] =\sum_i k_i \log k_i - \sum_i k_i\log {\kappa_{g_i}}+\sum_t\sum_{g,h} m_{ght}\log m_{ght} -\tilde{m}\,,\end{align*}

where $\tilde{m}$ is the total number of edges in the network. We do a similar calculation as one in Karrer and Newman (Reference Karrer and Newman2011) and write

\begin{align*} \sum_{i} k_i\log \kappa_{g_i}&=\sum_t\sum_{g} \sum_{i\in g} k_{it} \log\kappa_g\\&=\sum_t \sum_{g}\sum_{i\in g}\left( k_{\text{in},it} \log\kappa_g + k_{\text{out},it} \log\kappa_g\right)\\&=\sum_t\sum_g\kappa_{\text{out},gt}\log\kappa_g+\sum_t\sum_h\kappa_{\text{in},ht}\log \kappa_h\\&=\sum_t\sum_g\sum_h m_{ght}\log\kappa_g+\sum_t\sum_g\sum_h m_{ght}\log \kappa_h\\&=\sum_{t}\sum_{g,h} m_{ght}\log \kappa_g\kappa_h\,,\end{align*}

where $k_{\text{in},it}$ and $k_{\text{out},it}$ are the in-degrees and out-degrees of node i in layer t, the quantity $\kappa_{\text{in},gt}=\sum_{i\in g} k_{\text{in},it}$ is the number of edges that arrive at g in hour t, and $\kappa_{\text{out},gt}=\sum_{i\in g} k_{\text{out},it}$ is the number of edges that depart from g in hour t. Including only the terms that depend on block assignments yields the objective function

(7) \begin{align} \sum_t\sum_{g,h} m_{ght}\log\left(\frac{m_{ght}}{\kappa_g\kappa_h}\right)\,,\end{align}

which we maximize to obtain the MLE of the block assignments.

Unlike the directed SBM of Zhu et al. (Reference Zhu, Yan and Moore2013), our TDD-SBM does not have two degree-correction parameters for each station. Their directed SBM includes a parameter $\theta_i^{\text{in}}$ to modify in-degree and a parameter $\theta_i^{\text{out}}$ to modify out-degree. The validity of using undirected degree-correction parameters $\{\theta_i\}_{i\in\mathcal{N}}$ to model directed networks depends on the presence of a large correlation between the time-aggregated in-degrees and time-aggregated out-degrees of the nodes. In this situation, including distinct $\theta_i^{\text{in}}$ and $\theta_i^{\text{out}}$ parameters would lead to overfitting. The Pearson correlations between in-degree and out-degree are larger than $0.9$ for all of our bicycle-sharing networks (see Section 2). To increase the applicability of our work to other networks (which may not have such large correlations), we provide estimation code in our implementation materials (see Section 4.2) for a version of our TDD-SBM with directed degree-correction parameters $\theta_i^{\text{in}}$ and $\theta_i^{\text{out}}$.

4.1 Inference using the TDMM-SBM

Let $\mathbf{\omega}=\{\omega_{ght}\}$ be the $K\times K\times T$ array of block-connectivity parameters, and let $\mathbf{C}=\{C_{ig}\}$ be the matrix of node-strength parameters. We estimate the model parameters using a two-step gradient descent.^{Footnote 4} First, we move in the direction of the gradient with respect to $\mathbf{\omega}$ and update the block-connectivity parameters. Second, we move along the direction of the gradient with respect to $\mathbf{C}$ and update the node-strength parameters.

In our algorithm, $\omega^{(n)}$ and $\mathbf{C}^{(n)}$ denote the $n^{\text{th}}$ updates of the block-connectivity and node-strength parameters, respectively. We initialize the algorithm with random values $\mathbf{\omega}^{(0)}$ and $\mathbf{C}^{(0)}$ with entries that are distributed according to $\exp(X)$, where X is a Gaussian random variable with mean 0 and variance 1. (That is, we draw random values from a log-normal distribution.) We denote the mean activity along edge (i,j) with initial parameters $\mathbf{\omega}^{(0)}$ and $\mathbf{C}^{(0)}$ by $\mu_{ijt}^{(0)}$. We scale the initial parameters so that the TDMM-SBM at the starting point of the optimization procedure has the same mean total number of trips as the data. Specifically, we multiply the block-connectivity parameters $\omega^{(0)}_{ght}$ by $\left(\sum_{t}\omega^{(0)}_{ght}\right)^{-1}\big(\sum_{i,j,t} A_{ijt}\big)\big/K^2$ and normalize $\mathbf{C}^{0}$ to satisfy the constraint $\sum_i C_{ig}^{(0)}=1$ for each block g. This yields $\sum_{ijt} \mu_{ijt}^{(0)}=\sum_{i,j,t} \sum_{g,h} C_{ig}^{(0)}\omega_{ght}^{(0)} C_{jh}^{(0)}=\sum_{g,h,t} \omega_{ght}^{(0)}=\sum_{g,h}\left(\sum_{i,j,t} A_{ijt}\right)\big/K^2=\sum_{i,j,t} A_{ijt}$, which is the total number of edges in the network. Without this scaling, the initial parameters can have very small magnitudes, such that the mean total number of trips from the TDMM-SBM with these initial parameters is much smaller than the total number of trips in the data. In that situation, the initial likelihood is much smaller than the likelihood at the MLE.

To ensure that our estimated parameters are nonnegative, we use the following change of variables: $\exp(\tilde{\omega}^{(n)})=\omega^{(n)}$ and $\exp(\tilde{\mathbf{C}}^{(n)})=\mathbf{C}^{(n)}$, where the exponential of a matrix or a vector signifies componentwise exponentiation of each of its entries. We can then write the gradient descent as

\begin{align*} \tilde{\mathbf{\omega}}^{(n+1)} &= \tilde{\mathbf{\omega}}^{(n)}+\eta^{(n)}\nabla_{\mathbf{\omega}} \ell\left(\mathbf{C}^{(n)},\mathbf{\omega}^{(n)}\right)\exp\left(\tilde{\mathbf{\omega}}^{(n)}\right)\,, \\ \tilde{\mathbf{C}}^{(n+1)} &= \tilde{\mathbf{C}}^{(n)}+h^{(n)}\nabla_{\mathbf{C}}\ell\left(\mathbf{C}^{(n)},\mathbf{\omega}^{(n+1)}\right)\exp\left(\tilde{\mathbf{C}}^{(n)}\right)\,,\end{align*}

where $h^{(n)}$ and $\eta^{(n)}$ are small positive numbers. From the definitions of $\tilde{C}^{(n)}$ and $\tilde{\omega}^{(n)}$, we write

\begin{align*} \mathbf{\omega}^{(n+1)} &= \mathbf{\omega}^{(n)}\exp\left(\eta^{(n)}\nabla_{\mathbf{\omega}} \ell\left(\mathbf{C}^{(n)},\mathbf{\omega}^{(n)}\right) \mathbf{\omega}^{(n)}\right)\,,\\\mathbf{C}^{(n+1)} &= \mathbf{C}^{(n)}\exp\left(h^{(n)}\nabla_{\mathbf{C}} \ell(\mathbf{C}^{(n)},\mathbf{\omega}^{(n+1)}) \mathbf{C}^{(n)}\right)\,.\end{align*}

Let $h^{(0)}=\eta^{(0)}=\Delta>0$ be the fixed initial step size. For our application, we choose $\Delta =10^{-4}$. We generate two candidate updates for $\mathbf{\omega}^{(n+1)}$ for the first step in our algorithm using $h^{(n+1)}=1.2\,h^{(n)}$ and $h^{(n+1)}=0.8\,h^{(n)}$, and we choose the one that gives the $\mathbf{\omega}^{(n+1)}$ value that yields the larger $\ell(\mathbf{C}^{(n)},\mathbf{\omega}^{(n+1)})$. Similarly, we choose the one of $\eta^{(n+1)}=1.2\,\eta^{(n)}$ and $\eta^{(n+1)}=0.8\,\eta^{(n)}$ that gives the $\mathbf{C}^{(n+1)}$ value that yields the larger $\ell(\mathbf{C}^{(n+1)},\mathbf{\omega}^{(n+1)})$.

We compute the gradient of the log-likelihood function $\ell$ using the chain rule. Recall that we compute the log-likelihood in two parts. One part is the computation of the mean $\mu_{ijt}=\sum_{g,h} C_{ig}\omega_{gh}^t C_{jh}$ of the number of trips from node i to node j at time t. We then insert the expression for the mean into $\ell=\sum_t\sum_{i,j}\left(A_{ijt}\log(\mu_{ijt})-\mu_{ijt}\right)$. We compute the derivative of $\ell$ with respect to $\mu_{ijt}$ to obtain

\begin{align*} \frac{\partial \ell}{\partial \mu_{ijt}} = \frac{A_{ijt}}{\mu_{ijt}} - 1\,.\end{align*}

The derivative of $\mu_{ijt}$ with respect to $C_{kg}$ is

\begin{align*} \frac{\partial \mu_{ijt}}{\partial C_{kg}} =& \delta_{ki}\sum_{h} \omega_{ght} C_{jh} + \delta_{kj} \sum_{h}C_{ih} \omega_{hgt} \,,\end{align*}

where $\delta_{ab}$ denotes the Kronecker delta (i.e., $\delta_{ab}=1$ if $a=b$ and $\delta_{ab}=0$ if $a \neq b$). The derivative of $\mu_{ijt}$ with respect to $\omega_{ght}$ is

\begin{align*} \frac{\partial \mu_{ijt}}{\partial\omega_{ght}}=C_{ig}C_{jh}\,.\end{align*}

Using the above calculations, we see that the derivatives of $\ell$ with respect to $C_{kg}$ and $\omega_{ght}$ are

\begin{align*} \frac{\partial \ell} {\partial C_{kg}}&=\sum_{t=0}^{23}\sum_{i,j\in\mathcal{N}}\frac{\partial \ell}{\partial \mu_{ijt}}\frac{\partial \mu_{ijt}}{\partial C_{kg}}\\[3pt]&=\sum_{t=0}^{23}\left(\sum_{j\in\mathcal{N}}\left(\frac{A_{kjt}}{\mu_{kjt}}-1\right) \sum_{h}\omega_{ght}C_{jh}+\sum_{i\in\mathcal{N}}\left(\frac{A_{ikt}}{\mu_{ikt}} -1\right) \sum_{h}C_{ih}\omega_{hgt}\right) \,,\\[3pt] \frac{\partial\ell}{\partial \omega_{ght}} &= \sum_{i,j\in\mathcal{N}}\frac{\partial \ell}{\partial \mu_{ijt}}\frac{\partial \mu_{ijt}}{\partial \omega_{ght}} \\[3pt] &= \sum_{i,j\in\mathcal{N}} \left(\frac{A_{ijt}}{ \mu_{ijt}}-1\right) C_{ig} C_{jh}\,.\end{align*}

We run gradient descent until four significant digits of the base-10 floating-point representation of the log-likelihood (3) do not change for 600 consecutive steps. For the networks that we examine, this usually takes between 600 and 5,000 iterations, with models with more blocks generally needing more iterations to reach this stopping criterion. Because of the nonconvexity of the log-likelihood function (3), we are not guaranteed to reach a global optimum. Most of the time, our method converges to an interesting local optimum (which may also be a global optimum) that appears to reveal functional roles of bicycle-sharing stations (see Section 6). Our results produce recognizable block-connectivity parameters $\omega_{ght}$ (e.g., reflecting home–work commute patterns and leisure patterns) and the parameters $C_{ig}$ indicate known spatial divisions of the stations (e.g., residential versus commercial districts). In some cases, however, our algorithm converges to an uninteresting local optimum; for example, sometimes the block-assignment parameters $C_{ig}$ for each station appear as if they were assigned independently at random. To improve our results, we run our algorithm several times (specifically, 10 times for each network) and store the estimate with the largest likelihood.

For several examples, we compare the parameters that we obtain from gradient descent to those that we obtain by running a Hamiltonian Monte Carlo (HMC) sampling method in a Bayesian framework with weak priors. (We implement this sampling method in $\textsf{Stan}$ (Carpenter et al. Reference Carpenter, Gelman, Hoffman, Lee, Goodrich, Betancourt, Brubaker, Guo, Li and Riddell2017; Stan Development Team 2018).) The log-likelihoods that result from our gradient-descent method are as good or better than those that we obtain with the HMC method. The HMC method is more computationally intensive and memory intensive than our gradient-descent method, although it may be preferable in applications in which one has meaningful prior information about parameters. Improving our optimization method and investigating trade-offs between accuracy and efficiency are worthwhile topics for future work. For instance, it likely will be beneficial to adapt optimization methods for related time-dependent SBMs (Xing et al. Reference Xing, Fu and Song2010; Ho et al. Reference Ho, Song and Xing2011; Yang et al. Reference Yang, Chi, Zhu, Gong and Jin2011; Xu and Hero Reference Xu and Hero2014; Matias et al. Reference Matias, Rebafka and Villers2018) to the optimization of our SBMs.

The repository at https://github.com/jcarlen/tdsbm_supplementary_material has our $\textsf{Python}$ implementation of our gradient-descent method and code for our inference using <monospace>Stan</monospace> in $\textsf{R}$.

4.2 Inference using the TDD-SBM

To fit our TDD-SBM, we use a Kernighan–Lin-type (KL-type) algorithm (Kernighan and Lin Reference Kernighan and Lin1970) that we base on the one in Karrer and Newman (Reference Karrer and Newman2011). (Karrer and Newman (Reference Karrer and Newman2011) noted for their time-independent degree-corrected SBM that their KL-type algorithm is faster than node-switching algorithms and that it achieves results that are similar to those from such algorithms.) Given K blocks, we initialize the KL-type algorithm by assigning each node to a block uniformly at random, so each node has probability $1/K$ of being assigned to a given block. The algorithm then calculates the best possible block reassignment for any single node with respect to the associated change in log-likelihood (either the largest increase or the smallest decrease). It then chooses a new node uniformly at random and determines its best reassignment with respect to change in log-likelihood. The algorithm cycles through all nodes; a key feature of the algorithm is that a node that has been reassigned cannot change blocks again until all other nodes have been reassigned. A sequential reassignment of all nodes constitutes one step of the algorithm. The algorithm then searches all of the states that occurred during the step, and it selects the state with the maximum log-likelihood of any state in that step. This state is the starting point of the next step of the algorithm. A single run of the algorithm finishes when a step does not increase the log-likelihood beyond a preset tolerance value near 0. (In practice, we use $1 \times 10^{-4}$.) To find block assignments that are as good as possible, we do many runs of the algorithm for each network; this process is easily parallelizable. In our examples, we use 50 runs per network. We initialize each run randomly as described above.

Another key feature of the KL-type algorithm is that changes in the block memberships of nodes affect only the terms of the objective function that involve the origin and destination blocks of the changes. The simplified objective function (7) is a sum over block-pair terms over T time layers. Consequently, we do not need to recalculate the full objective function for each node reassignment. This estimation method is practical for networks with up to a few thousand nodes and a small number of blocks, depending on available computing resources.

We implement our KL-type algorithm for our TDD-SBM in an $\textsf{R}$ package using Rcpp (Eddelbuettel and FranÇois Reference Eddelbuettel and FranÇois2011; R Core Team 2018). The back-end calculations are in $\textsf{C++}$ for speed, and we return results in $\textsf{R}$ to enable visualization and other downstream tasks. One can also use our $\textsf{R}$ package to estimate the parameters of time-independent SBMs (including directed and/or degree-corrected ones), a variant of our TDD-SBM without degree-correction, and a variant of our TDD-SBM with directed degree-correction parameters $\theta_i^{\text{in}}$ and $\theta_i^{\text{out}}$. This facilitates comparisons of inferences from various related SBMs. See https://github.com/jcarlen/sbm for our $\textsf{R}$ package $\textbf{\texttt{sbmt}}$ for parameter estimation for our TDD-SBM.

6.1 Downtown Los Angeles

In Figure 3, we show the mixed-membership (TDMM-SBM) and discrete-membership (TDD-SBM) block assignments of two-block models for the downtown Los Angeles bicycle-sharing network. For the TDMM-SBM, we scale the size of each node i in our plots based on its value of $\sum_gC_{ig}$. We refer to these sums as “C total” values. These values correlate strongly with node degree (specifically, with the sum of the in-degrees and out-degrees of all layers), which is evident in the similarity of node sizes in the left and right panels of Figure 3. For both models, we observe that home and work blocks are interspersed geographically. As part of our discussion in this subsection, we will describe our method for determining the block labels in Figure 3.

The TDMM-SBM output reveals a group of stations in the left panel of Figure 3 that are neither strongly home-identified nor strongly work-identified. Instead, they possess a roughly even mixture of the two types. For this network, the TDD-SBM output is very similar to what we obtain from a discretization of the TDMM-SBM output (where we discretize by assigning each node i to the block with the maximum value of its $C_{ig}$ parameter), but this is not true for all of our bicycle-sharing networks.

Our model does not yield “home” and “work” labels for each block on its own, so we use the time-dependent block-connectivity parameter estimates $\hat{\omega}_{ght}$ to assign these labels. We assign the labels heuristically under the assumption that the “home” block is the origin of many trips to the work block in the morning and the “work” block is the origin of many trips to the home block in the evening. Figure 4, which shows $\hat{\omega}_{ght}$ for each possible value of g and h, with the hour t on the horizontal axis, supports our labeling. Given our labeling, we observe a clear peak in home-to-work traffic in the mornings and work-to-home traffic in the evenings. We make similar “home” and “work” assignments for San Francisco and New York City. In Los Angeles, the traffic within the work block peaks in the middle of the day. This may represent lunchtime errands, leisure activity, or touristic activity. (There are many tourist attractions in the downtown area.) The traffic from the home block to itself has a mild evening peak and is the lowest overall among the block pairs.

Figure 4. Estimated time-dependent block-connectivity parameters $\hat{\omega}_{ght}$ for (left) the two-block TDMM-SBM and (right) the two-block TDD-SBM for the downtown Los Angeles bicycle-sharing network.

To further evaluate our block labels, we use the zoning map for downtown Los Angeles from the Los Angeles Department of City Planning (Department of City Planning, City of Los Angeles 2015).^{Footnote 6} Zoning ordinances determine the allowable uses of city land. These ordinances distinguish between land that is available for commercial uses, industrial uses, residential uses, park districts, and other uses. In the background of Figure 5, we show a simplified version of the underlying zoning map in which we have grouped similar designations. The industrial areas include a mixture of manufacturing and commercial uses. Public facilities include government buildings, public schools, parking under freeways, and police and fire stations (Department of City Planning, City of Los Angeles 2022). In downtown Los Angeles, manufacturing and industrial areas are split cleanly from residential areas, whereas residential and commercial areas are intermixed across the bicycle-sharing system’s coverage area.

Figure 5. Two-block TDMM-SBM assignments of downtown Los Angeles bicycle-sharing stations overlaid on a simplified LA zoning map. The industrial zones include manufacturing and commercial areas. As in Figure 3, we scale the size of each node based on its value of $\sum_gC_{ig}$.

Figure 5 illustrates that most stations that are strongly home-identified are in or near zones for pure residential use or mixed residential and commercial use. We find that many stations that are not predominantly home-identified or work-identified align with mixed-use residential/commercial zones. The discrete-role plot (see the right panel of Figure 3) has a stripe of “home” stations that cut diagonally through the “work” stations. In Figure 5, we see that this aligns roughly with areas that are zoned for purely residential use. By contrast, industrial and public-facility zones tend to have stations that are primarily work-identified, although some of the most strongly work-identified stations are in mixed-used areas.

One station that seems to deviate from the overall pattern is the heavily trafficked station at Union Station. Although it is adjacent to a public facility zone with many government buildings, it is also strongly home-identified. This may seem surprising on its surface, but this classification is consistent with those of other home-identified stations because Union Station is a major transit hub for the Los Angeles metropolitan area. Many morning trips originate there as commuters transition from other modes of transportation, and many evening trips conclude there. This type of activity pattern is sensibly associated with home-identified stations. Such idiosyncrasies of transit hubs also arise in our results for San Francisco and New York City.

6.2 San Francisco

In Figure 6, we compare our two-block TDMM-SBM and two-block TDD-SBM for San Francisco. As we saw for Los Angeles, the San Francisco blocks are interspersed geographically and stations vary from strongly home-identified to strongly work-identified. The two most strongly home-identified stations are located at a major transit hub (the San Francisco Caltrain Station on 4th Street).

Figure 6. Classification of bicycle-sharing stations in San Francisco using (left) a two-block TDMM-SBM and (right) a two-block TDD-SBM. The sizes of the nodes take continuous values. In the left panel, we scale the size of each node based on its value of $\sum_gC_{ig}$. In the right panel, we scale the size of each node based on its degree.

In Figure 7, we show the estimated traffic between the “home” and “work” blocks for the two-block TDMM-SBM and the two-block TDD-SBM. As in the downtown Los Angeles network, we observe inter-block commuting. However, unlike in downtown LA, the TDD-SBM reveals intra-block morning and evening peaks in both the home and work blocks. This may be due to last-mile commuting, such as using bicycle-sharing facilities to get to or from a transit station. The 10 most-popular trips within blocks (i.e., home-to-home travel or work-to-work travel when there are two blocks) that originate during the heavy commuting hours of 7–9 am and 4–6 pm all begin at major public transportation stops in the morning and end at those stops in the evening. Recognizing last-mile usage is important for integrating bicycle-sharing systems with nearby public transportation. One possible reason that we do not observe a similar phenomenon in downtown LA is that San Franciscans are more likely than residents of Los Angeles to use public transportation (Weikel and Karlamangla Reference Weikel and Karlamangla2015). The intra-block morning and evening peaks may also arise from the intermixing of commercial and residential uses of land. In other words, some travel within blocks may also constitute commuting.

Figure 7. Estimated time-dependent block-connectivity parameters $\hat{\omega}_{ght}$ for (left) the two-block TDMM-SBM and (right) the two-block TDD-SBM for the San Francisco bicycle-sharing network.

6.2.1 Comparison to SBM variants

We briefly compare our above results for the Los Angeles and San Francisco bicycle-sharing networks to results that we obtain for these networks using a time-independent SBM and a time-dependent SBM without degree-correction. We do this to illustrate that both time-dependence and degree-correction are important for achieving our goal of detecting functional roles of bicycle-sharing stations.

To fit the time-independent SBM, we construct the time-independent adjacency matrix with entries $\tilde{A}_{ij}=\sum_{t=0}^{23}A_{ijt}$ and fit a two-block TDD-SBM to the associated time-aggregated network with a single layer. For the downtown Los Angeles network, we observe a clear geographical division in the output of the time-independent SBM (see Figure 8). For the San Francisco network, the differences between the blocks of the time-dependent SBM and time-independent SBM are less noticeable, although they are still present. This confirms that our time-dependent SBMs detect usage patterns that are not evident in the time-aggregated data.

Figure 8. Estimated blocks from two-block TDD-SBMs for time-aggregated bicycle-sharing networks in (left) downtown Los Angeles and (right) San Francisco.

We now fit the variant of our two-block TDD-SBM without degree-correction (see Section 5) to the Los Angeles and San Francisco bicycle-sharing networks. In Figure 9, we see that stations are now divided into one group of small-degree stations and one group of large-degree stations, where the “degree” is the sum over all layers of the in-degrees and out-degrees. These blocks differ significantly from those that we obtained using the TDD-SBM with degree-correction (see Figures 3 and 6). For example, for Los Angeles, the TDD-SBM with degree-correction classifies the high-traffic node at Union Station as a “home” station because of its frequent departures in the morning and frequent arrivals in the evening, whereas the TDD-SBM without degree-correction groups Union Station with other high-traffic stations that are classified predominantly as “work” stations in the degree-corrected TDD-SBM.

Figure 9. Estimated blocks from two-blocks TDD-SBMs without degree-correction for bicycle-sharing networks in (left) downtown Los Angeles and (right) San Francisco. We scale the size of each node based on its degree.

6.3 New York City

We fit our TDMM-SBM and our TDD-SBM to the New York City bicycle-sharing network with 2–6 blocks. Our results suggest that there is a limit to the size of networks for which our time-dependent SBMs can recover blocks that are based on functional roles of stations (rather than on geographical proximity). Therefore, in the present subsection, we apply our models to a subset of the New York City bicycle-sharing network to gain insight about the functional roles of stations in this subset. We select this subset based on the results of applying our three-block TDD-SBM to the entire New York City network. We examine the subnetwork that is induced by the stations that are in one of these three blocks. We explore three-block SBMs of the entire New York City network in Appendix A.3. In those results, we infer that stations are assigned to blocks along geographical lines, rather than to blocks with distinct functions.

6.3.1 Manhattan subnetwork

We examine the New York City bicycle-sharing network on a scale that is smaller than the full system and fit models to the subset of stations and trips that lie within the Manhattan (home) block of the three-block TDD-SBM (see Figure A.4) that we examine in Appendix A.3. We refer to this subnetwork as the “Manhattan” network, even though it does not include all of the stations in the borough of Manhattan. This network consists of 166 stations and 256,840 trips. In Figure 10, we show our results for a five-block TDMM-SBM and five-block TDD-SBM for the Manhattan network. The area without stations in the middle of each panel of Figure 10 is Central Park, which has stations on its perimeter but not in its interior.

Figure 10. Comparison of the blocks that we infer in the Manhattan subnetwork (i.e., the “Manhattan (home)” block in the right panel of Figure A.4) of the New York City bicycle-sharing network using (left) a five-block TDMM-SBM and (right) a five-block TDD-SBM.

The blocks that we infer using the five-block TDMM-SBM and five-block TDD-SBM outline similar regions. These models yield block-membership parameters that capture residential and commercial zones much better than the three-block TDMM-SBM and TDD-SBM do for the full New York City network. One can see this by comparing the five-block results with the underlying zoning map of the area in Figure A.6. The stations in residential zones tend to have larger block-membership parameters for “home” blocks than for “work” blocks, and the opposite is true for stations in commercial zones. We label the five detected blocks as (clockwise from the top left in each panel of Figure 10) “home (west),” “park,” “home (east),” “work,” and “mixed.” We base these labels on the land use of the underlying areas and the time-dependent block-connectivity parameters $\hat{\omega}_{ght}$ (see Figure 11). When we use the TDMM-SBM, the regions tend to blend together; they have sharp boundaries in the TDD-SBM (see Figure 10).

Figure 11. Estimated time-dependent block-connectivity parameters $\hat{\omega}_{ght}$ for (left) the five-block TDMM-SBM and (right) the five-block TDD-SBM for the Manhattan subnetwork of the New York City bicycle-sharing network.

We highlight the appearance of the “park” block, which we have not observed in previous examples and has distinctive usage patterns. The park block is similar to a residential block in its spike in morning traffic to the work block and its spike in evening traffic from the work block, but it has distinctive intra-block activity that peaks in the afternoon. The intra-block activity resembles weekend activity in the New York City bicycle-sharing system as a whole (see Figure 1); this reflects leisure use of the bicycles. Bicycles near Central Park (which also places them near several major museums) are likely to be used by tourists and other non-commuters during the day for leisure or for travel to nearby attractions.

In Figure 11, we show the values of the block-connectivity parameters $\hat{\omega}_{ght}$ for the five-block TDMM-SBM and the five-block TDD-SBM. Our estimates of $\hat{\omega}_{ght}$ for these models illustrate important differences in the usage patterns of different blocks that we observe only when we employ a time-dependent model.^{Footnote 7} We observe some overlap in the time-dependent usage patterns of the blocks, so there is potential overfitting. For example, the home (east), home (west), and mixed blocks have traffic that is similar to that in blocks other than their own. However, examining the Manhattan network using models with fewer than five blocks does not clearly distinguish the “park” block of stations from other stations in residential zones.

One reason that our time-dependent SBMs for the Manhattan network perform better (with respect to detecting functionally meaningful blocks) than models that we apply to the entire New York City network (see Appendix A.3) is the dependence of station-to-station trip counts on the distances between stations. Although our SBMs correct for the overall activity of each station, they do not normalize expected edge values by the distances between stations. In a small geographical area, such as the coverage areas of the Los Angeles and San Francisco networks, this is a reasonable choice, as all stations are within “biking distance” of each other. However, when examining a system as large as New York City’s, the lack of distance correction weakens the functional roles that we obtain using our time-dependent SBMs. Intra-block trips dwarf inter-block trips (see Figure A.5), and it seems more reasonable to interpret each block as its own ecosystem.