Data collected from a bike-sharing system exhibit complex temporal and spatial features. We analyze shared-bike usage data collected in Seoul, South Korea, at the level of individual stations while accounting for station-specific behavior and covariate effects. We adopt a penalized regression approach with a multilayer network fused Lasso penalty. The proposed fusion penalties are imposed on networks which embed spatio-temporal linkages, and capture the homogeneity in bike usage that is attributed to intricate spatio-temporal features without arbitrarily partitioning the data. We demonstrate that the proposed approach yields competitive predictive performance and provides a new interpretation of the data.

[1]Department of Statistics, University of Seoul, South Korea. Email<EMAIL_ADDRESS>
[2]School of Mathematics, University of Bristol, UK.
Email<EMAIL_ADDRESS>
[3]School of Economics, University of Seoul, South Korea.
Email<EMAIL_ADDRESS>
Keywords: bike sharing system, fused Lasso, high dimensionality, network analysis

§ INTRODUCTION

bike-sharing systems (BSSs) have become increasingly popular in urban areas and have successfully  complemented public transportation systems in dense metropolitan cities.
In addition to its utility to bike users, the installation of BSSs has been found
to reduce the usage of automobiles [10] and thus traffic congestion and possibly green house emissions [14].
To fully realize these benefits, efficient allocation of docking stations and bike docks is essential,
which in turn requires understanding the user behaviour and characteristics of the system
based on the abundant data collected on the BSS, in addition to
other urban and environmental factors that are known to influence bike usage.

In line with the increasing popularity of BSSs,
there exists a vast literature on the analysis of bike usage patterns (for and overview, see [21];
and [9])
Below we provide a brief summary of the literature on quantitative or statistical analysis of BSS usage data which is categorized into two, following [7].
The first branch addresses the problem of clustering stations in a system based on usage patterns according to some measure of similarity [12, 27, 7].
Regarding the BSS as a network, community detection algorithms have also been adopted for this purpose [1, 3, 30].
[13] model the demand for bikes as a multivariate temporal point process
and cluster stations based on functional canonical correlations
of log-intensity functions.

The second line of research concerns the problem of
predicting the station occupancy or the state of the system at a given time.
[8] model incoming and outgoing traffic at multiple stations as a panel with variables accounting for spatial and temporal autoregressive structures.
[18] model inter-station bike transitioning
for improving the effectiveness of rebalancing operations
by predicting the station drop-off demand.
[26] adopt functional linear regression
to model the directed flow between pairs of administrative divisions that aggregate multiple stations.

In all above, it is well-documented that BSS data show temporal and spatial patterns.
To address such patterns, some previous works pre-process the datasets by, for example, aggregating stations into administrative regions [26], partitioning the data using subject-specific knowledge [8], or separately analyzing the data collected on weekdays and at weekends [18].
In complex urban environments, however, it may be difficult to find a single clustering of the data that comprehensively accounts for the patterns in the usage of the BSS, since there exist multiple approaches to produce geographical or temporal divisions according to socioeconomic characteristics, land zones, traffic infrastructure or population composition.
Then, two data points that can belong to the same division for one aspect (for example, according to administrative division) may belong to different divisions in another aspect (for example, business district versus tourist area).
Besides, collecting in-depth information about the multifaceted nature of a large metropolitan city is typically costly or even impossible.

In this paper, we analyze the hourly bike rental data collected from a BSS in Seoul, South Korea, by adopting a penalized regression modelling approach. The proposed method utilizes the whole dataset without partitioning that still addresses the complex features inherent in the temporal and spatial properties of the data.
We model the BSS at the granularity of individual stations by including station-specific parameters as well as trends and variables related to air quality and precipitation.
Such a model enjoys considerable flexibility and captures time-dependent usage patterns at individual stations, but it potentially suffers from the risk of overfitting as the number of parameters ($> 4 \times 10^4$) increases linearly with the number of stations.

To address these issues, we propose a multilayer network fused Lasso penalty which extends the fused Lasso penalty [24].
In the absence of a natural ordering among the stations, the proposed penalty imposes the penalization using a multilayer network and promotes fusion of the parameters linked by edges in the network.
In doing so, we view the BSS as a multilayer network with the stations serving as its nodes.
In the network, the layers correspond to different hours of a day and within each layer, a pair of stations are connected based on their geographical distance as a proxy for shared nodal features.
The cross-layer edges come from that the usage patterns tend to vary over the course of a day.

The model fitted from the penalized regression method adaptively captures spatial and temporal homogeneity in bike usage, without (arbitrarily) partitioning the data which potentially leads to information loss.
The degree of homogeneity is determined by the data-driven choice of the penalty parameter which does not involve the researcher's subjective decision.
Our data analysis shows the superior predictive performance of the proposed multilayer fused Lasso penalty over alternative penalization approaches.
Also, we propose a new network-based model complexity measure which reveals that while similarities exist, the stations exhibit fair amount of heterogeneity.
This conclusion supports that partitioning the stations into a handful of clusters may be inappropriate for such large-scale urban transportation systems.
The proposed idea of treating the underlying structure as a multilayer network, can be useful in problems beyond that of BSS modelling whenever the data exhibits network-like features.

The remainder of this paper is organised as follows.
In Section <ref>, we describe the notable properties of BSSs such as station-specific temporal patterns, spatial homogeneity, and covariate effects,
through an exploratory analysis of
bike rental data collected in Seoul, South Korea.
Section <ref> introduces our proposed Poisson regression model and
the accompanying penalization strategy that uses of fusion penalties,
and provides a network-based interpretation of the latter.
Section <ref> demonstrates the effectiveness of our proposed approach
on the examined BSS dataset.
We conclude the paper in Section <ref>, and describe the algorithms for handling large-scale datasets in Appendix.

§ EXPLORATORY ANALYSIS OF BSS DATA

In this section, we describe data collected from “Ddareungi”, a public BSS in Seoul, South Korea.
We use the hourly rental records collected from individuals between April 1, 2019 and May 30, 2019. [The dataset is available at <https://data.seoul.go.kr/>.]
By selecting the temperate months of April and May, we avoid dealing with seasonality or possibly abnormal observations due to extreme weather conditions.
We also exclude three public holidays falling in this period from our analysis.
Thus, our final dataset comprises observations from $T = 57$ days and $S = 1,505$ stations.

Since the launch of the BSS in 2010, the numbers of subscribed users, stations, and available bikes have steadily increased, as has the number of trips.
In addition to this overall trend, the bike usage data exhibits substantial temporal variations across localities, as well as dependence on precipitation and air quality, which is in line with the observations of previous BSSs literature.

Hourly bike rentals from four selected stations over eight weeks, adjusted by the capacity (number of docking stations) of individual stations. Here, stations ST-767 and ST-816 are adjacent to one another and so are ST-107 and ST-108. We note that Week 5 contains three public holidays and show some anomalies.

Bike rental patterns exhibit substantial heterogeneity across stations.
In Figure <ref>, we plot the hourly bike rental frequencies from selected pairs of stations adjacent to one another over the examined two-month period.
There are clear station-specific patterns over the course of a day and a week, and nearby stations tend to display similar usage patterns.
Specifically, as stations ST-767 and ST-816 are located in a commercial district with government agencies and large firms, a spike in bike rental frequency is observed during weekday commuting-time.
On the other hand, stations ST-107 and ST-108 are located along a major riverside park in Seoul, and hence, are typically used for leisure activities, as evidenced by the large number of rentals concentrated on weekends.
Additionally, although daily peaks and troughs in usage can be observed, these temporal patterns do not undergo abrupt changes in the sense that the number of bikes rented out between 9am and 10am is reasonably close to that between 10am and 11am.

Average daily bike rentals at individual stations on weekdays (left), and weekends (middle) and their station-wise differences (right). The differences are taken in log-scale for better visualisation.

Difference in average daily bike rentals (in log-scale) at individual stations due to precipitation (left) and air quality (right).

Figure <ref> plots the daily rental frequencies averaged over weekdays and weekends as well as their differences (in log-scale) at individual stations.
The figure shows that bike usage behavior differs between weekdays and weekends, and that the degree of variation differs across stations.
Additionally, the last panel of Figure <ref> exhibits the presence of local clusters that share similar weekday/weekend variations.

Figure <ref> shows that bike usage depends on the weather condition and air quality.
We record each day as “rainy” if positive precipitation is recorded in any part of the city at any time of the day.
Comparing the average rental frequencies on days with and without precipitation (in log-scale, see the left panel of Figure <ref>), precipitation reduces bike usage, as expected.
Similar observations can be made with respect to the air quality: We plot the differences in average rental frequencies (in log-scale) on days when the air quality was “good” and “very bad”, in the right panel of Figure <ref>.
The qualitative categorization into “good”, “average”, “bad,” and “very bad” follows  classification system of the  Korean Ministry of Environment based on PM10 and PM2.5 dust concentrations.
The relative rental frequencies across the BSS display a similar pattern on “good” and “very bad” days, but the volume of bike rentals is lower on “very bad” days.

In summary, the BSS dataset exhibits the following characteristics:

*  There exists an overall increasing trend.

*  Bike rentals show station-specific temporal patterns and the transition is reasonably smooth.

*  These patterns are shared across stations that are geographically adjacent to one another.

*  Bike rentals are influenced by the weather condition and the air quality.

§ MODEL AND ESTIMATION

Motivated by the observations made in <ref>–<ref> in Section <ref>,
we propose a Poisson regression model
for hourly rental frequencies collected from the entire BSS (Section <ref>).
Section <ref> presents the accompanying estimation strategy
and introduces the multilayer fused Lasso penalty,
which is designed to capture the characteristics <ref>–<ref>.
In Section <ref>, we provide a network interpretation of the proposed penalization technique which aids in understanding and visualizing the penalty.

§.§ Poisson regression model

Let $Y_i$ denote the $i$th observation representing the hourly rental frequency
at station $S(i) \in \mc S = \{1, \ldots, S\}$
and hour $H(i) \in \mc H = \{0, \ldots, 23\}$,
on day $D(i) \in \mc D = \{\text{Mo, Tu, \ldots, Su}\}$
for $i = 1, \ldots, n$, with $n = 2,058,840$ denoting the number of observations in the bike usage dataset.
We denote
the time index of the $i$th observation by $t(i) \in \mc T = \{0, \ldots, T - 1\}$.
Since we consider the period of two months, as described in Section <ref>
(excluding three public holidays), we have $T = 57$.
Additionally, $z^{\rain}_i \in \{0, 1\}$ and $\mbf z^{\air}_i \in \{0, 1\}^4$
denote the variables representing
the precipitation and air quality statuses associated with the $i$th observation, respectively.
Consistent with the literature on BSSs and with the observations made in Section <ref>,
bike usage depends on the aforementioned variables
and we collect the variables relevant for the $i$th observation
in $\mbf x_i = (S(i), t(i), D(i), H(i), z^{\rain}_{i}, (\mbf z^{\air}_{i})^\top)^\top$ (for an overview, see Table <ref>).

Variables in the BSS dataset.

Variable     Range     Description

(lr)1-1 (lr)2-2 (lr)3-3
$S(i)$     $\mc S = \{1, \ldots, S\}$     Station node for the $i$th observation (categorical)

$H(i)$     $\mc H = \{0, \ldots, 23\}$     Hour of a day for the $i$th observation (categorical)

$D(i)$     $\mc D = \{\text{Mo, \ldots, Su}\}$     Day of a week for the $i$th observation (categorical)

$z^{\rain}_i$     $\{0, 1\}$     Whether it rained or not for the $i$th observation (categorical)

$\mbf z^{\air}_i$     $\{0, 1\}^4$     Whether air quality was very bad, bad, average

or good for the $i$th observation (categorical)

$t(i)$     $\{0, \ldots, T - 1\}$     Time point for the $i$th observation (numeric)

Each station in the BSS of Seoul, South Korea, has a fixed number of docks, but unlike other BSSs, this does not determine the capacity of a station, as it is possible to leave bikes even if the docks are fully occupied through chaining them to existing bikes.
Therefore, we adopt a Poisson distribution for modelling $Y_i$, the hourly count of the bikes rented out at station $S(i)$, and model the relationship between $Y_i$ and $\mbf x_i$
via Poisson regression <cit.> as

\begin{align*}
Y_i \vert \mbf x_i & \sim_{\iid} \text{Poisson}(\mu_i) \quad \text{with} \quad
\mu_i \equiv \mu(\mbf x_i) = \E(Y_i \vert \mbf x_i).
\end{align*}

One way to model the link between $\mu_i$ and $\mbf x_i$ is to use the following linear model:

\begin{align}
\log\l(\frac{\mu_i}{C_{S(i)}}\r) = \theta_{S(i)} + \alpha\, t(i) +
\beta^{\rain}\, z^{\rain}_{i} + \langle \bm\beta^{\air}, \mbf z^{\air}_{i} \rangle +
\theta^{\hod}_{H(i)} + \theta^{\dow}_{D(i)}
\label{eq:model:no}
\end{align}

where $\bm\beta^{\air} = (\beta^{\air}_j, \, 0 \le j \le 3)^\top$.
Here, $\theta^{\hod}_h$ and $\theta^{\dow}_d$, respectively, contain “hour of a day” and “day of a week” effects common to all stations.
An offset term $C_s$ relates to the capacity of station $s$, that is, the number of docks,
so that (<ref>) can be interpreted as modelling
the expected rental frequency per hour per station capacity.
This model, referred to as the no-interaction model,
does not permit the temporal effects to be station-specific
and thus is too simple to address <ref>.

Allowing for interactions between the station and temporal effects,
we consider the following full-interaction model:

\begin{align}
\log\l(\frac{\mu_i}{C_{S(i)}}\r) =& \ \theta_{S(i)} + \alpha\, t(i)
+ \beta^{\rain}\, z^{\rain}_{i}
+ \langle \bm\beta^{\air}, \mbf z^{\air}_{i} \rangle
\nn \\
& + \theta^{\hod}_{H(i)} + \theta^{\dow}_{D(i)}
+ \theta^{\hod}_{S(i), H(i)} + \theta^{\dow}_{S(i), D(i)}.
\label{eq:model}
\end{align}

The station-hour and station-day interaction terms $\theta^{\hod}_{s, h}$ and $\theta^{\dow}_{s, d}$
allow each station to exhibit individual temporal patterns.
For model identifiability, we set the baseline parameters to zero;
in other words, $\beta^{\air}_0 = \theta^{\hod}_0 = \theta^{\dow}_{\text{Mo}} = \theta^{\hod}_{s, 0} = \theta^{\dow}_{s, \text{Mo}} = \theta^{\hod}_{1, h} = \theta^{\dow}_{1, d} = 0$.

Mmodel (<ref>) accounts for <ref>, <ref>, and to a certain extent, <ref>
by including the parameters $\alpha$ capturing the overall trend,
$\beta^{\rain}$ and $\bm\beta^{\air}$ capturing the effects of precipitation and air quality,
and $\theta^{\hod}_{s, h}$ and $\theta^{\dow}_{s, d}$ addressing station-specific temporal patterns.
In doing so, we take a different approach from those taken in previous studies in which,
after (arbitrarily) partitioning the dataset according to temporal or spatial variables, or both,
individual partitions are separately modelled (e.g., [1]; [18]).
Instead, by including the interaction terms, we use the full dataset and avoid any information loss from data partitioning.

While the full-interaction model enjoys considerably more flexibility than the no-interaction model in (<ref>),
it suffers the risk of overfitting the data with the number of parameters to be estimated amounting to
$p = 34 + S + (S-1)\cdot(23 + 6) = 45,155$.
In addition, the characteristic identified in <ref>,
that is,  stations geographically close to one another tending to exhibit similar usage patterns,
is not adequately accounted for by the model fitted without any constraint.
Such an approach does not benefit from the temporal ordering inherent in the parameters $\theta^{\hod}_{s, h}$, and thus does not fully account for <ref>.
In the next section, we propose a penalized maximum likelihood estimation (MLE) methodology for model in (<ref>) with a multilayer network fused Lasso which explicitly sets out to address these issues.

§.§ Penalized MLE via multilayer network fused Lasso

We bridge the two models (<ref>) and (<ref>) at extreme ends, by adopting a penalized MLE strategy with a Lasso penalty mposed on the interaction parameters, and the fusion penalties designed to capture spatial and temporal homogeneity in bike usage patterns observed in the data, that is, <ref> and <ref>.

First proposed by [24], the fused Lasso augments
the $\ell_1$-penalized least squares estimation method–or Lasso [23]–
with a penalty that takes advantage of a meaningful ordering of the variables when such is available.
This fusion penalty is imposed on the $\ell_1$-norm of successive differences in the parameters corresponding to the ordered variables and encourages local constancy therein.
Its use has been extended beyond regression problems
such as trend estimation [25],
change point detection [15] and graphical modelling [6] among others.

Under full-interaction model in (<ref>),
we partition the parameters into
$\Theta = \{\theta_s, \, s \in \mc S\}$,
$\Theta_H = \{\theta^{\hod}_{h}, \, \theta^{\hod}_{s, h}, \, h \in \mc H \setminus \{0\}, \, s \in \mc S \setminus \{1\} \}$,
$\Theta_D = \{\theta^{\dow}_{d}, \, \theta^{\dow}_{s, d}, \, d \in \mc D \setminus \{\text{Mo}\}, \,  s \in \mc S \setminus \{1\} \}$ and
$\Delta = \{\alpha, \beta^{\rain}, \beta^{\air}_1, \beta^{\air}_2, \beta^{\air}_3\}$.
We adopt the fusion penalty to pool the information

*  from adjacent stations for the estimation of $\Theta$, $\Theta_H$ and $\Theta_D$, and

*  over the course of a day for the estimation of $\Theta_H$.

We also impose a standard lasso penalty to encourage sparsity. Subsequently, our objective is to minimize the following penalized negative log-likelihood

\begin{align}
- \ell(\Theta, \Theta_H, \Theta_D, \Delta)
&+ \lambda
\sum_{s \in \mc S} \l( \sum_{h \in \mc H} \vert \theta^{\hod}_{s, h} \vert
+ \sum_{d \in \mc D} \vert \theta^{\dow}_{s, d} \vert \r)
+ \lambda_N \cdot p_N(\Theta, \Theta_H, \Theta_D)
+ \lambda_H \cdot p_H(\Theta_H)
\label{eq:org_obj_admm}
\end{align}

with respect to $\Theta$, $\Theta_H$, $\Theta_D$ and $\Delta$.
Here, $\ell(\Theta, \Theta_H, \Theta_D, \Delta)$ denotes
the log-likelihood

\begin{align*}
\ell(\Theta, \Theta_H, \Theta_D, \Delta)
= -\sum_{i = 1}^n \mu_i + \sum_{i = 1}^n y_i \log(\mu_i) + \text{constant},
\end{align*}

and $\mu_i = \mu_i(\Theta, \Theta_H, \Theta_D, \Delta)$
is a function of the parameter vectors under (<ref>).
The parameter $\lambda$ controls the degree of the Lasso penalization
and $\lambda_N$ and $\lambda_H$ control that of the fusion penalization;
below we discuss the construction of
the functions $p_N$ and $p_H$
to achieve the above goals <ref> and <ref>.

To capture the similarities between geographically adjacent stations,
we define a set of neighboring stations for each station $s$
as $\mc N_r(s) = \{s^\prime \in \mc S \setminus \{s\}: \, d(s, s^\prime) < r\}$ and its cardinality by $\vert \mc N_r(s) \vert$,
where $d(s, s^\prime)$ denotes the distance between the two stations $s$ and $s^\prime$
measured as the crow flies, and $r$ denotes a pre-specified distance.
We define

\begin{align}
& p_N(\Theta, \Theta_H, \Theta_D) = \sum_{s \in \mc S} \sqrt{ \vert \mc N_r(s) \vert \, \sum_{s^\prime \in \mc N_r(s)} \l[ \sum_{h \in \mc H} (\phi^{\hod}_{s, h} - \phi^{\hod}_{s', h})^2 +
\sum_{d \in \mc D} (\phi^{\dow}_{s, d} - \phi^{\dow}_{s^\prime, d})^2 \r] }
\label{def:admm_net_penalty}
\end{align}

with $\phi^{\hod}_{s, h} = \theta_s + \theta^{\hod}_{h} + \theta^{\hod}_{s, h}$ and $\phi^{\dow}_{s, d} = \theta_s + \theta^{\dow}_{d} + \theta^{\dow}_{s, d}$
which fully encode the station-specific hourly and daily patterns under (<ref>).

The penalty $p_N$ forces the pairs of parameters $(\phi^{\hod}_{s, h}, \phi^{\hod}_{s', h})$ and $(\phi^{\dow}_{s, d}, \phi^{\dow}_{s^\prime, d})$ to fuse neighboring stations $s$ and $s'$,
which encourages them to exhibit similar hourly and daily patterns and thus addresses the behavior noted in <ref>.

Through $p_N$, the estimates of $\phi^{\hod}_{s, h}$ and $\phi^{\hod}_{s', h}$ may fuse even when $s' \notin \mc N_s(r)$
as long as they are connected via intermediate station nodes,
and the same applies to the estimates of $\phi^{\dow}_{s,d}$ and $\phi^{\dow}_{s',d}$.
This may be interpreted that by adopting $p_N$, we pool information across the BSS to estimate the interaction parameters (for a further discussion, see Section <ref>).
The penalty in (<ref>) is related to the spatial fusion penalties [22, 17, 20] adopted in the spatial regression literature, but is clearly distinguished as we adopt the group Lasso [29] in its formulation.
Every parameter inside the square root is considered to belong to the same group,
and the weighting applied with the size of $\mc N_r(s)$
follows the convention in the group Lasso literature.
Then, the parameters associated with the stations with a large number of neighbors receive more penalization.

We choose to decompose the station-specific hourly and daily effects as $\phi^{\hod}_{s, h} = \theta_s + \theta^{\hod}_h + \theta^{\hod}_{s, h}$ and $\phi^{\dow}_{s, d} = \theta_s + \theta^{\dow}_d + \theta^{\dow}_{s, d}$.
Then, the proposed method applies the Lasso penalty to the $\ell_1$-norm of $\theta^{\hod}_{s, h}$ and $\theta^{\dow}_{s, d}$ only, which gives it the interpretation of bridging between no-interaction and full-interaction models.
We choose not to impose the Lasso penalization on $\theta_s$ to avoid cancelling out station-specific usage patterns,
or $\theta^{\hod}_h$ and $\theta^{\dow}_d$, which represent overall hourly and daily patterns shared across all stations.

Recall that, for model identifiability, we set $\theta^{\hod}_{s, 0} = 0$.
Thus, we can write

\begin{align*}
\sum_{h \in \mc H} (\phi^{\hod}_{s, h} - \phi^{\hod}_{s', h})^2
= (\theta_s - \theta_{s'})^2 + \sum_{h = 1}^{23} (\theta^{\hod}_{s, h} - \theta^{\hod}_{s', h})^2
\end{align*}

and similarly, since $\theta^{\dow}_{s, \text{Mo}} = 0$, we have

\begin{align*}
\sum_{d \in \mc D} (\phi^{\dow}_{s, d} - \phi^{\dow}_{s', d})^2
= (\theta_s - \theta_{s'})^2 + \sum_{d \in \mc D \setminus \{\text{Mo}\}} (\theta^{\dow}_{s, d} - \theta^{\dow}_{s', d})^2.
\end{align*}

From these observations, we can re-write $p_N$ as

\begin{align*}
& p_N(\Theta, \Theta_H, \Theta_D) =
\sum_{s \in \mc S} \sqrt{ \vert \mc N_r(s) \vert \, \sum_{s^\prime \in \mc N_r(s)} \l[ 2(\theta_s - \theta_{s^\prime})^2 + \sum_{h = 1}^{23} (\theta^{\hod}_{s, h} - \theta^{\hod}_{s', h})^2 +
\sum_{d \in \mc D \setminus \{\text{Mo}\}} (\theta^{\dow}_{s, d} - \theta^{\dow}_{s^\prime, d})^2 \r] },
\end{align*}

which shows that $p_N$ implicitly encourages the parameters
$(\theta_s, \theta_{s'})$ to take similar values when stations $s$ and $s'$ are neighbors.

There is a natural temporal ordering inherent in $\Theta_H$ that gives rise to
the fusion penalty in its canonical form imposed on $\phi^{\hod}_{s, h}$:

\begin{align}
p_H(\Theta_H) =  \sum_{s \in \mc S}
\sum_{h = 0}^{23} \l\vert \phi^{\hod}_{s, h} - \phi^{\hod}_{s, h + 1} \r\vert
=  \sum_{s \in \mc S}
\sum_{h = 0}^{23} \l\vert (\theta^{\hod}_{h} +\theta^{\hod}_{s, h})- (\theta^{\hod}_{h+1} + \theta^{\hod}_{s, h + 1}) \r\vert
\label{def:admm_hour_penalty}
\end{align}

with $\theta^{\hod}_{s, 24} = \theta^{\hod}_{s, 0}$.
Imposing a penalty on $p_H(\Theta_H)$ encourages the consecutive station-specific hourly effects
$\phi^{\hod}_{s, h}$ and $\phi^{\hod}_{s, h+1}$ at a given station $s$,
to become close to one another and suppress abrupt changes in usage.
This reflects the fact that rental counts at each station rarely undergo radical shifts
during the course of the day.

Jointly, $p_N$ and $p_H$ comprise the proposed multilayer network fusion penalty.
We devote Section <ref> to its interpretation with the description of the underlying multilayer network.
The impact of the fusion penalty is determined by the sizes of $\lambda_N$ and $\lambda_H$,
which we select via cross-validation, as described in Section <ref>.
The Alternating Direction Method of Multipliers (ADMM) algorithm [4] is employed to solve the convex optimization problem in (<ref>).
Efficient implementation of the algorithm requires careful re-parametrization of  model (<ref>), which makes use of the specific structure of the data, and we discuss this in detail in Appendix <ref>.

§.§ multilayer network interpretation of the fusion penalty

We provide a network interpretation of the proposed penalized regression method.

First, we introduce the following networks that underpin the penalty functions $p_N$ and $p_H$:

\begin{align}
\mathfrak{N}_{\text{single}}(r) =& (\mc S, \mc E_{\text{single}}(r))
\text{ \ with \ }
\mc E_{\text{single}}(r) = \cup_{s \in \mc S} \{ (s, s'), \, s' \in \mc N_r(s)\},
\label{eq:net:single} \\
\mathfrak{N}_{\text{multi}}(r) =&
(\mc S \times \mc H, \, \mc E_{\text{multi}}(r)) \text{ \ with \ } \nn
\\
\mc E_{\text{multi}}(r) =&
\l[ \cup_{h \in \mc H} \cup_{s \in \mc S} \l\{ \l( (s, h), (s', h) \r), \, s' \in \mc N_r(s) \r\} \r] \bigcup \l[\cup_{s \in \mc S} \l\{ \l( (s, h), (s, h + 1) \r), \, h \in \mc H \r\}\r].
\label{eq:net:multi}
\end{align}

(For an illustrative example of $\mathfrak{N}_{\text{single}}$ and $\mathfrak{N}_{\text{multi}}$, see Figure <ref>).
While both $\mathfrak{N}_{\text{single}}$ and $\mathfrak{N}_{\text{multi}}$ as well as their edge sets
depend on the choice of $r$, we suppress this dependency for simplicity when it does not cause any confusion.

Illustration of a single-layer network $\mathfrak{N}_{\text{single}}$ (left) and a multilayer network $\mathfrak{N}_{\text{multi}}$ (right). Each dot represents a node (station) and a solid line represents an edge connecting the nodes within each layer and across adjacent layers.

Network $\mathfrak{N}_{\text{single}}$ is a single-layer, undirected network
that is solely determined by the sets of neighbors $\mc N_r(s)$ of the stations.
In this network, a pair of “day of a week” parameters $(\phi^{\dow}_{s, d}, \phi^{\dow}_{s', d})$ for each given day $d$, are encouraged to take values close to one another by the penalty function $p_N$, provided that the stations belong to the same connected component of $\mathfrak{N}_{\text{single}}$.
On the other hand, $\mathfrak{N}_{\text{multi}}$ is a multilayer, undirected network;
we follow the notational convention of [16], in which the set $\mc H$ serves as a set of elementary layers for the hourly aspect,
and each edge connects a pair of node-layer tuples $(s, h)$ and $(s', h')$ for some $s, s' \in \mc S$ and $h, h' \in \mc H$.
Each pair of the “hour of a day” parameters
$\phi^{\hod}_{s, h}$ and $\phi^{\hod}_{s', h'}$

is encouraged to fuse with one another by the penalty functions $p_N$ and $p_H$,
if the corresponding pair of nodes are connected in $\mathfrak{N}_{\text{multi}}$.

Next, we define networks whose edges
are determined by the coefficient estimates from the penalized MLE.
$\wh\phi^{\hod}_{s, h}$ and $\wh\phi^{\dow}_{s, d}$ denote the estimates of the parameters $\phi^{\hod}_{s, h}$ and $\phi^{\dow}_{s, d}$, respectively;

their values depend on the tuning parameters
$(r, \lambda, \lambda_N, \lambda_H)$; however,
we omit this dependence for simplicity.
Then, the networks associated with these estimates are defined as

\begin{align}
\wh{\mathfrak{N}}_{D, d} = \l(\mc S, \, \mc E_{D, d}\r) \text{ \ with \ }
\mc E_{D, d} =& \l\{(s, s'), \, s \ne s': \, \wh\phi^{\dow}_{s, d} = \wh\phi^{\dow}_{s', d} \r\} \text{ \ for each \ } d \in \mc D,
\nn \\
\wh{\mathfrak{N}}_H = \l(\mc S \times H, \mc E_H\r)  \text{ \ with \ }
\mc E_H =&  \cup_{h \in \mc H} \cup_{s \in \mc S}
\l\{((s, h), (s', h')), \, (s,h) \ne (s',h'): \, \wh\phi^{\hod}_{s, h} = \wh\phi^{\hod}_{s', h'} \r\}
\label{eq:net:hod}
\end{align}

As with $\mathfrak{N}_{\text{single}}$, the networks $\wh{\mathfrak{N}}_{D, d}$ are single-layer networks
and an edge joins two station nodes $s$ and $s'$ when their node features (i.e.  parameter estimates of  $\phi^{\dow}_{s,d}$ and $\phi^{\dow}_{s',d}$ for a given $d \in \mc D$) are identical, possibly due to the fusion penalty but not necessarily so.
The network $\wh{\mathfrak{N}}_H$, as with $\mathfrak{N}_{\text{multi}}$,
is a multilayer network with the hourly layer given by $\mc H$,
and an edge is formed between a pair of nodes $(s,h)$ and $(s',h')$ when the estimates of $\phi^{\hod}_{s, h}$ and $\phi^{\hod}_{s', h'}$ agree at $(s,h) \ne (s',h')$.
(For an illustrative example of $\wh{\mathfrak{N}}_H$, see Figure <ref>).

Illustration of multilayer networks.
Left: Each layer of $\mathfrak{N}_{\text{multi}}$ (see (<ref>)) embeds the linkages between the stations determined by their geographical distances at a given hour $h$, and the station at layer $h$ is linked to itself at layers $h - 1$ and $h + 1$,
which underpins how $\phi^{\hod}_{s, h}$ are encouraged to be fused by $p_N$ and $p_H$.
Middle: $\wh{\mathfrak{N}}_H$ (see (<ref>)) is determined by the values of the estimates of $\phi^{\hod}_{s, h}$ with an edge indicating that the connected estimates share the identical values.
Right: A multilayer network formed with its edges obtained as an intersection of the edge sets of $\mathfrak{N}_{\text{multi}}$ and $\wh{\mathfrak{N}}_H$,
which contains four connected components.

For two networks (either single- or multilayer)
$\mathfrak{N}_i = (\mc V, \mc E_i), \, i = 1, 2$, sharing the same node set $\mc V$,
denote by $\mathfrak{N}_1 \cap \mathfrak{N}_2 = (\mc V, \mc E_1 \cap \mc E_2)$
the network formed by taking the intersection of their edge sets.
Our proposed penalized regression method takes as an input the observable networks $\mathfrak{N}_{\text{single}}$ and $\mathfrak{N}_{\text{multi}}$, which inform spatial and temporal proximity.
Then, it outputs the networks capturing the homogeneity between the stations nodes, namely, $\mathfrak{N}_{\text{single}} \cap \wh{\mathfrak{N}}_{D, d}$ (on a given day of a week $d$) and $\mathfrak{N}_{\text{multi}} \cap \wh{\mathfrak{N}}_H$ (along the hourly layer),
which can provide insights into BSS management and urban planning.
Depending on the choice of penalty parameters, the output networks are not necessarily sparse; in fact, this is the case in our data analysis reported in Section <ref>.
This distinguishes our approach from the existing literature on clustering or partitioning the dataset using spatial or temporal variables prior to analysis.

§ DATA ANALYSIS

§.§ Tuning parameter selection

For the selection of $r$ that determines the set of neighbors $\mathcal{N}_r(s)$ for each station $s$, we examine the distance between each station and its $10$ nearest stations,
that is, $d_{s, i}, \, i = 1, \ldots, 10$, for all $s \in \mc S$.
The median of the first quartile

of $d_{s, i}$ over $s \in \mc S$ is $762$ meters,

from which we choose to consider $r \in \{750, 1500\}$.
The network $\mathfrak{N}_{\text{single}}(r)$ (see (<ref>)) with
$r = 750$ has $34$ connected components, whereas the network with $r = 1500$ has two connected components.
Later, we show that our penalized regression approach
is not sensitive to $r$ on the BSS dataset
through the adaptive selection of the penalty parameters discussed below.

The parameters $\lambda$, $\lambda_N$, and $\lambda_H$
control the overall complexity of the fitted model, and
the latter two, in particular control the degree of spatial and temporal homogeneity
induced by the parameter estimates of $\Theta$, $\Theta_H$ and $\Theta_D$.

We propose the selection of the tuning parameters via cross validation (CV).
In the penalized regression literature, CV is typically performed by randomly partitioning the data into five or ten folds,

which, in the case of the BSS dataset, ignores the temporal structure therein.
Instead, we make use of the fact that our dataset covers seven weekends,
and adopt seven-fold CV, in which each fold is ensured to
include a balanced number of all seven days of the week.

As a CV measure, we adopt the mean squared Pearson residuals (MSPR):

\begin{align}
\label{eq:cv}
\text{CV}(r, \lambda, \lambda_N, \lambda_H) =
\frac{1}{7} \sum_{j = 1}^7 \frac{1}{n_j} \sum_{i = 1}^{n_j}
\frac{(Y^{(j)}_i - \wh\mu^{(j)}_i(r, \lambda, \lambda_N, \lambda_H))^2}{\wh\mu^{(j)}_i(r, \lambda, \lambda_N, \lambda_H)}
\end{align}

where for the $j$th fold,
$n_j$ denotes the total number of observations,
$Y_i^{(j)}$ the $i$th observation and
$\wh\mu_i^{(j)}(r, \lambda, \lambda_N, \lambda_H)$ is the corresponding estimate of the mean from the model fitted to the remaining data with the given tuning parameters.
We evaluate $\text{CV}(r, \lambda, \lambda_N, \lambda_H)$ on grids of equispaced values on a log-scale for $(\lambda, \lambda_N, \lambda_H)$
and select the combination that returns the smallest CV.

§.§ Analysis of the BSS data

§.§.§ Effects of penalization

To assess the effect of penalization, we compare our multilayer network fused Lasso-based penalized regression approach (referred to as “fused Lasso”) with the methods that adopt either the fusion penalty or the Lasso penalty only, which are referred to as “fusion-only” and “Lasso-only”.
We also consider models (<ref>) and (<ref>) fitted without any penalization,
referred to as “no-interaction” and “full-interaction”, respectively.

Predictive performance

Due to small sample size and inherent non-exchangeability of the data, we adopt CV measures as an indicator of the predictive performance.
In Table <ref> and Figure <ref>,
we report the fold-wise MSPRs involved in the seven-fold CV in (<ref>) as well as the overall CV error, all evaluated at the penalty parameters selected to minimize the overall CV error for the respective methods.

When comparing no-interaction and full-interaction methods,
the flexibility afforded by allowing for station-specific temporal effects proves useful in enhancing the predictive performance as the latter model consistently attains a considerably smaller MSPR.
We observe further improvement when appropriate penalization is applied to the interaction parameters.
In particular, adopting the proposed fused Lasso penalty returns the minimum CV error across the seven folds compared to the other methods regardless of the choice of $r$.
Between the fusion-only and the Lasso-only methods, the former outperforms the latter, indicating that capturing across-station homogeneity pays off by accounting for the stylized features of the BSS data that neighboring stations exhibit similar usage patterns, (see <ref>).
The MSPR attained with $r = 1500$ is consistently smaller, albeit by a small margin, than that attained with $r = 750$.
Recalling that the choice $r = 1500$ reduces the number of connected components in $\mathfrak{N}_{\text{single}}(r)$ compared to when $r = 750$ (from $54$ to $2$), this result is indicative of the benefit of pooling information across the wider spatial networks of BSSs.

Mean squared Pearson residuals (MSPRs) from each fold used in the seven-fold CV and the overall CV error as measured in (<ref>), with the penalty parameters selected to minimize the latter for the respective methods where relevant.

7cFold-wise CV     2c

Method     1     2     3     4     5     6     7     Average

(lr)1-1 (lr)2-8 (lr)9-10
Fused Lasso (750m)     1.558     1.423     1.424     1.497     1.449     1.831     1.683     1.552

Fused Lasso (1500m)     1.553     1.420     1.420     1.494     1.443     1.826     1.675     1.547

(lr)1-1 (lr)2-8 (lr)9-10
Fusion-only (750m)     1.572     1.443     1.445     1.512     1.473     1.861     1.701     1.569

Fusion-only (1500m)     1.568     1.440     1.438     1.511     1.463     1.847     1.700     1.567

Lasso-only     1.609     1.479     1.473     1.537     1.493     1.902     1.728     1.603

Full-interaction     1.729     1.601     1.556     1.629     1.596     2.016     1.841      1.710

(lr)1-1 (lr)2-8 (lr)9-10
No-interaction     1.991     1.857     1.791     1.918     1.838     2.264     2.118     1.969

Mean squared Pearson residuals (MSPRs) from each fold used in the seven-fold CV and the overall CV error as measured in (<ref>),
with the penalty parameters selected to minimize the CV error for the respective methods.

Sparsity and model complexity

We examine the reduction in complexity brought about by the Lasso penalty and the fusion penalties $p_N$ and $p_H$ in (<ref>).

The sparsity induced by the Lasso penalty is easily measured by the proportion of non-zero coefficient estimates (see the last column of Table <ref>).
The model fitted with the fused Lasso penalty is approximately $13$% (resp. $20$%) less complex than the full-interaction model when $r = 750$ (resp. $r = 1500$),
yet the former has much better predictive power than the latter as  Table <ref> shows.
At the same time, further sparsity induced by adopting the Lasso-only penalty does not improve the prediction performance.
The fusion-only method does not set any parameter estimate to be exactly zero.

Model complexity (MC) of the models fitted the penalty parameters selected to minimize the CV error for the respective methods (where relevant).
We report the overall MC out of all parameters ($\Delta \cup \Theta \cup \Theta_H  \cup \Theta_D$) as well as that out of the interaction  parameters representing hour-of-a-day ($\phi^{\hod}_{s,h}$) and day-of-a-week ($\phi^{\dow}_{s,d}$) effects. Additionally, the proportion of non-zero coefficient estimates is presented.

3cParameter set     Proportion of

Method     All     $\Theta_H \setminus \{\theta^{\hod}_h, \, h \in \mc H\}$     $\Theta_D \setminus \{\theta^{\dow}_d, \, d \in \mc D\}$     non-zeros

(lr)1-1 (lr)2-2 (lr)3-3 (lr)4-4 (lr)5-5
Fused Lasso (750m)     0.828     0.790     0.977     0.867

Fused Lasso (1500m)     0.848     0.827     0.927     0.795

(lr)1-1 (lr)2-2 (lr)3-3 (lr)4-4 (lr)5-5
Fusion-only (750m)     0.903      0.885     0.974     1.000

Fusion-only (1500m)     0.853     0.836     0.916     1.000

(lr)1-1 (lr)2-2 (lr)3-3 (lr)4-4 (lr)5-5
–     –     –     0.632

We use the networks introduced in Section <ref> to define a model complexity (MC) measure that evaluates the effect of the fusion penalties.
Simply put, the model represents the proportion of coefficient estimates that are not fused by the penalties $p_N$ and $p_H$,
out of the total number of parameters in model in (<ref>), as

\begin{align}
& \text{MC}(r, \lambda, \lambda_N, \lambda_H)
=  \frac{1}{p}
\l( 34 + \mc C(\mathfrak{N}_{\text{multi}}(r) \cap \wh{\mathfrak{N}}_H) + \sum_{d \in \mc D  \setminus \{\text{Mo} \}} \mc C(\mathfrak{N}_{\text{single}}(r) \cap \wh{\mathfrak{N}}_{D, d}) \r)
% {34  + \vert \Theta_H  \vert + \vert \Theta_D  \vert },
\label{eq:mc}
\end{align}

where we suppress the dependence on $(r, \lambda, \lambda_N, \lambda_H)$ for notational simplicity.
Recall that $34$ represents the number of parameters that are not penalized, namely $\alpha$, $\beta^{\text{rain}}$, $\{\beta^{\text{air}}_j, \, 1 \le j \le 3\}$, $\{\theta^{\hod}_h, \, h \in \mc H\}$ and $\{\theta^{\dow}_d, \, d \in \mc D\}$.
Denoting the number of connected components in a network $\mathfrak{N}$ by $\mc C(\mathfrak{N})$,
the number of unique parameter estimates of $\phi^{\hod}_{s,h}$, which are not fused by penalization,
is given by $\mc C(\mathfrak{N}_{\text{multi}}(r) \cap \wh{\mathfrak{N}}_H)$.
Figure <ref> shots that the network $\mathfrak{N}_{\text{multi}} \cap \wh{\mathfrak{N}}_H$ contains four connected components; consequently,
$\mc C(\mathfrak{N}_{\text{multi}} \cap \wh{\mathfrak{N}}_H) = 4$. Similarly,
We can find the number of unique parameter estimates for $\phi^{\dow}_{s,d}$.
In (<ref>), we exclude the intersection network $\mathfrak{N}_{\text{single}}(r) \cap \wh{\mathfrak{N}}_{D, \text{Mo}}$ from the numerator since, due to model identifiability constraints, we have $\phi^{\dow}_{s,\text{Mo}} = \phi^{\hod}_{s, 0} = \theta_{s}$.
That is, the fusion among the station-specific intercept parameters $\theta_s$ has already been accounted for by $\mathfrak{N}_{\text{multi}}(r) \cap \wh{\mathfrak{N}}_H$ at layer $h = 0$.
We refer to Appendix <ref> for the efficient calculation of MC which requires some care owing to the multilayered nature of the network in (<ref>) underlying the penalties.
We analogously evaluate the complexity of the model fitted using the fusion-only method (for the full results, see Table <ref>).

The MC ranges between $0$ and $1$, and when its value is closer to $0$, it implies that most stations exhibit homogeneous behavior with their neighboring stations.
The resultant intersection networks
$\mathfrak{N}_{\text{multi}}(r) \cap \wh{\mathfrak{N}}_H$ and
$\mathfrak{N}_{\text{single}}(r) \cap \wh{\mathfrak{N}}_{D, d}$
are not highly homogeneous, as evidenced by the MC being close to one, that is,  many stations exhibit individual behavior.
In other words, station-specific parameters account for a large portion of the variation in bike usage, which supports modelling the data at the individual station level.
The high degree of heterogeneity across the BSS can be attributed to the fact that Seoul is a large city with high population density, so that each station is associated with multiple aspects of usages.
This indicates that partitioning stations into a handful of clusters may ignore the complex nodal features that drive the usage of bikes at each station.
Instead, our proposed multilayer fused Lasso method imposes fusion penalties through the networks $\mc N_{\text{single}}(r)$ and $\mc N_{\text{multi}}(r)$, and learns the degree of homogeneity in bike usage patterns across the BSS in a data-driven manner.

Generally, more homogeneity is observed in the hour-of-a-day effects than in the day-of-a-week effects.
Also, since the choice of $r = 1500$ leads to a highly connected $\mathfrak{N}_{\text{single}}(r)$, it leads to a fitted model with a larger MC than the choice of $r = 750$.
Overall, the fused Lasso method tends to further reduce the MC compared to the fusion-only method.
We attribute this to the fact that the Lasso penalty enhances the effect of the fusion penalty by coercing estimated values to be zeros.

§.§.§ Parameter estimates

There exist methods for performing inference in high-dimensional generalized linear models [2, 5], but they do not easily apply to our setting because of the presence of a fusion penalty.
Instead, we examine the coefficient estimates obtained from different folds used in the cross validation step in (<ref>) (for the estimates of the parameters associated with the overall trend and the effects of the precipitation and the air quality, obtained from each fold as well as from the full data, see Table <ref>).
We only report the results when $r = 1500$ in the main text because we obtain nearly identical estimates when $r = 750$ as shown in Table <ref> in Appendix <ref>).
Table <ref> shows that, while the values of the estimates vary slightly from one fold to another, their signs and overall magnitude do not change.
This confirms the observations <ref> and <ref> made in Section <ref>, that is, the variables have meaningful effects on overall bike usage across the system.

Estimated coefficients for the covariate effects by the proposed fused Lasso regression method from each fold used in the seven-fold CV and from the full data when $r = 1500$.
For comparison, we also report the estimates obtained with fusion Lasso-only, full-interaction and no-interaction methods.

7cFold     5cFull data

1     2     3     4     5     6     7     Fused     Fusion-only     Lasso-only     Full     No

(lr)1-1 (lr)2-8 (lr)9-9 (lr)10-10 (lr)11-11 (lr)12-12 (lr)13-13
$\alpha$     0.057     0.052     0.051     0.050     0.050     0.054     0.050      0.052     0.052     0.052     0.052     0.052

$\beta^{\text{rain}}$     -2.301     -2.356     -2.380     -2.376     -2.238     -2.182     -2.382      -2.324     -2.325     -2.321     -2.323     -2.310

$\beta^{\text{air}}_1$     0.098     0.141     0.100     0.177     0.124     0.169     0.168     0.140     0.141     0.138     0.138     0.137

$\beta^{\text{air}}_2$     0.117     0.174     0.099     0.177     0.133     0.159     0.188      0.149     0.150     0.145     0.144     0.146

$\beta^{\text{air}}_3$     0.218     0.285     0.251     0.169     0.259     0.321     0.310     0.275     0.278     0.271     0.276     0.243

Parameter estimates for $\phi^{\hod}_h, \, h \in \mc H$ (left) and $\phi^{\dow}_d, \, d \in \mc D$ (right) from each fold used in the seven-fold CV and from the full data when $r = 1500$.

Estimated station-specific bike demands in log-scale (given by $\wh\theta_s + \wh\theta^{\hod}_h + \wh\theta^{\dow}_d + \wh\theta^{\hod}_{s, h} + \wh\theta^{\dow}_{s, d}$) from the model fitted with $r = 1500$ at 8am on Tuesdays (left), at 8pm on Sundays (middle) and their differences (right).

In Figure <ref>, we plot the estimates of $\theta^{\hod}_h$ and $\theta^{\dow}_d$ which are shared by all the stations belonging to the BSS.
We observe that the smooth transition over the course of a day and a week is well captured across the seven folds, and there are peaks corresponding to the high demand by commuters.
In addition, we plot the combined effects of temporal variables on the mean bike demand (in log-sale), namely $\wh\theta_s + \wh\theta^{\hod}_h + \wh\theta^{\dow}_d + \wh\theta^{\hod}_{s, h} + \wh\theta^{\dow}_{s, d}$ under (<ref>), for all stations $s \in \mc S$ when $(h, d) = (8, \text{Tu})$ (8am on Tuesday) and $(20, \text{Su})$ (8pm on Sunday) (see Figure <ref>).
As expected the spatial distribution of bike usage concentration is markedly different when $(h, d) = (8, \text{Tu})$ and $(h, d) = (20, \text{Su})$, as  bikes are primarily used for commuting at 8am on Tuesdays and for leisure activities at 8pm on Sundays.
In Figures <ref> and <ref>, we report the results obtained with $r = 1500$ only. (For those obtained with $r = 750$, indicating that the parameter estimates are not sensitive to the choice of $r$,
we refer to Figures <ref> and <ref> in Appendix <ref>.)

§ CONCLUSIONS

In this study, we examine the problem of modelling usage data collected
from a BSS spanning a large metropolitan city.
We model the data at the granularity of individual stations
by incorporating covariate effects as well as spatial and temporal characteristics commonly observed in bike usage data.
The proposed multilayer fused Lasso penalty is imposed on the networks encoding  the geographical proximity of the stations with an additional hourly layer, and successfully captures the spatial and temporal homogeneity.
Combined with the data-driven choice of penalty parameters, our penalized regression approach strikes a good balance between a simplistic model that does not allow for station-specific behavior,
and a complex model suffering from high dimensionality,
by returning a fitted model with good predictive performance.

We envision that the proposed method is applicable to different datasets with network-like features,
such as those collected from large transportation, communication, or logistic systems.
In particular, when information about the factors driving link homophily
(such as the nodal features related to land use, slope of terrain, nearby landmarks, and other modes of transportation in the case of BSSs) is not readily available,
our penalized regression approach enables learning of the linkages in an unobservable network
from the fusion of parameters induced by the penalties defined on an observable network
(such as those determined by the distance between pairs of stations).

§.§.§ Acknowledgement

Yunjin Choi was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT, No. NRF-2020R1G1A1A01014088). Haeran Cho was supported by the Leverhulme Trust Research Project Grant (RPG-2019-390).

[1]
Austwick, M. Z., O’Brien, O., Strano, E., and Viana, M. (2013).
The structure of spatial networks and communities in bicycle sharing
PloS one, 8(9):e74685.

[2]
Belloni, A., Chernozhukov, V., and Wei, Y. (2016).
Post-selection inference for generalized linear models with many
Journal of Business & Economic Statistics, 34(4):606–619.

[3]
Borgnat, P., Robardet, C., Abry, P., Flandrin, P., Rouquier, J.-B., and
Tremblay, N. (2013).
A dynamical network view of lyon’s vélo’v shared bicycle
In Dynamics On and Of Complex Networks, Volume 2, pages
267–284. Springer.

[4]
Boyd, S., Parikh, N., and Chu, E. (2011).
Distributed Optimization and Statistical Learning via the
Alternating Direction Method of Multipliers.
Now Publishers Inc.

[5]
Cai, T. T., Guo, Z., and Ma, R. (2021).
Statistical inference for high-dimensional generalized linear models
with binary outcomes.
Journal of the American Statistical Association (to appear).

[6]
Danaher, P., Wang, P., and Witten, D. M. (2014).
The joint graphical Lasso for inverse covariance estimation across
multiple classes.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 76(2):373.

[7]
Etienne, C. and Latifa, O. (2014).
Model-based count series clustering for bike sharing system usage
mining: a case study with the vélib’system of paris.
ACM Transactions on Intelligent Systems and Technology (TIST),

[8]
Faghih-Imani, A. and Eluru, N. (2016).
Incorporating the impact of spatio-temporal interactions on bicycle
sharing system demand: A case study of new york citibike system.
Journal of Transport Geography, 54:218–227.

[9]
Fishman, E. (2016).
Bikeshare: A review of recent literature.
Transport Reviews, 36(1):92–113.

[10]
Fishman, E., Washington, S., and Haworth, N. (2014).
Bike share’s impact on car use: Evidence from the united states,
great britain, and australia.
Transportation Research Part D: Transport and Environment,

[11]
Friedman, J., Hastie, T., Tibshirani, R., Narasimhan, B., Tay, K., Simon, N.,
and Qian, J. (2020).
glmnet: Lasso and Elastic-Net Regularized Generalized Linear
R package version 4.0-2.

[12]
Froehlich, J. E., Neumann, J., and Oliver, N. (2009).
Sensing and predicting the pulse of the city through shared
In Twenty-First International Joint Conference on Artificial

[13]
Gervini, D. and Khanal, M. (2019).
Exploring patterns of demand in bike sharing systems via replicated
point process models.
Journal of the Royal Statistical Society: Series C (Applied
Statistics), 68(3):585–602.

[14]
Hamilton, T. L. and Wichman, C. J. (2018).
Bicycle infrastructure and traffic congestion: Evidence from dc's
capital bikeshare.
Journal of Environmental Economics and Management, 87:72–93.

[15]
Harchaoui, Z. and Lévy-Leduc, C. (2010).
Multiple change-point estimation with a total variation penalty.
Journal of the American Statistical Association,

[16]
Kivelä, M., Arenas, A., Barthelemy, M., Gleeson, J. P., Moreno, Y., and
Porter, M. A. (2014).
Multilayer networks.
Journal of Complex Networks, 2(3):203–271.

[17]
Li, F. and Sang, H. (2019).
Spatial homogeneity pursuit of regression coefficients for large
Journal of the American Statistical Association,

[18]
Liu, J., Sun, L., Chen, W., and Xiong, H. (2016).
Rebalancing bike sharing systems: A multi-source data smart
In Proceedings of the 22nd ACM SIGKDD International Conference
on Knowledge Discovery and Data Mining, pages 1005–1014.

[19]
McCullagh, P. and Nelder, J. A. (1989).
Generalized Linear Models.
Chapman and Hall.

[20]
Sass, D., Li, B., and Reich, B. J. (2021).
Flexible and fast spatial return level estimation via a spatially
fused penalty.
Journal of Computational and Graphical Statistics,

[21]
Shaheen, S. A., Guzman, S., and Zhang, H. (2010).
Bikesharing in europe, the americas, and asia: past, present, and
Transportation research record, 2143(1):159–167.

[22]
Sun, Y., Wang, H. J., and Fuentes, M. (2016).
Fused adaptive Lasso for spatial and temporal quantile function
Technometrics, 58(1):127–137.

[23]
Tibshirani, R. (1996).
Regression shrinkage and selection via the Lasso.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 58(1):267–288.

[24]
Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., and Knight, K. (2005).
Sparsity and smoothness via the fused Lasso.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 67:91–108.

[25]
Tibshirani, R. J. (2014).
Adaptive piecewise polynomial estimation via trend filtering.
The Annals of Statistics, 42(1):285–323.

[26]
Torti, A., Pini, A., and Vantini, S. (2021).
Modelling time-varying mobility flows using function-on-function
regression: Analysis of a bike sharing system in the city of Milan.
Journal of the Royal Statistical Society: Series C (Applied
Statistics), 70(1):226–247.

[27]
Vogel, P., Greiser, T., and Mattfeld, D. C. (2011).
Understanding bike-sharing systems using data mining: Exploring
activity patterns.
Procedia-Social and Behavioral Sciences, 20:514–523.

[28]
Wahlberg, B., Boyd, S., Annergren, M., and Wang, Y. (2012).
An admm algorithm for a class of total variation regularized
estimation problems.
IFAC Proceedings Volumes, 45(16):83–88.

[29]
Yuan, M. and Lin, Y. (2006).
Model selection and estimation in regression with grouped variables.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 68(1):49–67.

[30]
Zhou, X. (2015).
Understanding spatiotemporal patterns of biking behavior by analyzing
massive bike sharing data in chicago.
PloS one, 10(10):e0137922.

§ COMPUTATIONAL CONSIDERATIONS

§.§ ADMM algorithm for optimization of penalized MLE

§.§.§ ADMM framework

We adopt the alternating direction method of multipliers (ADMM) algorithm for fitting (<ref>) based on [28],
where the optimization of objective functions with fusion penalties is discussed.
We recall that the fused Lasso penalties lead to the objective function of the form

\begin{align}
\min_{\Theta, \Theta_H, \Theta_D, \Delta} & \sum_{i=1}^{n} \mu_i(\Theta, \Theta_H, \Theta_D, \Delta) - \sum_{i=1}^{n} y_i \log(\mu_i(\Theta, \Theta_H, \Theta_D, \Delta)) \nn \\
&+ \lambda
\sum_{s \in \mc S} \l( \sum_{h \in \mc H} \vert \theta^{\hod}_{s, h} \vert
+ \sum_{d \in \mc D} \vert \theta^{\dow}_{s, d} \vert \r)
\nn \\
&+ \lambda_N \sum_{s \in \mc S} \sqrt{ \vert \mc N_r(s) \vert \, \sum_{s^\prime \in \mc N_r(s)} \l[ 2(\theta_s - \theta_{s^\prime})^2 + \sum_{h = 1}^{23} (\phi^{\hod}_{s, h} - \phi^{\hod}_{s', h})^2 +
\sum_{d \in \mc D \setminus \{\text{Mo}\}} (\phi^{\dow}_{s, d} - \phi^{\dow}_{s^\prime, d})^2 \r]}  \nn \\
&+ \lambda_H \sum_{s \in \mc S}
\sum_{h = 0}^{23} \l\vert \phi^{\hod}_{s, h} - \phi^{\hod}_{s, h + 1} \r\vert.
% =  \sum_{s \in \mc S} \sum_{h = 0}^{23} \l\vert (\theta^{\hod}_{h} +\theta^{\hod}_{s, h})- (\theta^{\hod}_{h+1} + \theta^{\hod}_{s, h + 1}) \r\vert
\label{admmeq:org_obj_admm}
\end{align}

In this section, we treat $\Theta$, $\Theta_{H}$ and $\Theta_{D}$ as row-vectors without confusion:

\begin{align*}
\Theta =& (\theta_s, \, s \in {\mc S}), \\
\Theta_H =& \left((\theta_{s,1}^{\hod}, \, s\in {\mc S}), \ldots, (\theta_{s, 23}^{\hod}, \, s\in {\mc S}), \, (\theta_1^{\hod}, \ldots, \theta_{23}^{\hod}) \right), \\
\Theta_D =& \left((\theta_{s, \text{Mo}}^{\dow}, \, s\in {\mc S}), \ldots, (\theta_{s, \text{Su}}^{\dow}, \, s\in {\mc S}), \, (\theta_{\text{Mo}}^{\dow}, \ldots, \theta_{\text{Su}}^{\dow}) \right).
\end{align*}

By re-parametrizing the fused lasso penalty terms, we re-write (<ref>) as,

\begin{align}
\min_{\Theta, \Theta_H, \Theta_D, \Delta, \Gamma, \Psi} & \sum_{i=1}^{n} \mu_i(\Theta, \Theta_H, \Theta_D, \Delta) - \sum_{i=1}^{n} y_i \log(\mu_i(\Theta, \Theta_H, \Theta_D, \Delta)) \nn \\
&+ \lambda
\sum_{s \in \mc S} \l( \sum_{h \in \mc H} \vert \theta^{\hod}_{s, h} \vert
+ \sum_{d \in \mc D} \vert \theta^{\dow}_{s, d} \vert \r)
\nn \\
&+ \lambda_N \cdot \widetilde{P}_{N} (\Gamma) + \widetilde{P}_H(\Psi)
+ \mbf I_{\mathcal{C}}\left( \Theta, \Theta_{H}, \Theta_{D}, \Gamma, \Psi \right).
\label{admm:obj2}
\end{align}

Here, the penalty functions are defines as

\begin{align*}
&\widetilde{P}_N(\Gamma) =
\sum_{s \in \mc S} \sqrt{ \vert \mc N_r(s) \vert \, \sum_{s^\prime \in \mc N_r(s)} \l[
\left(\gamma_{s, s'}\right)^2 + \sum_{h=1}^{23} \left(\gamma_{s, s',h}^{\hod}\right)^2 + \sum_{d \in \mc D \setminus \{\text{Mo}\}} \left(\gamma_{s, s',d}^{\dow}\right)^2
\r]},
\\
&\text{where \ }
\Gamma = (\bm\gamma_s, \, s \in \mc S) \text{ \ with \ }
\bm\gamma_s = \l( \gamma_{s, s'}, (\gamma^{\hod}_{s, s', h}, \, h \in \mc H \setminus \{0\}),
(\gamma^{\dow}_{s, s', d}, \, d \in \mc D \setminus \{\text{Mo}\})), \, s' \in \mc N_r(s) \r),
% \left(\left(\begin{array}{c}
%   \left( \gamma_{s, s'}, \, s' \in \mc N_r(s) \right)^\top
%       \\
%      \left( (\gamma^{\hod}_{s, s', 1}, \ldots, \gamma^{\hod}_{s, s',23}), \, s' \in \mc N_r(s) \right)^\top
%             \\
%      \left((\gamma^{\dow}_{s, s',\text{Tu}}, \ldots, \gamma^{\dow}_{s, s', \text{Su}}), %\,  s' \in \mc N_r(s))\right)^\top
%\right), \; s \in {\mc S} \right)
\end{align*}

\begin{align*}
\widetilde{P}_H(\Psi) &= \sum_{s \in \mc S}\sum_{h \in \{0, \ldots, 23 \}} \vert \psi_{s,h} \vert
\quad \text{with} \quad
\Psi = \l( \bm\psi_s = (\psi_{s, 0}, \ldots, \psi_{s, 23}), \, s \in \mc S \r).
\end{align*}

In addition, $I_{\mathcal{C}}$ denotes an indicator function on the constraint set $\mathcal{C}$ which forces the newly defined penalty functions $\widetilde{P}_{N}(\Gamma)$ (resp.  $\widetilde{P}_H(\Psi)$) to be the same as ${P}_N(\Theta, \Theta_{H}, \Theta_{D})$ (resp. $P_H(\Theta, \Theta_{H}, \Theta_{D})$) such that

\begin{align*}
I_\mathcal{C}(\mbf z) =
\begin{cases}
0 & \text{ if  } \mbf z \in \mathcal{C} \\
\infty & \text{ otherwise}.
\end{cases}
\end{align*}

The constraint set $\mathcal{C}$ is a subset of $\vert (\Theta, \Theta_H, \Theta_D, \Gamma, \Psi) \vert$-dimensional space whose elements fulfil

\begin{align}
\gamma_{s, s'} &= \sqrt{2} (\theta_{s} - \theta_{s'}),  \nn \\
\gamma_{s, s',h}^{\hod} &= \theta_{s} - \theta_{s'} + \theta_{s,h}^{\hod} - \theta_{s', h}^{\hod} \: \text{ for } h \in \mathcal{H}\setminus \{0\}, \nn  \\
\gamma_{s, s',d}^{\dow} & =  \theta_{s} - \theta_{s'} + \theta_{s,d}^{\dow} - \theta_{s', d}^{\dow} \: \text{ for } d \in \mathcal{D}\setminus \{\text{Mo}\}
\label{admm:para1}
\end{align}

for all $s, s' \in \mc S$, and

\begin{align}
\label{admm:para2}
\psi_{s,h} &= \theta_{h+1}^{\hod}-\theta_{h}^{\hod} + \theta_{s,h+1}^{\hod} - \theta_{s,h}^{\hod} \: \text{ for } h \in \{1, \ldots 22\}, \nn \\
\psi_{s, 0} &= \theta_{1}^{\hod} + \theta_{s, 1}^{\hod}
\quad \text{and} \quad
\psi_{s, 23} = -\theta_{23}^{\hod} -\theta_{s, 23}^{\hod},
\end{align}

for all $s \in \mc S$,
where (<ref>) and (<ref>) account for that the baseline parameters are set to be zero.

To utilize ADMM, we rewrite the objective function (<ref>) as follows:

\begin{align}
\label{admm:obj3}
\min_{ \substack {\Theta, \Theta_H, \Theta_D, \Delta, \Gamma, \Psi, \\
Z_{\Theta}, Z_{\Theta_{H}}, Z_{\Theta_{D}}, S_{\Gamma}, S_{\Psi}
}} & \sum_{i=1}^{n} \mu_i(\Theta, \Theta_H, \Theta_D, \Delta) - \sum_{i=1}^{n} y_i \log(\mu_i((\Theta, \Theta_H, \Theta_D, \Delta))) \nn \\
&+ \lambda
\sum_{s \in \mc S} \l( \sum_{h \in \mc H} \vert \theta^{\hod}_{s, h} \vert
+ \sum_{d \in \mc D} \vert \theta^{\dow}_{s, d} \vert \r)
\nn \\
&+ \lambda_N \cdot \widetilde{P}_{N} (\Gamma) + \widetilde{P}_H(\Psi) + \mbf I_{\mathcal{C}}\left( Z_{\Theta}, Z_{\Theta_{H}}, Z_{\Theta_{D}}, S_{\Gamma}, S_{\Psi} \right)
\\
\text{subject to} &  \quad (\Theta, \Theta_{H}, \Theta_{D}) = (Z_{\Theta}, Z_{\Theta_H}, Z_{\Theta_D})
\text{ and } (\Gamma, \Psi) = (S_{\Gamma}, S_{\Psi}). \nn
\end{align}

The ADMM optimizes (<ref>) in three steps as follows:

* Update the primal variables as

\begin{align}
(\Theta^{k+1}, \Theta_{H}^{k+1}, \Theta_{D}^{k+1}, \Delta^{k+1})  =& \argmin_{\Theta, \Theta_H, \Theta_D, \Delta}
\sum_{i = 1}^{n} \mu_i(\Theta, \Theta_H, \Theta_D, \Delta) - \sum_{i=1}^{n} y_i \log(\mu_i((\Theta, \Theta_H, \Theta_D, \Delta))) \nn \\
&+ \lambda
\sum_{s \in \mc S} \l( \sum_{h \in \mc H} \vert \theta^{\hod}_{s, h} \vert
+ \sum_{d \in \mc D} \vert \theta^{\dow}_{s, d} \vert \r)
\nn \\
&+ \frac{\rho}{2} \l\Vert   (\Theta, \Theta_{H}, \Theta_{D}) -   (Z_{\Theta}^k, Z_{\Theta_{H}}^k, Z_{\Theta_{D}}^k) +  (U_{\Theta}^k, U_{\Theta_{H}}^k, U_{\Theta_{D}}^k) \r\Vert^2,
\label{admm:step:one:one} \\
\Gamma^{k+1} =& \, \argmin_{\Gamma} \lambda_N \widetilde{P}_N(\Gamma) + \frac{\rho}{2} \l\Vert \Gamma - S_{\Gamma}^{k} + T_{\Gamma}^{k} \r\Vert^{2},
\label{admm:step:one:two} \\
\Psi^{k+1} =& \, \argmin_{\Psi} \lambda_H \widetilde{P}_N(\Psi) + \frac{\rho}{2} \l\Vert \Psi - S_{\Psi}^{k} + T_{\Psi}^{k} \r\Vert^{2},
\label{admm:step:one:three}
\end{align}

where $(U_{\Theta}, U_{\Theta_H}, U_{\Theta_D})$, $T_{\Gamma}$, and $T_{\Psi}$ are dual variables associated with the constraints $(\Theta, \Theta_H, \Theta_D)=(Z_{\Theta}, Z_{\Theta_H}, Z_{\Theta_D})$, $\Gamma = S_{\Gamma}$ and $\Psi = S_{\Psi}$, respectively.

* Update $( Z_{\Theta}, Z_{\Theta_{H}}, Z_{\Theta_{D}}, S_{\Gamma}, S_{\Psi})$ by projecting $(\Theta^{k+1} + U_{\Theta}^{k}, \Theta_H^{k+1} + U_{\Theta_H}^{k}, \Theta_{D}^{k+1} + U_{\Theta_D}^{k}, \Gamma^{k+1} + T_{\Gamma}^{k}, \Psi^{k+1} + T_{\Psi}^{k})$ onto the constraint set $\mathcal{C}$, as

\begin{align}
& \left( Z_{\Theta}^{k+1}, Z_{\Theta_{H}}^{k+1}, Z_{\Theta_{D}}^{k+1}, S_{\Gamma}^{k+1}, S_{\Psi}^{k+1} \right)
\nn \\
=&  \prod_{\mathcal{C}} (\Theta^{k+1}+U_{\Theta}^{k}, \Theta_H^{k+1}+U_{\Theta_H}^{k}, \Theta_{D}^{k+1}+U_{\Theta_D}^{k}, \Gamma^{k+1}+T_{\Gamma}^{k}, \Psi^{k+1}+T_{\Psi}^{k}),
\label{admm:step:two}
\end{align}

with $\prod_{\mc C}$ denoting the projection operator.

* Update the dual variable as

\begin{align*}
(U_{\Theta}^{k+1},\; U_{\Theta_H}^{k+1}, \; U_{\Theta_D}^{k+1}) =& \, (U_{\Theta}^k, U_{\Theta_H}^k, U_{\Theta_D}^k) + (\Theta^{k + 1}, \Theta_H^{k + 1}, \Theta_D^{k + 1}) -
(Z_{\Theta}^{k+1}, Z_{\Theta_H}^{k+1}, Z_{\Theta_D}^{k+1}), \\
(T_{\Gamma}^{k+1}, T_{\Psi}^{k+1}) =& \,  (T_{\Gamma}^k, T_{\Psi}^k) +  (\Gamma^{k+1}, \Psi^{k+1})
-  (S_\Gamma^{k+1}, S_\Psi^{k+1}).
\end{align*}

While Step 3 is straightforward, Steps 1 and 2 involve relatively heavy computations. The detail of Step 1 and Step 2 are discussed in the following subsections.

§.§.§ Computational details of Step 1

Let ${\mc P} = (\Theta, \Theta_{H}, \Theta_{D})$.

* We update $(\Theta, \Theta_{H}, \Theta_{D}, \Delta)$ by minimizing the objective function in (<ref>).

It in turn can be optimized

via Iteratively Reweighted Least Square (IRLS) method with the Lasso penalty as below,
at some given fixed values for $(Z_{\Theta}^k, Z_{\Theta_H}^k, Z_{\Theta_D}^k)$ and
$(U_{\Theta}^k, U_{\Theta_H}^k, U_{\Theta_D}^k)$. The $j+1$th iteration of IRLS is as follows:

\begin{align}
\label{admm:irls1}
& (\Theta^{j + 1}, \Theta_{H}^{j + 1}, \Theta_{D}^{j + 1}, \Delta^{j + 1}) =  \argmin_{\Theta, \Theta_{H}, \Theta_{D}, \Delta}
\left(\mbf X
\begin{pmatrix}
\Theta^\top \\
\Theta_{H}^\top \\
\Theta_{D}^\top \\
\Delta^\top
\end{pmatrix}
- {\mbf z}^j \right)^\top
{\mbf W}^j
\left(\mbf X
\begin{pmatrix}
\Theta^\top \\
\Theta_{H}^\top \\
\Theta_{D}^\top \\
\Delta^\top
\end{pmatrix}
- {\mbf z}^j \right) \nn \\
& + \lambda
\sum_{s \in \mc S} \l( \sum_{h \in \mc H} \vert \theta^{\hod}_{s, h} \vert
+ \sum_{d \in \mc D} \vert \theta^{\dow}_{s, d} \vert \r)
+ \frac{\rho}{2} \l\Vert (\Theta, \Theta_{H}, \Theta_{D}) -   (Z_{\Theta}^{k}, Z_{\Theta_{H}}^{k}, Z_{\Theta_{D}}^{k}) +  (U_{\Theta}^{k}, U_{\Theta_{H}}^{k}, U_{\Theta_{D}}^{k}) \r\Vert^2,
\end{align}

where $\mbf X$ is a data matrix of dimension $n \times (\vert\Delta\vert + (\vert \mc D \vert + \vert \mc H\vert -1)\cdot\vert\mc S\vert)$.
Also, ${\mbf W}^j$ is an $n \times n$ diagonal matrix with its $i$th entry being $\wh\mu_i(\Theta^j, {\Theta_H}^j, {\Theta_D}^j, \Delta^j)$, the fitted value of the $i$th observation after the $j$th iteration, and ${\mbf z}^{j}$ is a length-$n$ vector defined as follows:

\begin{align*}
{\mbf z}^{j} = {\mbf X}
(\Theta^j, \Theta_{H}^j, \Theta_{D}^j, \Delta^j)^\top
+ \left({\mbf W}^{j}\right)^{-1}{\mbf y} - \mathbbm {1}_{n}.
\end{align*}

Here, $\mbf y$ denotes a length-$n$ vector, the $i$th entry of which is the $i$th response observation, and $\mathbbm{1}_n$ denots an all-one vector of length $n$. The objective function in (<ref>) can further be written as the $\ell_1$-penalized least squares estimation problem as follows:

\begin{align}
& (\Theta^{j+1}, \Theta^{j+1}_{H}, \Theta^{j+1}_{D}, \Delta^{j+1}) =
\nn \\
& \argmin_{\Theta, \Theta_{H}, \Theta_{D}, \Delta}
\left\Vert {\mbf X}_{\text{ext}}^{j}
(\Theta, \Theta_H, \Theta_D, \Delta)^\top
- {\mbf z}_{\text{ext}}^j\right\Vert^{2}  + \lambda
\sum_{s \in \mc S} \l( \sum_{h \in \mc H} \vert \theta^{\hod}_{s, h} \vert
+ \sum_{d \in \mc D} \vert \theta^{\dow}_{s, d} \vert \r),
\label{admm:lassoform}
\end{align}

where ${\mbf X}_{\text{ext}}^{j}$ and ${\mbf z}_{\text{ext}}^{j}$ are a matrix of size
$(n + \vert \mc P \vert) \times (\vert \mc P \vert + \vert\Delta\vert)$ and a vector of length $ (n + \vert \mc P \vert)$, respectively, such that

\begin{align*}
\mbf{X}_{\text{ext}}^{j}=&
\left(
\begin{array}{  c  c }
\multicolumn{2}{c}{ \left({\mbf W}^{j}\right)^{1/2} {\mbf X} } \\
\cdashline{1-2}
\multicolumn{1}{ p{10mm} : }{$\sqrt{\rho} I_{\vert \mc P \vert} $}  & \multicolumn{1}{ p{10mm} }{$\mathbf{0}_{\vert \mc P\vert \times \vert\Delta\vert} $}   \\
\end{array}
\right), \\
{\mbf z}^{j}_{\text{ext}}=&
\left(
\begin{array}{c}
\left({\mbf W}^{j}\right)^{1/2} {\mbf z}^{j} \\
\sqrt{\rho} (Z_{\Theta}^k, Z_{\Theta_{H}}^k, Z_{\Theta_{D}}^k)^\top - \sqrt{\rho} (U_{\Theta}^k, U_{\Theta_{H}}^k, U_{\Theta_{D}}^k)^\top
\end{array}
\right).
\end{align*}

We evaluated (<ref>) using the R package

* For (<ref>)–(<ref>),
we obtain $\Gamma^{k + 1}$ and $\Psi^{k + 1}$ using a soft-thershold operator $\mathbb{S}_{\lambda}$ that takes an input vector and outputs $\mathbb{S}_{\lambda}(\mbf v) = \left(1 - \lambda/\Vert \mbf v \Vert \right)_{+} \cdot \mbf v$ with $\mathbb{S}_{\lambda}(\mathbf 0) = \mbf 0$ and $c_{+} = \max\{0, c\}$.

\begin{align*}
\bm\gamma_s^{k+1} &= \mathbb{S}_{\rho^{-1} \sqrt{{\mc N}_r(s)}\lambda_{N}}  \l(S_{\bm\gamma_s}^{k} - T_{\bm\gamma_s}^{k}\r) \text{ \ for each \ } s \in {\mc S}, \\
\psi_{s, h}^{k+1} &= \mathbb{S}_{\rho^{-1}\lambda_{H}}(S_{\Psi, s, h}^k - T_{\Psi, s, h}^k)   \: \text{ for } \: s \in {\mc S} \text{, } \: h \in {\mc H}.
\end{align*}

§.§.§ Computational details of Step 2

In Step 2, the update of $\left( Z_{\Theta}, Z_{\Theta_{D}}, Z_{\Theta_{H}}, S_{\Gamma}, S_{\Psi} \right)$ is achieved via projection in (<ref>),
where ${\mc C}$ is the constraint set specified in (<ref>) and (<ref>). This step is the bottleneck of the computation due to the large number of variables to be updated. In this section, the details of the procedure is illustrated.
Throughout, we denote by $\mbf I$ and $\mbf 0$ an identity matrix and
a matrix of zeros, respectively, and their dimensions are determined by the context unless specified.

The projection $(A_1, A_2, A_3, B_1, B_2)= \prod_{\mc C} (E_1, E_2, E_3, F_1, F_2)$ is equivalent to the following minimization problem

\begin{align*}
\min_{A_1,A_2,A_3,B_1, B_2} & \Vert A_1 -E_1\Vert^2+\Vert A_2 -E_2\Vert^2 +\Vert A_3 -E_3\Vert^2+\Vert B_1 -F_1\Vert^2+\Vert B_2 -F_2\Vert^2 \\
\text{subject to} & \quad B_1 = (A_1, A_2, A_3) \mbf D_\Theta^\top \text{ and }
B_2 = (A_1, A_2, A_3) \mbf D_{\mc H}^\top,
\end{align*}

where $\mbf D_\Theta$ and $\mbf D_{\mc H}$ are matrices encoding the constraints (<ref>) and  (<ref>), respectively. Then, the above optimization problem can be re-written as

\begin{align*}
\min_{A_1,A_2,A_3} & \Vert A_1 -E_1\Vert^2+\Vert A_2 -E_2\Vert^2 +\Vert A_3 -E_3\Vert^2+\Vert (A_1,A_2,A_3) \mbf D_\Theta^\top - F_1\Vert^2+\Vert (A_1, A_2, A_3) \mbf D_{\mc H}^\top - F_2 \Vert^2,
\end{align*}

and its optimizer is the solution of its normal equation

\begin{align}
\label{projsol}
(A_1, A_2, A_3) \underbrace{\left( \mbf I + \mbf D_\Theta^\top \mbf D_\Theta + \mbf D_{\mc H}^\top \mbf D_{\mc H}
%\left( \begin{array}{c}
%     \mathbf{0} \\
%      \cdashline{1-1}
%      \mbf D_{\mc H}^\top \\
%      \cdashline{1-1}
%     \mbf 0
%\end{array} \right) \left( \begin{array}{c}
%     \mathbf{0} \\
%      \cdashline{1-1}
%      \mbf D_{\mc H}^\top \\
%      \cdashline{1-1}
%     \mbf 0
%\end{array} \right)^\top
\right)}_{\mbf P}
= (E_1, E_2, E_3) + F_1 \mbf D_\Theta + F_2 \mbf D_{\mc H}.
\end{align}

Once the inverse of $\mbf P$ is available,
the solution $(A_1, A_2, A_3)$ of (<ref>) can be calculated in a straightforward manner. Also, the inverse matrix remains the same throughout the iterations and thus no re-computation is required.
In our problem, however, as the size of the matrix $\mbf P$ is huge, its dimension reaching approximately $47,000 \times 47,000$, and inverting this matrix can be very demanding with the computational complexity of $\mathcal{O}(10^{11})$. Additionally, even if we can compute the inverse matrix, it is huge in size and occupies a large portion of memory space which hinders efficient computation.
Given the situation, we avoid direct computation of $\mbf P^{-1}$ and find the solution of (<ref>) by utilizing the specific structure of the matrix $\mbf P$.

We start by defining an $M \times \vert\mc S\vert$ matrix $\dnet$ to be a matrix that represents the network constructed from the neighborhood relations so that each row is associated with two connected stations where $M = \sum_{s \in {\mc S}} \vert{\mc N}_r(s)\vert$. Defining $\mc R_{\text{ind}}: \mc S \times \mc S \rightarrow \{1, \ldots, M\}$ to be the mapping that returns the row index of $s-s'
$ connection for $s' \in {\mc N}_r(s)$, the $\mc R_{\text{ind}}(s, s')$-th row of $\dnet$ is given by $\mbf e_s - \mbf e_{s'}$,
where $\mbf e_i$ is a standard basis vector of length $\vert\mc S\vert$.
That is, each row of $\dnet$ is composed of $\{-1,0,1\}$ with exactly one $1$ and one $-1$ and the rest of the entries are all $0$s. By construction, the Laplacian matrix, say $\lnet$, of the network can be represented by $\dnet$ as follows:

\begin{align*}
\lnet = \frac{1}{2}\dnet^\top \dnet.
\end{align*}

Then, we have

\begin{align*}
\dnet \Theta^\top &= (\theta_s - \theta_{s'}, \,  s' \in \mc N_r(s), \,  s \in {\mc S})^\top, \\
\dnet (\Theta^{\circ}_{H, h})^\top &= (\theta^{\hod}_{s, h} - \theta^{\hod}_{s', h}, \, s' \in \mc N_r(s), \, s \in \mc S)^\top, \\
% \left(\begin{array}{c}
% (\theta_{s,1}^{\hod} - \theta_{s',1}^{\hod} \vert \; (s,s') \text{ with } s \in % {\mc S} \;\& \; {\mc N}_r(s) ) \\
% \cdashline{1-1}
% \vdots \\
% \cdashline{1-1}
% (\theta_{s,23}^{\hod} - \theta_{s',23}^{\hod} \vert \; (s,s') \text{ with } s \in %{\mc S} \;\& \; {\mc N}_r(s) )
%\end{array} \right) \\
\dnet (\Theta^{\circ}_{D, d})^\top &= (\theta^{\dow}_{s, d} - \theta^{\hod}_{s', d}, \, s' \in \mc N_r(s), \, s \in \mc S)^\top, \\
% \left(
% \begin{array}{c}
% (\theta_{s,\text{Tu}}^{\dow} - \theta_{s',\text{Tu}}^{\dow} \vert \;  (s,s') \text{ with } s \in {\mc S} \;\& \; {\mc N}_r(s) )       \\
% \cdashline{1-1}
% \vdots \\
% \cdashline{1-1}
%  (\theta_{s,\text{Su}}^{\dow} - \theta_{s',\text{Su}}^{\dow} \vert\; (s,s') \text{ with } s \in {\mc S} \;\& \; {\mc N}_r(s) )
\end{align*}

where $\Theta_{H, h}^{\circ} = (\theta^{\hod}_{s, h}, \, s \in \mc S)$ and $\Theta_{D, d}^{\circ} = (\theta^{\dow}_{s, d}, \, d \in \mc D)$.

The matrix $\mbf D_\Theta$ can be written with $\dnet$ as

[matrix of math nodes,left delimiter=(,right delimiter=),row sep=0.5cm,column sep=0.5cm] (m)
0      0       0      -     0     0     0      0

0     ⋱     0      0     ⋮     0       0      0      0

0      0          0    -     0      0      0      0

0      0      0      0     -          0      0      0

0      0      0      0     ⋮     0      ⋱     0      0

0      0      0      0    -     0     0           0

0      0      0      0    -√(2)     0      0      0      0

[dashed] ($0.5*(m-1-4.north east)+0.5*(m-1-3.north west)$) –
($0.5*(m-7-4.south east)+0.5*(m-7-3.south west)$);

[dashed] ($0.5*(m-1-4.north east)+0.5*(m-1-5.north west)$) –
($0.3*(m-7-4.south east)+0.7*(m-7-5.south west)$);

[dashed] ($0.5*(m-1-5.north east)+0.5*(m-1-6.north west)$) –
($0.8*(m-7-5.south east)+.2*(m-7-6.south west)$);
[dashed] ($0.5*(m-1-8.north east)+0.5*(m-1-9.north west)$) –
($0.5*(m-7-8.south east)+0.5*(m-7-9.south west)$);

[dashed] ($0.5*(m-3-1.south west)+0.5*(m-4-1.north west)$) –
($0.5*(m-3-9.south east)+0.5*(m-4-9.north east)$);
[dashed] ($0.5*(m-6-1.south west)+0.5*(m-7-1.north west)$) –
($0.5*(m-6-9.south east)+0.5*(m-7-9.north east)$);

[above=3pt of m-1-1] (top-1) ;
[above=3pt of m-1-2] (top-2) ;
[above=3pt of m-1-3] (top-3) ;
[above=3pt of m-1-4] (top-4) ;
[above=3pt of m-1-5] (top-5) ;
[above=3pt of m-1-6] (top-6) ;
[above=3pt of m-1-7] (top-7) ;
[above=3pt of m-1-8] (top-8) ;
[above=3pt of m-1-9] (top-9) ;

[left=13pt of m-1-1] (left-1) ;
[left=13pt of m-2-1] (left-2) ;
[left=13pt of m-3-1] (left-3) ;
[left=13pt of m-4-1] (left-4) ;
[left=13pt of m-5-1] (left-5) ;
[left=13pt of m-6-1] (left-6) ;
[left=17pt of m-7-1] (left-7) ;

[rectangle,above delimiter={] (del-top-1) at ($0.5*(top-1.south) +0.5*(top-3.south)$) (0,0) rectangle (3.5,0);;
[above=10pt] at (del-top-1.north) $\scard \cdot \dcard$;

[rectangle,above delimiter={] (del-top-2) at ($1.0*(top-4.south)$) (0,0) rectangle (0.7,0);;
[above=10pt] at (del-top-2.north) $\dcard$;

[rectangle,above delimiter={] (del-top-3) at ($1.0*(top-5.south)$) (0,0) rectangle (1.6,0);;
[above=10pt] at (del-top-3.north) $\scard$;

[rectangle,above delimiter={] (del-top-4) at ($1.0*(top-7.south)$) (0,0) rectangle (3.5,0);;
[above=10pt] at (del-top-4.north) $\scard \cdot \hcard $;

[rectangle,above delimiter={] (del-top-5) at ($1.0*(top-9.south)-(0.15,0)$) (0,0) rectangle (0.6,0);;
[above=10pt] at (del-top-5.north) $\hcard$;

[rectangle,left delimiter={] (del-left-1) at ($0.5*(left-1.east) +0.5*(left-3.east)$) (left-1.north east) rectangle (left-3.south west);;
[left=10pt] at (del-left-1.west) $M \cdot \dcard$;

[rectangle,left delimiter={] (del-left-2) at ($0.5*(left-4.east) +0.5*(left-6.east)$) (left-4.north east) rectangle (left-6.south west);;
[left=10pt] at (del-left-2.west) $M \cdot \hcard$;

[rectangle,left delimiter={] (del-left-2) at ($1.0*(left-7.east)$) (0,.5) rectangle (0,1.1);;
[left=10pt] at (del-left-2.west) $M$;

where ${\mc D}^{\circ}={\mc D}\setminus\{\text{Mo}\}$ and ${\mc H}^{\circ}={\mc H}\setminus\{0\}$ so that the baseline parameters are removed.

Denoting $\frac{1}{2}\lnet$ by $\wtlnet$,  the form of $\mbf D_\Theta^\top \mbf D_\Theta$ is as below:

[matrix of math nodes,left delimiter=(,right delimiter=),row sep=0.5cm,column sep=0.5cm] (m)
0      0       0      -     0     0     0      0

0     ⋱     0      0     ⋮     0       0      0      0

0      0          0    -     0      0      0      0

0      0     0      0     0      0      0      0      0

-     ⋯     -       0     (+)    -     ⋯     -     0

0      0      0      0     -          0      0      0

0      0      0      0     ⋮     0      ⋱     0      0

0      0      0      0    -     0     0           0

0      0      0      0     0      0      0      0      0

[dashed] ($0.5*(m-1-4.north east)+0.5*(m-1-3.north west)$) –
($0.5*(m-9-4.south east)+0.5*(m-9-3.south west)$);

[dashed] ($0.6*(m-1-4.north east)+0.4*(m-1-5.north west)$) –
($0.7*(m-9-4.south east)+0.3*(m-9-5.south west)$);

[dashed] ($0.5*(m-1-5.north east)+0.5*(m-1-6.north west)$) –
($0.4*(m-9-5.south east)+.6*(m-9-6.south west)$);

[dashed] ($0.5*(m-1-8.north east)+0.5*(m-1-9.north west)$) –
($0.5*(m-9-8.south east)+0.5*(m-9-9.south west)$);

[dashed] ($0.5*(m-3-1.south west)+0.5*(m-4-1.north west)$) –
($0.5*(m-3-9.south east)+0.5*(m-4-9.north east)$);
[dashed] ($0.5*(m-4-1.south west)+0.5*(m-5-1.north west)$) –
($0.5*(m-4-9.south east)+0.5*(m-5-9.north east)$);
[dashed] ($0.5*(m-5-1.south west)+0.5*(m-6-1.north west)$) –
($0.5*(m-5-9.south east)+0.5*(m-6-9.north east)$);
[dashed] ($0.5*(m-8-1.south west)+0.5*(m-9-1.north west)$) –
($0.5*(m-8-9.south east)+0.5*(m-9-9.north east)$);

[above=3pt of m-1-1] (top-1) ;
[above=3pt of m-1-2] (top-2) ;
[above=3pt of m-1-3] (top-3) ;
[above=3pt of m-1-4] (top-4) ;
[above=3pt of m-1-5] (top-5) ;
[above=3pt of m-1-6] (top-6) ;
[above=3pt of m-1-7] (top-7) ;
[above=3pt of m-1-8] (top-8) ;
[above=3pt of m-1-9] (top-9) ;

[rectangle,above delimiter={] (del-top-1) at ($0.5*(top-1.south) +0.5*(top-3.south)$) (0,0) rectangle (3.5,0);;
[above=10pt] at (del-top-1.north) $\scard \cdot \dcard$;

[rectangle,above delimiter={] (del-top-2) at ($1.0*(top-4.south)$) (0,0) rectangle (0.7,0);;
[above=10pt] at (del-top-2.north) $\dcard$;

[rectangle,above delimiter={] (del-top-3) at ($1.0*(top-5.south)$) (0,0) rectangle (1.6,0);;
[above=10pt] at (del-top-3.north) $\scard$;

[rectangle,above delimiter={] (del-top-4) at ($1.0*(top-7.south)$) (0,0) rectangle (3.5,0);;
[above=10pt] at (del-top-4.north) $\scard \cdot \hcard $;

[rectangle,above delimiter={] (del-top-5) at ($1.0*(top-9.south)-(0.15,0)$) (0,0) rectangle (0.6,0);;
[above=10pt] at (del-top-5.north) $\hcard$;

Now, we define $\mbf D_{\mc H}$, the difference matrix that addresses the association between consecutive hours of each of the stations so that

\begin{align*}
\mbf D_{\mc H} (\Theta, \Theta_D, \Theta_{H})^\top = \left(
\begin{array}{c}
(- (\theta_{1}^{\hod} + \theta_{s,1}^{\hod}), \,  s \in {\mc S} )^\top\\
\cdashline{1-1}
((\theta_{1}^{\hod} + \theta_{s,1}^{\hod})- (\theta_{2}^{\hod} + \theta_{s,2}^{\hod}), \,  s \in {\mc S} )^\top \\
\vdots \\
((\theta_{22}^{\hod} + \theta_{s,22}^{\hod})- (\theta_{23}^{\hod} + \theta_{s,23}^{\hod}), \,  s \in {\mc S} )^\top \\
\cdashline{1-1}
((\theta_{23}^{\hod} + \theta_{s,23}^{\hod})^\top, \,  s \in {\mc S} )^\top
\end{array}
\right).
\end{align*}

The specific form of $\mbf D_{\mc H}$ is as follows:

[matrix of math nodes,left delimiter=(,right delimiter=),row sep=0.5cm,column sep=0.4cm] (m)
0     0     0     0      0       0      0      0    0     -1      0     0     0

0     0      0     0      0       0      0      0     0     1     -1     0     0

0     0    0     0      0       0      0      0     0       0     ⋱    ⋱    0

0     0     0     0      0       0      0      0     0     0     0     1     -1

0     0     0     0      0       0      0      0     0    0     0     0     1

0     0     0     0      0       -I_      0      0     0     -1_      0     0     0

0     0     0     0      0       I_      -I_      0     0     1_     -1_     0     0

0     0     0     0      0      0      ⋱    ⋱     0     0     ⋱    ⋱    0

0     0     0     0      0      0      0     I_      -I_     0     0    1_     -1_

0     0     0     0      0      0      0     0      I_     0    0     0     1_

[dashed] ($0.5*(m-1-4.north east)+0.5*(m-1-3.north west)$) –
($0.5*(m-10-4.south east)+0.5*(m-10-3.south west)$);
[dashed] ($0.5*(m-1-5.north east)+0.5*(m-1-4.north west)$) –
($0.5*(m-10-5.south east)+0.5*(m-10-4.south west)$);
[dashed] ($0.5*(m-1-6.north east)+0.5*(m-1-5.north west)$) –
($0.5*(m-10-6.south east)+0.5*(m-10-5.south west)$);

[dashed] ($0.5*(m-1-10.north east)+0.5*(m-1-9.north west)$) –
($0.6*(m-10-10.south east)+0.4*(m-10-9.south west)$);

[dashed] ($0.5*(m-5-1.south west)+0.5*(m-6-1.north west)$) –
($0.5*(m-5-13.south east)+0.5*(m-6-13.north east)$);

[above=3pt of m-1-1] (top-1) ;
[above=3pt of m-1-2] (top-2) ;
[above=3pt of m-1-3] (top-3) ;
[above=3pt of m-1-4] (top-4) ;
[above=3pt of m-1-5] (top-5) ;
[above=3pt of m-1-6] (top-6) ;
[above=3pt of m-1-7] (top-7) ;
[above=3pt of m-1-8] (top-8) ;
[above=3pt of m-1-9] (top-9) ;
[above=3pt of m-1-10] (top-10) ;
[above=3pt of m-1-11] (top-11) ;
[above=3pt of m-1-12] (top-12) ;
[above=3pt of m-1-13] (top-13) ;

[left=13pt of m-1-1] (left-1) ;
[left=13pt of m-2-1] (left-2) ;
[left=13pt of m-3-1] (left-3) ;
[left=13pt of m-4-1] (left-4) ;
[left=13pt of m-5-1] (left-5) ;
[left=13pt of m-6-1] (left-6) ;
[left=13pt of m-7-1] (left-7) ;
[left=13pt of m-8-1] (left-8) ;
[left=13pt of m-9-1] (left-9) ;
[left=13pt of m-10-1] (left-10) ;

[rectangle,above delimiter={] (del-top-1) at ($1.0*(top-2.south)$) (0,0) rectangle (2.2,0);;
[above=10pt] at (del-top-1.north) $\scard \cdot \dcard$;

[rectangle,above delimiter={] (del-top-2) at ($1.0*(top-4.south)$) (0,0) rectangle (0.5,0);;
[above=10pt] at (del-top-2.north) $\dcard$;

[rectangle,above delimiter={] (del-top-3) at ($1.0*(top-5.south)$) (0,0) rectangle (0.5,0);;
[above=10pt] at (del-top-3.north) $\scard$;

[rectangle,above delimiter={] (del-top-4) at ($0.5*(top-7.south)+0.5*(top-8.south)$) (0,0) rectangle (4.7,0);;
[above=10pt] at (del-top-4.north) $\scard \cdot \hcard $;

[rectangle,above delimiter={] (del-top-5) at ($0.5*(top-11.south)+0.5*(top-12.south)$) (0,0) rectangle (4.7,0);;
[above=10pt] at (del-top-5.north) $ \hcard $;

[rectangle,left delimiter={] (del-left-1) at ($0.5*(left-1.east) +0.5*(left-5.east)$) (left-6.north east) rectangle (left-10.south west);;
[left=10pt] at (del-left-1.west) $ \vert \mc H \vert$;

[rectangle,left delimiter={] (del-left-1) at ($0.5*(left-6.east) +0.5*(left-10.east)$) (left-1.north east) rectangle (left-5.south west);;
[left=10pt] at (del-left-1.west) $ \vert \mc H \vert \cdot \scard$;

Thus, $\mbf D_{\mc H}^\top \mbf D_{\mc H}$ has the form

\begin{align}
D_{\mc H}^\top D_{\mc H} =&
\left( \begin{array}{c : c : c}
{\mbf 0} & {\mbf 0} & {\mbf 0}\\
\cdashline{1-3}
{\mbf 0}^\top & \mbf M_1 & \mbf M_2 \\
{\mbf 0}^\top & \mbf M_2^\top & \mbf M_3
% {\mbf 0}_{ 1  } & {\mbf 0}_{ 2  } & {\mbf 0}_{ 3  }\\
% \cdashline{1-3}
% {\mbf 0}_{2}^\top & \mbf M_1 & \mbf M_2 \\
%  {\mbf 0}_{3}^\top & \mbf M_2^\top & \mbf M_3
\end{array} \right), \quad \text{where}
\nn \\
% where ${\mbf 0}_1$, ${\mbf 0}_2$, and ${\mbf 0}_3$ are zero-matrices of dimensions $ (\scard\cdot \dcard + \scard+ {\vert \mc D \vert})  \times ((\scard\cdot \dcard) + \scard + {\vert \mc D \vert})$, $(\scard\cdot \dcard+ \scard+ {\vert \mc D \vert}) \times \scard \cdot \hcard $, and $(\scard\cdot \dcard+ \scard + {\vert \mc D \vert}) \times \hcard$ respectively, and $\mbf M_1$, $\mbf M_2$, and $\mbf M_3$ are matrices of dimensions $\scard \cdot \hcard \times \scard \cdot \hcard $, $\scard \cdot \hcard \times \hcard $, and $\hcard  \times  \hcard $ respectively, such that
\mbf M_1=& \begin{pmatrix}
2 \mbf I_{\scard} & -\mbf I_{\scard} & {\mbf 0} & {\mbf 0} & {\mbf 0} \\
-\mbf I_{\scard} &  2 \mbf I_{\scard} & -\mbf I_{\scard} & {\mbf 0} & {\mbf 0} \\
{\mbf 0} & \ddots & \ddots & \ddots  &  {\mbf 0}\\
{\mbf 0} &  {\mbf 0} &   -\mbf I_{\scard} &  2 \mbf I_{\scard} & -\mbf I_{\scard} \\
{\mbf 0} &    {\mbf 0} &  {\mbf 0} &   -\mbf I_{\scard} &  2 \mbf I_{\scard}
\end{pmatrix}, \label{proj_def_m1} \\
\mbf M_2= & \begin{pmatrix}
2 \mathbbm{1}_{\scard} & -\mathbbm{1}_{\scard} & {\mbf 0} & {\mbf 0} & {\mbf 0} \\
-\mathbbm{1}_{\scard} &  2 \mathbbm{1}_{\scard} & -\mathbbm{1}_{\scard} & {\mbf 0} & {\mbf 0} \\
{\mbf 0} & \ddots & \ddots & \ddots  &  {\mbf 0}\\
{\mbf 0} &  {\mbf 0} &   -\mathbbm{1}_{\scard} &  2 \mathbbm{1}_{\scard} & -\mathbbm{1}_{\scard} \\
{\mbf 0} &    {\mbf 0} &  {\mbf 0} &   -\mathbbm{1}_{\scard} &  2 \mathbbm{1}_{\scard}
\end{pmatrix}, \nn \\
\mbf M_3= & \begin{pmatrix}
2 (\scard + 1) & -(\scard + 1) & 0 & 0 & 0 \\
-(\scard + 1) &  2 (\scard + 1) & -(\scard + 1) & 0 & 0 \\
0 &   \ddots & \ddots & \ddots & 0 \\
0 &  0 &   -(\scard + 1) &  2 (\scard + 1) & -(\scard + 1) \\
0 &    0 &  0 &   -(\scard + 1) &  2 (\scard + 1)
\end{pmatrix}. \nn
\end{align}

From these, we can illustrate the structure of $\mbf P$.
One notable thing is that both $\mbf D_{\Theta}$ and $\mbf D_{\mc H}$ have
their columns corresponding to $(\theta_{\text{Tu}}^{\dow}, \ldots, \theta_{\text{Su}}^{\dow})$ set exactly to be zero.
Thus, in computing $\mbf P^{-1}$, we remove their corresponding columns so that the inverse matrix is applied only to $(\Theta, \Theta_{D^{\circ}}, \Theta_{H})$. The structure of
$\bar{\mbf P}$, a sub-matrix of $\mbf P = \mbf I + \mbf D_{\Theta}^\top \mbf D_{\Theta}+ \mbf D_{\mc H}^\top \mbf D_{\mc H}$ without the columns corresponding to the daily parameters, is as follows:

\begin{align}
% \label{inverse1}
% &\mbf I + \mbf D_{\Theta}^\top \mbf D_{\Theta}+ \mbf D_{\mc H}^\top \mbf D_{\mc H} =
& \bar{\mbf P} = \nn \\
\begin{array}{ccc:c:ccccc:c}
\wtlnet + \mbf I & {\mbf 0} & {\mbf 0} & -\wtlnet & \multicolumn{5}{c:}{{\mbf 0} } & {\mbf 0}\\
{\mbf 0}  & \ddots & {\mbf 0} & \vdots & \multicolumn{5}{c:}{ \vdots } & {\mbf 0} \\
{\mbf 0}  &  {\mbf 0} &\wtlnet + \mbf I & -\wtlnet & \multicolumn{5}{c:}{{\mbf 0}} &  {\mbf 0} \\
\cdashline{1-10}
-\wtlnet & \cdots & -\wtlnet & (\vert \mc D \vert + \vert \mc H\vert)\wtlnet + \mbf I & -\wtlnet & \multicolumn{3}{c}{{\cdots }} & -\wtlnet & {\mbf 0}  \\
\cdashline{1-10}
{\mbf 0} &   \cdots &   {\mbf 0} & -\wtlnet & \wtlnet + 3\mbf I & -\mbf I &   {\mbf 0} & {\mbf 0} & {\mbf 0} & \multirow{5}{*}{$\mbf M_2$} \\
\multirow{3}{*}{$\vdots$} &  \multirow{3}{*}{$\ddots$} &  \multirow{3}{*}{$\vdots$} &  \multirow{3}{*}{$\vdots$} & -\mbf I & \wtlnet + 3\mbf I & -\mbf I &  {\mbf 0} & \\
& & &  &   {\mbf 0} & \ddots & \ddots & \ddots &  {\mbf 0} & \\
&  &  &  &  {\mbf 0} &  {\mbf 0} & -\mbf I & \wtlnet +3\mbf I & -\mbf I & \\
{\mbf 0} & \cdots & {\mbf 0} & -\wtlnet & {\mbf 0} & {\mbf 0} & {\mbf 0} & -\mbf I & \wtlnet + 3\mbf I & \\
\cdashline{1-10}
{\mbf 0} & \cdots & {\mbf 0} & {\mbf 0} & \multicolumn{5}{c:}{\mbf M_2^\top} & \mbf M_3+\mbf I_{\hcard}
\end{array}
\right) \nn
\end{align}

where $\mbf M_2$ and $\mbf M_3$ are as in (<ref>). Here, the size of each block is $(\scard \cdot \dcard)$, $\scard$, $(\scard \cdot \hcard)$, and $\hcard$ from left to right and top to bottom, which match that of $\mbf D_{\Theta}^\top \mbf D_{\Theta}$ and $\mbf D_{\mc H}^\top \mbf D_{\mc H}$ specified above after removing corresponding columns of $(\theta_{\text{Tu}}^{\dow}, \ldots, \theta_{\text{Su}}^{\dow})$.

Let $\wtlnet = \mbf E \bm\Lambda \mbf E^\top$ denote the eigenvalue decomposition of $\wtlnet$ with $\bm\Lambda$ denoting a diagonal matrix having the eigenvalues as its diagonal entries.

\begin{align}
\label{inverse2}
& {\mbf W}^\top \bar{\mbf P} {\mbf W}= \nn \\
\begin{array}{ccc:c:ccccc:c}
\bm\Lambda + \mbf I & {\mbf 0} & {\mbf 0} & -\bm\Lambda & \multicolumn{5}{c:}{{\mbf 0} } & {\mbf 0}\\
{\mbf 0}  & \ddots & {\mbf 0} & \vdots & \multicolumn{5}{c:}{ \vdots } & {\mbf 0} \\
{\mbf 0}  &  {\mbf 0} &\bm\Lambda + \mbf I & -\bm\Lambda & \multicolumn{5}{c:}{{\mbf 0}} &  {\mbf 0} \\
\cdashline{1-10}
-\bm\Lambda & \cdots & -\bm\Lambda & (\vert \mc D \vert + \vert \mc H\vert)\bm\Lambda + \mbf I & -\bm\Lambda & \multicolumn{3}{c}{{\cdots }} & -\bm\Lambda & {\mbf 0}  \\
\cdashline{1-10}
{\mbf 0} &   \cdots &   {\mbf 0} & -\bm\Lambda & \bm\Lambda + 3\mbf I & - \mbf I &   {\mbf 0} & {\mbf 0} & {\mbf 0} & \multirow{5}{*}{$\wt{\mbf M}_2$} \\
\multirow{3}{*}{$\vdots$} &  \multirow{3}{*}{$\ddots$} &  \multirow{3}{*}{$\vdots$} &  \multirow{3}{*}{$\vdots$} & - \mbf I & \bm\Lambda + 3\mbf I & - \mbf I &  {\mbf 0} &  {\mbf 0} \\
& & &  &   {\mbf 0} & \ddots & \ddots & \ddots &  {\mbf 0} & \\
&  &  &  &  {\mbf 0} &  {\mbf 0} & - \mbf I & \bm\Lambda +3\mbf I & - \mbf I & \\
{\mbf 0} & \cdots & {\mbf 0} & -\bm\Lambda & {\mbf 0} & {\mbf 0} & {\mbf 0} & - \mbf I & \bm\Lambda + 3\mbf I & \\
\cdashline{1-10}
{\mbf 0} & \cdots & {\mbf 0} & {\mbf 0} & \multicolumn{5}{c:}{\wt{\mbf M}_2^\top} & \mbf M_3 + \mbf I_{ \hcard }
\end{array}
\right)
\end{align}

where $\wt{\mbf M}_2$ is

\begin{align*}
\wt{\mbf M}_2= & \begin{pmatrix}
2 \mbf E^{(\text{r})} & -\mbf E^{(\text{r})} & {\mbf 0} & {\mbf 0} & {\mbf 0} \\
-\mbf E^{(\text{r})} &  2 \mbf E^{(\text{r})} & -\mbf E^{(\text{r})} & {\mbf 0} & {\mbf 0} \\
{\mbf 0} & \ddots & \ddots & \ddots  &  {\mbf 0}\\
{\mbf 0} &  {\mbf 0} &   -\mbf E^{(\text{r})} &  2 \mbf E^{(\text{r})} & -\mbf E^{(\text{r})} \\
{\mbf 0} &    {\mbf 0} &  {\mbf 0} &   -\mbf E^{(\text{r})} &  2 \mbf E^{(\text{r})}
\end{pmatrix} \text{ \ with \ } \mbf E^{(\text{r})} = \mbf E \mathbbm{1}_{\scard}, \text{ \ and \ }
{\mbf W} =& \left(\begin{array}{ccc :c }
\mbf E & {\mbf 0} &  {\mbf 0} & {\mbf 0} \\
{\mbf 0} & \ddots & {\mbf 0} &  \vdots \\
{\mbf 0} & {\mbf 0} & \mbf E & {\mbf 0} \\
\cdashline{1-4}
\multicolumn{3}{c:}{ {\mbf 0}} & \mbf I_{\hcard}
\end{array}
\right).
\end{align*}

Thus, the solution of $\bar{\mbf P} {\mbf z} = {\mbf b}$ can be achieved by solving

\begin{align}
\label{inverse3}
\left( {\mbf W}^\top \bar{\mbf P} {\mbf W}\right)  {\widetilde {\mbf z}} = {\wt{\mbf b}}
\end{align}

with $\wt{\mbf z} = {\mbf W}^\top {\mbf z}$ and $ {\widetilde {\mbf b}} = {\mbf W}^\top {\mbf b}$, and then finally setting $\mbf z = {\mbf W} \widetilde{\mbf z}$.
Solving (<ref>) can be achieved efficiently by taking into account the structure of ${\mbf W}$ and ${\mbf W}^\top \bar{\mbf P} {\mbf W}$.
Specifically, $\mbf W$ is block diagonal with repeated blocks which reduces the matrix multiplication complexity. The matrix ${\mbf W}^\top \bar{\mbf P} {\mbf W}$ is very sparse as described in (<ref>) with all blocks involving $\bm\Lambda$ being diagonal, which facilitates efficient computation.
Therefore, (<ref>) can be solved in an iterative manner using LU decomposition. The details of this procedure is illustrated later in this section. Having the full illustration of the structures, our suggested method for the projection step is described in Algorithm <ref>.

Projection in Step 2 of ADMM

[1]
\begin{array}{l}
\left(\Theta^{k+1}, \; \Theta_{D}^{k+1}, \; \Theta_{H}^{k+1},  \;\Gamma^{k+1}, \;\Psi^{k+1}\right), \:
\left( U_{\Theta}^{k}, \;, U_{\Theta_{D}}^{k},\; U_{\Theta_{H}}^{k}, \; T_{\Gamma}^{k}, \;, T_{\Psi}^{k}\right), \\
{\mbf D}_{\Theta},\; {\mbf D}_{\mc H}, \; {\mbf W}
\end{array}

Compute ${\mbf b}_1 \leftarrow   \left(\Theta^{k+1 }, \; \Theta_{D}^{k+1}, \; \Theta_{H}^{k+1}\right)^\top +  \left( U_{\Theta}^{k}, \;, U_{\Theta_{D}}^{k},\; U_{\Theta_{H}}^{k}\right)^\top$

Compute $ {\mbf b}_2 \leftarrow \mbf D_{\Theta} \left( \Gamma^{k+1} + T_{\Gamma}^{k}\right)^\top + \mbf D_{\mc H} \left( \Psi^{k+1} + T_{\Psi}^{k}\right)^\top$

Compute $\widetilde {\mbf b}  \leftarrow {\mbf W}^\top \left({\mbf b}_1 + {\mbf b}_2 \right) $

Solve   $\left( {\mbf W}^\top\left( \mbf I + \mbf D_{\Theta}^\top \mbf D_{\Theta}+ \mbf D_{\mc H}^\top \mbf D_{\mc H} \right){\mbf W}\right)^{-1}  {\widetilde {\mbf z}} = {\widetilde {\mbf b}}$ for $\widetilde {\mbf z}$

Compute ${\mbf z} \leftarrow  {\mbf W} \widetilde{{\mbf z}}$

Set $(Z_{\Theta}^{k+1},  Z_{\Theta_{D}}^{k+1},  Z_{\Theta_{H}}^{k+1}) \leftarrow ({\mbf z}_{\Theta}, {\mbf z}_{\Theta_{D}}, {\mbf z}_{\Theta_{H}})$ where ${\mbf z} = ({\mbf z}_{\Theta}, {\mbf z}_{\Theta_{D}}, {\mbf z}_{\Theta_H})$

Set $(S_{\Gamma}^{k+1}, S_{\Psi}^{k+1}) \leftarrow (\mbf D_{\Theta} {\mbf z}, \; \mbf D_{\mc H} {\mbf z}_{\Theta_H} )$

$(Z_{\Theta}^{k+1},  Z_{\Theta_{D}}^{k+1},  Z_{\Theta_{H}}^{k+1}, S_{\Gamma}^{k+1}, S_{\Psi}^{k+1})$

In order to solve (<ref>) in line 5 of Algorithm <ref>, we can utilize the LU decomposition of the matrix in (<ref>). Writing the LU decomposition of the matrix (<ref>) by $\mbf L\mbf L^\top$, the matrix $\mbf L$ has the banded block diagonal structure as follows:

\begin{align}
\label{Lmatrix}
\mbf L = \left(\begin{array}{c :c :c :c }
\mbf A_{11} & {\mbf 0} &{\mbf 0} & {\mbf 0}\\
\cdashline{1-4}
\mbf A_{21} & \mbf A_{22}& {\mbf 0} & {\mbf 0} \\
\cdashline{1-4}
{\mbf 0} & \mbf A_{32} & \mbf A_{33} & {\mbf 0} \\
\cdashline{1-4}
{\mbf 0} & {\mbf 0} & \mbf A_{43} & \mbf A_{44}
\end{array}
\right).
\end{align}

The blocks match the size of their counterparts in (<ref>). Here, the relatively large matrices $\mbf A_{11}$, $\mbf A_{21}$, $\mbf A_{22}$, $\mbf A_{32}$, and $\mbf A_{33}$ have specific repetitive sparse structures which facilitate memory saving and efficient computation.
Specifically, these contain repeated sub-matrices of dimension $\scard \times \scard$ which are of the following forms:

\begin{align*}
\mbf A_{11}=& \, \left(
\begin{array}{c c c }
\mathbb{D}_{11} &  {\mbf 0} &  {\mbf 0}\\
{\mbf 0} & \ddots &  {\mbf 0}\\
{\mbf 0} & {\mbf 0} & \mathbb{D}_{11}
\end{array}
\right), \quad
\mbf A_{21}= \, \left( \mathbb{D}_{21} \ldots \mathbb{D}_{21} \right), \quad
\mbf A_{22}= \, \mathbb{D}_{22}, \quad
\mbf A_{32}= \, \left(
\mathbb{D}_{32}
\cdots
\mathbb{D}_{32}
\right)^\top, \\
\mbf A_{33} = & \, \left(
\begin{array}{cccccc}
d_{1,1} &  {\mbf 0} &  {\mbf 0} &  {\mbf 0} &  {\mbf 0} & {\mbf 0}\\
g_{2,1}  & d_{2,2} &  {\mbf 0}  & \vdots  & \multirow{2}{*}{$\vdots$}  & \multirow{2}{*}{$\vdots$}\\
f_1 & g_{3,2} & d_{3,3} & {\mbf 0}  &   &  \\
\vdots & f_2 & \ddots & \ddots &  {\mbf 0}  &  \\
\vdots & \vdots &  & \ddots & \ddots  & {\mbf 0} \\
f_1 & f_2 & \cdots & f_{\hcard -1 } & g_{\hcard,\hcard-1} & d_{\hcard,\hcard} \\
\end{array}
\right).
\end{align*}

Using this banded diagonal structure of the LU-decomposition, solving (<ref>) can be done in two steps, (i) solving $\mbf L \mbf y = \widetilde{\mbf b}$, and then (ii) solving $\mbf L^\top \widetilde{\mbf z} = {\mbf y}$. The details are presented below.

* Solve $\mbf L \mbf y = \widetilde{\mbf b}$ with $\mbf y = (\mbf y^{(1)}, \ldots, \mbf y^{(4)})^\top$:

* Solve $\mbf A_{11}{\mbf y}^{(1)} = {\widetilde{\mbf b}}^{(1)}$:

$ {\mbf y}^{(1)}_{i}=  \mathbb{D}_{11}^{-1} {\widetilde{\mbf b}}^{(1)}_i \text{ for } i=1, \ldots, \dcard. $

* Solve $\mbf A_{22}{\mbf y}^{(2)} = {\widetilde{\mbf b}}^{(2)} - \mbf A_{21}{\mbf y}^{(1)}$:

{\mbf y}^{(2)}=  \mathbb{D}_{22}^{-1} \left( {\widetilde{\mbf b}}^{(2)} - \mathbb{D}_{21} \sum_{i=1}^{\dcard} {\mbf y}^{(1)}_{i} \right).$

* Solve $\mbf A_{33} {\mbf y}^{(3)} = {\widetilde{\mbf b}}^{(3)} - \mbf A_{32}{\mbf y}^{(2)}$:

\begin{align*}
{\mbf y}^{(3)}_1 &= d_{1,1}^{-1}\left( {\widetilde{\mbf b}}^{(3)}_1  - \mathbb{D}_{32} {\mbf y}^{(3)}_1 \right), \\
{\mbf y}^{(3)}_2 &= d_{2,2}^{-1} \left(\left( {\widetilde{\mbf b}}^{(3)}_2  - \mathbb{D}_{32} {\mbf y}^{(3)}_2 \right) - g_{2,1}  {\mbf y}^{(3)}_1 \right), \\
{\mbf y}^{(3)}_i &=d_{i,i}^{-1} \left(\left( {\widetilde{\mbf b}}^{(3)}_i  - \mathbb{D}_{32} {\mbf y}^{(3)}_i \right) - g_{i,i-1}  {\mbf y}^{(3)}_{i-1} - \sum_{j=1}^{i-2} f_{j} {\mbf y}^{(3)}_{j} \right)
\text{ for } i =3, \ldots, \hcard.
\end{align*}

* ${\mbf y}^{(4)} = \mbf A_{44}^{-1} \left({\widetilde{\mbf b}}^{(4)} - \mbf A_{43}{\mbf y}^{(3)}\right)$.

* Solve $\mbf L^\top \widetilde{\mbf z} = \mbf y$ with $\wt{\mbf z} = (\wt{\mbf z}^{(1)}, \ldots, \wt{\mbf z}^{(4)})^\top$:

* ${\widetilde{\mbf z}}^{(4)} = \left(\mbf A_{44}^\top\right)^{-1} {\mbf y}^{(4)}$

* Solve $\mbf A_{33}^\top {\widetilde{\mbf z}}^{(3)} = {\mbf y}^{(3)} - \mbf A_{43}^\top {\widetilde{\mbf z}}^{(4)}$:

\begin{align*}
{\widetilde{\mbf z}}^{(3)}_{\hcard} &= d_{ \hcard, \hcard}^{-1} \left(  {\mbf y}^{(3)}_{\hcard} - \left(\mbf A_{43}^\top\right)_{\hcard} {\widetilde{\mbf z}}^{(4)} \right), \\
{\widetilde{\mbf z}}^{(3)}_{\hcard-1} &= d_{ \hcard-1, \hcard-1}^{-1} \left(  {\mbf y}^{(3)}_{\hcard-1}  - \left(\mbf A_{43}^\top\right)_{\hcard-1} {\widetilde{\mbf z}}^{(4)} -g_{\hcard, \hcard-1} {\widetilde{\mbf z}}^{(3)}_{\hcard} \right),
\\
{\widetilde{\mbf z}}^{(3)}_{i} &= d_{i, i}^{-1} \left(  {\mbf y}^{(3)}_{i}  - \left(\mbf A_{43}^\top\right)_{i} {\widetilde{\mbf z}}^{(4)} -g_{i+1, i} {\widetilde{\mbf z}}^{(3)}_{i+1} - (\hcard -i - 1)f_{i} \sum_{j=i+2}^{\hcard} {\widetilde{\mbf z}}^{(4)}_{j}\right)
\end{align*}

for $i = \hcard -2, \ldots, 1$.

* Solve $\mbf A_{22}^\top {\widetilde{\mbf z}}^{(2)} = {\mbf y}^{(2)} - \mbf A_{32}^\top {\widetilde{\mbf z}}^{(3)}$:

${\widetilde{\mbf z}}^{(2)} = \mathbb{D}_{22}^{-1}\left({\mbf y}^{(2)} - \mathbb{D}_{32} \left( \sum_{i=1}^{\hcard} {\widetilde{\mbf z}}^{(3)}_i \right)\right)$

* Solve $\mbf A_{11}^\top {\widetilde{\mbf z}}^{(1)} = {\mbf y}^{(1)} - \mbf A_{21}^\top {\widetilde{\mbf z}}^{(2)}$:

${\widetilde{\mbf z}}^{(1)}_i = \mathbb{D}_{11}^{-1}\left({\mbf y}^{(1)}_i - \mathbb{D}_{21} {\widetilde{\mbf z}}^{(2)}\right)$  for  $i = 1, \ldots, \dcard$.

Here, $(\mbf A_{43})_i$ denotes its $i$th row. The vector arguments are partitioned as

\begin{align*}
{\mbf x}^\top = \l(   (\mbf x^{(1)}_1)^\top , \ldots , (\mbf x^{(1)}_{\dcard})^\top,  \: (\mbf x^{(2)})^\top, \:  (\mbf x^{(3)}_{1})^\top, \ldots,  (\mbf x^{(3)}_{\hcard})^\top, \: (\mbf x^{(4)})^\top \r)^\top,
\end{align*}

where ${\mbf x}^{(j)}_i$ for $j=1, 3$ and ${\mbf x}^{(2)}$ are of length $\scard$ and ${\mbf x}^{(4)}$ is of length $\hcard$.

§.§ Computation of model complexity

The numerator of MC in (<ref>) counts the number of connected components of the given graphs. Typically, the number of connected components of a graph can be found by counting the number of zero-eigenvalue of the graph's Laplacian matrix.

The Laplacian matrix of $\mathfrak{N}_{\text{single}}(r) \cap \wh{\mathfrak{N}}_{D, d}$ has the dimension of $\scard \times \scard$ for each $d$, and it is feasible to compute its eigenvalues and calculate $\mc C (\mathfrak{N}_{\text{single}}(r) \cap \wh{\mathfrak{N}}_{D, d})$.

However, computing $\mc C (\mathfrak{N}_{\text{multi}}(r) \cap \wh{\mathfrak{N}}_{H})$ is not so straightforward since the Laplacian matrix of the graph $\mathfrak{N}_{\text{multi}}(r) \cap \wh{\mathfrak{N}}_H $ has the dimension of $\scard \vert \mc H \vert \times \scard \vert \mc H \vert $, it is not practical to compute its eigendecomposition.

Instead, we propose to obtain $\mc C (\mathfrak{N}_{\text{multi}}(r) \cap \wh{\mathfrak{N}}_{H})$ by first finding the layer-specific connected components, and then coalescing the components of two consecutive layers if they are connected transversely over the layers.

Before describing the proposed method, we introduce some notations relevant to $\mathfrak{N}_{\text{single}}(r) \cap \wh{\mathfrak{N}}_{H, h}$, a single-layer network for each $h \in \mc H$:

\begin{align*}
\wh{\mathfrak{N}}_{H, h} =& \, \l(\mc S, \, \mc E_{H, h}\r) \text{ \ with \ } \mc E_{H, h} = \l\{(s, s'), \, s \ne s': \, \wh\phi^{\hod}_{s, h} = \wh\phi^{\hod}_{s', h} \r\},
\\
\mathfrak{G}_h =& \, \mathfrak{N}_{\text{single}}(r) \cap \wh{\mathfrak{N}}_{H, h},
\\
N_{h} =& \, \mc C (\mathfrak{G}_h),
\\
\mbf{M}^{(0)} =& \, \text{$\vert \mc H \vert \times \scard $ matrix that encodes the layer-specific connected components of }\mathfrak{G}_h,
\\
\mbf L_h =& \, \text{the Laplacian matrix of }\mathfrak{G}_h
\\
\mbf{c}_{h,k}(\mbf M) =& \, \text{the index vector of entries that corresponds to the elements in the $k$th cluster}
\\
&\text{of the $h$th row of a cluster label matrix $\mbf M$} .
\end{align*}

Specifically, for each connected component, the corresponding $(h, s)$ elements of $\mbf M^{(0)}$ take a unique value, and thus the number of unique values in $\mbf M^{(0)}$ agrees with the number of total connected components $\sum_{h=0}^{23} N_h$.

We begin by describing how to construct $\mbf M^{(0)}$. Recall that the number of zero eigenvalues of $\mbf L_h$ corresponds to $N_h$. The eigenvector associated with the zero eigenvalue provides some information of the connected components.
Namely, it is a linear combination of the connected component indicator vectors $\mbf c_{h, k}(\mbf M^{(0)})$.
For an arbitrary index vector $\mbf g = (g_i: \, g_i \in \{1, \ldots, \scard\})$, we denote the indicator vector of $\mbf g$ by $\mbf e_{\mbf g}$:
It has its $g_i$th entry to be one for $i = 1, \ldots, \vert \mbf g \vert$, and all the rest are zeros.
Denote by $\mbf v^{(0)}_{h, i}$ the $i$th eigenvector associated with the zero eigenvalues of $\mbf L_h$ for $i = 1, \ldots, N_h$.
Then, we have $\mbf v^{(0)}_{h,i} = \sum_{k = 1}^{N_h} a^{(h,i)}_{k} \mbf e_{\mbf c_{h, k}(\mbf M^{(0)})}$
for some constants $a^{(h,i)}_k$s. Thus, we utilize the eigenvectors to identify the cluster label of each stations and construct $\mbf M^{(0)}$. Precisely, we find the partition $\mc P$ of the index set $\{1, \ldots, \scard\}$ with the smallest cardinality
such that for each $P$ belonging to the partition,
all the elements of $\mbf v^{(0)}_{h, i}$ located at $P$ take the same value, for all $i = 1, \ldots, N_h$.
This procedure is described in Algorithm <ref>.

Construction of $\mbf M^{(0)}$.

[1]

The Laplacian matrix $\mbf L_h$ for $h \in \mc H$

$\mbf M^{(0)} \leftarrow$ a $\vert {\mc H} \vert \times \scard$-matrix of zeros

$\texttt{cluster\_label} \leftarrow 1$

$h \in \mc H$
Perform eigenvalue decomposition of $\mbf L_h$
$N_h \leftarrow$ the number of zero eigenvalues
$\mbf V \leftarrow$ a $\scard \times N_h$ matrix having the eigenvector $\mbf v^{(0)}_{h,i}$ as its $i$th column
$\texttt{unlabeled} \leftarrow \{1, 2, \ldots, \scard \}$
$s \in \{1, \ldots, \scard\}$
$s \notin \texttt{unlabeled}$
$\texttt{cluster} \leftarrow \{1, \ldots, \scard\}$
$i \in \{1, \ldots, N_h\}$
$ \texttt{value} \leftarrow V_{s,i}$, the $s$th entry of the $i$th eigenvector $\mbf v^{(0)}_{h,i}$
$\texttt{value\_set} \leftarrow \{s': V_{s',i} = \texttt{value}\}$
$\texttt{cluster} \leftarrow \texttt{cluster} \cap \texttt{value\_set}$
$M^{(0)}_{h, s'}\leftarrow \texttt{cluster\_label}$ for $s' \in \texttt{cluster}$
$\texttt{unlabeled} \leftarrow \texttt{unlabeled}\setminus \texttt{cluster}$
$\texttt{cluster\_label} \leftarrow \texttt{cluster\_label}+1$
$\mbf M^{(0)}$

Once the layer-specific cluster label matrix $\mbf M^{(0)}$ is provided, we link the connected components throughout the hourly layers as described in Algorithm <ref>. Within the procedure, the sub-routine given in Algorithm <ref> is utilized, which sequentially links connected components lying in two consecutive layers. In each run, Algorithm <ref> links a component in evaluation to exactly one component in another layer while there can be more than one component that are supposed to be linked. Thus, Algorithm <ref> keeps running the sub-algorithm until there exists no more component left to be linked. One notable feature is that the sub-algorithm is executed twice with different order vector $\mbf{r}$ each time. This is to account for the circular feature of hour-of-a-day. The function $\texttt{is}(\cdot)$ in the sub-algorithm <ref> is defined to return the Boolean of the input statement.

Counting the connected components of $\mathfrak{N}_{\text{multi}}(r) \times \wh{\mathfrak{N}}_H$.

[1]
The layer-specific cluster label matrix $\mbf M^{(0)}$

The estimates $\wh\phi^{\hod}_{s, h}$ for $s \in \mc S$ and $h \in \mc H$
$\mbf M \leftarrow \mbf M^{(0)}$
$\texttt{changed} \leftarrow \texttt{true}$

$\texttt{changed} \leftarrow \texttt{false}$
Make links across the layers in a forward manner by running Algorithm <ref>
with $\mbf r = (0,1,\ldots,23)$
Connect the layers of $h=0$ and $h=23$ in a forward manner by running Algorithm <ref>
with $\mbf r = (23,0,1,\ldots,22)$
the number of unique values in $\mbf M$

Link across ordered layers.

[1]

The estimates $\wh\phi^{\hod}_{s, h}$ for $s \in \mc S$ and $h \in \mc H$

An order vector $\mbf{r}$ of length $\vert \mc H \vert$

An $\vert \mbf r \vert \times \vert \mc S \vert$ cluster label matrix $\mbf{M}$

A Boolean variable

$\mbf M^{\text{old}} \leftarrow \mbf M$
$i \in \{1, \ldots, \vert\mbf{r} \vert\}$
// link the $r_i$th and the $r_{i+1}$th layers:
$\mbf M^{\text{new}} \leftarrow \mbf M^{\text{old}}$
$k \in \{1, \ldots, N_{r_i}\}$
// investigate whether each component is subject to further connection
$\text{ind} \leftarrow \mbf{c}_{r_i,k}(\mbf M^{\text{old}})$
$j \leftarrow 0$
$\texttt{matched} \leftarrow \texttt{false}$
$\texttt{exhausted} \leftarrow \texttt{is}(k \ge \vert \mbf{c}_{r_i,k}(\mbf M^{\text{old}}) \vert)$
$=$  and  $=$
$j \leftarrow j+1$
$c_1 \leftarrow \wh{\phi}^{\hod}_{\text{ind}_j, r_i}$
$c_2 \leftarrow \wh{\phi}^{\hod}_{\text{ind}_{j}, r_{i+1}}$
$\text{ind}^{(j,r_{i+1})} \leftarrow$  the index vector of the cluster to which the $j$th entry of  belongs, in the $r_{i+1}$th row of $\mbf M^{\text{new}}$ i.e. $\mbf c_{r_{i+1}, {k'}}(\mbf M^{\text{ new}})$ for some $k'$ such that $\texttt{ind}_j \in c_{r_{i+1}, {k'}}(\mbf M^{\text{new}})$
$m_1 \leftarrow \min\; {\mbf M^{\text{old}}}_{r_i,\text{ind}}$
$m_2 \leftarrow \min\; {\mbf M}^{\text{new}}_{r_{i+1}, \text{ind}^{(j,r_{i+1})}}$
$\texttt{same\_estimate}\leftarrow \texttt{is}(c_1=c_2)$
$\texttt{same\_label} \leftarrow \texttt{is}(m_1=m_2)$
$\texttt{should\_connect} \leftarrow \texttt{is}(\texttt{same\_estimate} \text{ and not } \texttt{same\_label})$
// connect the to components and update their labels
$m = \min \{m_1, m_2\}$
${\mbf M}^{\text{new}}_{r_i, \text{ind}} \leftarrow m$
${\mbf M}^{\text{new}}_{r_{i+1}, \text{ind}^{(j,r_{i+1})} } \leftarrow m$
$\texttt{changed} \leftarrow \texttt{true}$
$\texttt{matched} \leftarrow \texttt{true}$
$\texttt{exhausted} \leftarrow \texttt{is}(k \ge \vert \mbf{c}_{r_i,k}(\mbf M^{\text{old}}) \vert)$
$\mbf M^{\text{old}} \leftarrow \mbf M^{\text{new}}$
$\mbf M \leftarrow \mbf M^{\text{new}}$
The updated cluster label matrix $\mbf M$

The updated Boolean variable

§ ADDITIONAL DATA ANALYSIS RESULTS

Table <ref> presents the results obtained with $r = 750$,
see Table <ref> in the main text for the results obtained with $r = 1500$.

Estimated coefficients for the covariate effects by the proposed fused Lasso regression method from each fold used in the $7$-fold CV and from the full data when $r = 750$.
For comparison, we also report the estimates obtained with fusion Lasso-only, full-interaction and no-interaction methods.

7cFold     5cFull data

1     2     3     4     5     6     7     Fused     Fusion-only     Lasso-only     Full     No

(lr)1-1 (lr)2-8 (lr)9-9 (lr)10-10 (lr)11-11 (lr)12-12 (lr)13-13
$\alpha$     0.057     0.052     0.051     0.050     0.050     0.054     0.050     0.052     0.052     0.052     0.052     0.052

$\beta^{\text{rain}}$     -2.301     -2.356     -2.380     -2.377     -2.238     -2.182     -2.383     -2.324     -2.325     -2.321     -2.323     -2.310

$\beta^{\text{air}}_1$     0.098     0.141     0.100     0.177     0.124     0.169     0.168     0.140     0.141     0.138     0.138     0.137

$\beta^{\text{air}}_2$     0.116     0.173     0.098     0.177     0.133     0.159     0.187     0.149     0.149     0.145     0.144     0.146

$\beta^{\text{air}}_3$     0.218     0.285     0.250     0.175     0.259     0.321     0.310     0.275     0.277     0.271     0.276     0.243

Figures <ref> and <ref> plot the results obtained with $r = 750$,
see Figures <ref> and <ref> in the main text for the results obtained with $r = 1500$.

Parameter estimates for $\phi^{\hod}_h, \, h \in \mc H$ (left) and $\phi^{\dow}_d, \, d \in \mc D$ (right) from each fold used in the $7$-fold CV and from the full data when $r = 750$.

Estimated station-specific bike demands in log-scale ($\wh{\theta}_s+\wh{\theta}_d + \wh{\theta}_{h} + \wh{\theta}_{s,d}+\wh{\theta}_{s,h} $) from the model fitted with $r = 750$ at 8am on Tuesdays (left), at 8pm on Sundays (middle) and their differences (right).