DANR: Discrepancy-aware Network Regularization

Since existing network-based regularizers rely on known weights on edges, the corresponding solutions get misled when the edge weights are erroneous or unknown. Directly learning the unknown edge weights by using the same regularizer as Eq. (1.1) would lead to an optimization problem:

(2.2)

$\displaystyle\underset{\bm{x},\bm{\omega}}{\operatorname{min~{}}}\quad\sum_{i\in\mathcal{V}}f_{i}(\bm{x}_{i})+\lambda\sum_{(j,k)\in\mathcal{E}}\omega_{jk}\left\|\bm{x}_{j}-\bm{x}_{k}\right\|_{2}$

which yields the trivial all-zero solution of $\bm{\omega}$ and thus becomes unsatisfactory. Instead of imposing more sophisticated model-based or problem-dependent regularization as suggested in [33], we consider an alternative formulation that explicitly accounts for discrepancies between models on adjacent nodes:

(2.3)

$$\underset{\bm{x},\bm{\alpha}}{\operatorname{min~{}}}\quad\sum_{i\in\mathcal{V}}f_{i}(\bm{x}_{i})+\lambda\cdot\mathcal{R}_{S}(\bm{x},\bm{\alpha})\\ \mathcal{R}_{S}(\bm{x},\bm{\alpha})=\mu\sum_{(j,k)\in\mathcal{E}}\omega_{jk}\left\|\bm{x}_{j}+\bm{\alpha}_{jk}-\bm{x}_{k}\right\|_{2}+(1-\mu)\|\bm{\alpha}\|_{1,p}$$

where we denote $\|\bm{\alpha}\|_{1,p}=\sum_{(j,k)\in\mathcal{E}}\|\bm{\alpha}_{jk}\|_{p}$ as the sum of $p$-norms of all $\bm{\alpha}_{jk}$’s. We set a penalty parameter $\lambda\geq 0$ to control the overall strength of network regularization, and a portion parameter $0<\mu<1$ to control the emphasis between the two terms in $\mathcal{R}_{S}(\bm{x},\bm{\alpha})$. We name this formulation in Eq. (2.3) the discrepancy-aware network regularization (DANR).

In the first term of $\mathcal{R}_{S}$, we define a discrepancy-buffering (DB) variable $\bm{\alpha}_{jk}\in\mathbb{R}^{d}$ to denote the preserved discrepancy between local model parameters $\bm{x}_{i}$ and $\bm{x}_{j}$. More specifically, when an edge weight $\omega_{jk}$ is given but potentially imprecise or otherwise corrupted, $\bm{\alpha}_{jk}$ would compensate for abnormally large differences between models $\bm{x}_{j}$ and $\bm{x}_{k}$ to reduce the magnitude of $\omega_{jk}$-weighted edge penalty term in $\mathcal{R}_{S}$. The DB variable provides additional flexibility to its associated models, allowing them to stay adequately close to their own local solutions w.r.t. minimizing local loss functions, and avoid over-penalized consensus. When all $\omega_{jk}$’s are not given, $\bm{\alpha}_{jk}$’s enable solving Eq. 2.3 under anisotropic regularizations even with homogeneous weights.

The second term of $\mathcal{R}_{S}$ is the $\ell_{1,p}$-regularizer (with $1<p<\infty$) [28] that ensures sparsity at the group level. Note that we regard each variable vector $\bm{\alpha}_{jk}$ as a group ($|\mathcal{E}|$ groups in total). The first key intuition here is that the regularization on $\bm{\alpha}$ helps to exclude trivial solutions, that is $\bm{\alpha}_{jk}=\bm{x}_{k}-\bm{x}_{j}$. Second, it allows us to identify a succinct set of non-zero vectors $\bm{\alpha}_{jk}$’s, compensating for possible intrinsic discrepancies between two adjacent nodes.

To sum up, discrepancy-buffering variable $\bm{\alpha}_{jk}$’s are designed to elaborately adjust network regularization strength on all edges, and thus reduce negative effects of the unsquared norm regularizer. We will show in later sections that this modified formulation remains convex and tractable via parallel optimization algorithms on large-scale networks.

In this section, we propose an ADMM-based algorithm for solving the ST-DANR problem in Eq. (3.12) and present the convergence and complexity of the proposed algorithm. The algorithm can be easily adapted to the spatial-only DANR problem in Eq. (2.3) by omitting temporal-related updates.

The ADMM method was originally derived in [8] and has been reformulated in many contexts including optimal control and image processing [23]. The method can be considered as combining augmented Lagrangian methods and the method of multipliers [4]. It aims to solve optimization problems with two-block separable convex objectives in the following form of

(2.4)

$\displaystyle\operatorname{min}~{}f(\bm{x})+g(\bm{z})\qquad~{}s.t.~{}A\bm{x}+B\bm{z}=\bm{c}$

Observe that our proposed objective in Eq. (3.12) has a separable convex objective function: it can be reorganized into two-block separable convex objectives as in Eq. (2.4), for which ADMM methods guarantee convergence to global optimal solutions [7].

More precisely, to fit our problem into the ADMM framework, we first define $\bm{x}$=$\{\bm{x}_{j}\}^{j\in V}$ for model parameter vectors on the given undirected network $\mathcal{G}$, and define consensus variables $\bm{u}$ = $[\{\bm{u}_{jk},\bm{u}_{kj}\}^{(j,k)\in\mathcal{E}}]$ as copies of $\bm{x}$ in the spatial penalty term $\mathcal{R}_{S}(\bm{x},\bm{\alpha})$. Then we rewrite Eq. (3.12) as follows:

(2.5)		$\displaystyle\underset{\bm{x},\bm{u},\bm{\alpha}}{\operatorname{min}}\quad$	$\displaystyle\sum_{j\in\mathcal{V}}f_{j}(\bm{x}_{j})+\lambda_{1}(1-\mu_{1})\|\bm{\alpha}\|_{1,p}$
		$\displaystyle+\lambda_{1}\mu_{1}\sum_{(j,k)\in\mathcal{E}}\omega_{jk}\left\|\bm{u}_{jk}+\bm{\alpha}_{jk}-\bm{u}_{kj}\right\|_{2}$
	s.t.	$\displaystyle\left[\bm{u}_{jk},\bm{u}_{kj}\right]~{}~{}~{}=\left[\bm{x}_{j},\bm{x}_{k}\right],\quad(j,k)\in\mathcal{E}$

in which the first term corresponds to the $f$ block in Eq. (2.4), and the remaining terms correspond to the $g$ block. The equality constraints on Eq. (3.12) are used to force consensus between variables $\bm{x}$ and $\bm{u}$. Next, we derive the augmented Lagrangian of Eq. (E.1):

(2.6)			$\displaystyle\mathcal{L}_{\rho_{1}}(\bm{x},\bm{u},\bm{\alpha},\bm{\delta})=\sum_{j\in V}f_{j}(\bm{x}_{j})+\lambda_{1}(1-\mu_{1})\|\bm{\alpha}\|_{1,p}$
		$\displaystyle\hskip 27.00005pt+\sum_{(j,k)\in E}\Big{(}\lambda_{1}\mu_{1}\omega_{j,k}\|\bm{u}_{jk}+\bm{\alpha}_{jk}-\bm{u}_{kj}\|_{2}$
		$\displaystyle\hskip 54.00009pt+\frac{\rho_{1}}{2}\|\bm{x}_{j}-\bm{u}_{jk}+\bm{\delta}^{u}_{jk}\|_{2}^{2}-\frac{\rho_{1}}{2}\|\bm{\delta}^{u}_{jk}\|_{2}^{2}$
		$\displaystyle\hskip 54.00009pt+\frac{\rho_{1}}{2}\|\bm{x}_{k}-\bm{u}_{kj}+\bm{\delta}^{u}_{kj}\|_{2}^{2}-\frac{\rho_{1}}{2}\|\bm{\delta}^{u}_{kj}\|_{2}^{2}\Big{)}$

where $\bm{\delta}$ = $[\{\bm{\delta}^{u}_{jk},\bm{\delta}^{u}_{kj}\}^{(j,k)\in\mathcal{E}}]$ are scaled dual variables for each equality constraint on the elements of $\bm{u}$. The parameters $\rho_{1}>0$ penalize the violation of equality constraints in the spatial and temporal domain [21] respectively. The iterative scheme of ADMM under the above setting can be written as follows, with $l$ denoting the iteration index:

(2.7)

$\displaystyle\begin{rcases}&\bm{x}^{(l+1)}=\underset{\bm{x}}{\operatorname{argmin~{}}}\mathcal{L}_{\rho_{1}}(\bm{x},(\bm{u},\bm{\alpha},\bm{\delta})^{(l)})\\ &(\bm{u},\bm{\alpha})^{(l+1)}=\underset{\bm{u},\bm{\alpha}}{\operatorname{argmin~{}}}\mathcal{L}_{\rho_{1}}(\bm{x}^{(l+1)},\bm{u},\bm{\alpha},\bm{\delta}^{(l)})\\ &\bm{\delta}^{(l+1)}=\bm{\delta}^{(l)}+(\tilde{\bm{x}}^{(l+1)}-\bm{u}^{(l+1)})\end{rcases}$

where $\tilde{\bm{x}}$=$[\{\bm{x}_{j},\bm{x}_{k}\}^{(j,k)\in\mathcal{E}}]$ is composed of replicated elements in $\bm{x}$, and thus has a one-to-one correspondence with elements in $\bm{u}$.

Because the augmented Lagrangian in Eq. (2.6) has a separable structure as well, we can further split the optimization above over each univariate element in $\bm{x}$ and $\bm{z}$. Next, we provide details for each ADMM update step.

$\bm{x}$-Update. In the first update step of our ADMM updating scheme, we can decompose the problem of minimizing $\bm{x}$ into separately minimizing $\bm{x}_{j}$ for each node $j$ :

$\displaystyle\bm{x}_{j}^{(l+1)}=\underset{\bm{x}_{j}}{\operatorname{argmin~{}}}\Big{(}f_{j}(\bm{x}_{j})+\frac{\rho_{1}}{2}\sum_{k\in N_{j}}\|\bm{x}_{j}-\bm{u}_{jk}^{(l)}+\bm{\delta}_{jk}^{u(l)}\|_{2}^{2}\Big{)}$

($\bm{u},\bm{\alpha}$)-Update. We can further decompose the ($\bm{z},\bm{\alpha}$)-update step in Eq. (E.5) into subproblems on each edge (and its related variables $\bm{u}_{jk},\bm{u}_{kj},\bm{\alpha}_{jk}$) as follows, which can be solved in parallel:

(2.8)

$$(\bm{u}_{jk}^{(l+1)},\bm{u}_{kj}^{(l+1)},\bm{\alpha}_{jk}^{(l+1)})=\\ \operatorname{argmin~{}}\Big{(}\lambda_{1}(1-\mu_{1})\|\bm{\alpha}_{jk}\|_{p}+\lambda_{1}\mu_{1}\omega_{jk}\|\bm{u}_{jk}+\bm{\alpha}_{jk}-\bm{u}_{kj}\|_{2}\\ \hskip 9.00002pt+\frac{\rho_{1}}{2}\|\bm{x}_{j}^{(l+1)}-\bm{u}_{jk}+\bm{\delta}^{u(l)}_{jk}\|_{2}^{2}+\frac{\rho_{1}}{2}\|\bm{x}_{k}^{(l+1)}-\bm{u}_{kj}+\bm{\delta}^{u(l)}_{kj}\|_{2}^{2}\Big{)}$$

Notice that the sum of convex functions which are defined on different sets of variables preserves the convexity. To minimize the objective with $\bm{u}_{jk},\bm{u}_{kj}$ and $\bm{\alpha}_{jk}$ simultaneously, the convexity of Eq. (2.2) motivates us to adopt an alternating descent algorithm, which minimizes each component iteratively with respect to $(\bm{u}_{jk},\bm{u}_{kj})$ and $\bm{\alpha}_{jk}$ while holding the other fixed. In detail, we solve Eq. (2.9) and Eq. (2.10) iteratively until convergence is achieved:

(2.9)			$\displaystyle(\bm{u}_{jk}^{(l^{\prime}+1)},\bm{u}_{kj}^{(l^{\prime}+1)})={\operatorname{argmin}}\Big{(}\lambda_{1}\mu_{1}\omega_{jk}\|\bm{u}_{jk}+\bm{\alpha}_{jk}^{(l^{\prime})}-\bm{u}_{kj}\|_{2}$
		$\displaystyle+\frac{\rho_{1}}{2}\|\bm{x}_{j}^{(l+1)}-\bm{u}_{jk}+\bm{\delta}^{u(l)}_{jk}\|_{2}^{2}+\frac{\rho_{1}}{2}\|\bm{x}_{k}^{(l+1)}-\bm{u}_{kj}+\bm{\delta}^{u(l)}_{kj}\|_{2}^{2}\Big{)}$

(2.10)		$\displaystyle\bm{\alpha}_{jk}^{(l^{\prime}+1)}$	$\displaystyle={\operatorname{{argmin}}}\Big{(}(1-\mu_{1})\|\bm{\alpha}_{jk}\|_{p}$
		$\displaystyle\hskip 18.00003pt+\mu_{1}\omega_{jk}\|\bm{u}_{jk}^{(l^{\prime}+1)}+\bm{\alpha}_{jk}-\bm{u}_{kj}^{(l^{\prime}+1)}\|_{2}\Big{)}$

$\bm{\delta}$-Update. In the last step of our ADMM updating scheme, we have fully independent update rules for each scaled dual variable $\bm{\delta}_{jk}^{u}$ as follows:

(2.11)

$\displaystyle{\bm{\delta}_{jk}^{u}}^{(l+1)}$

$\displaystyle={\bm{\delta}_{jk}^{u}}^{(l)}+(\bm{x}_{j}^{(l+1)}-\bm{u}_{jk}^{(l+1)})$

Finally, we present the pseudo-code of our DANR approach in Appendix D.

Stopping Criterion and Global Convergence. For our ADMM iterative scheme in Eq. (E.5), we use the norm of primal residual $\bm{r}^{(l)}=\tilde{\bm{x}}^{(l)}-\bm{u}^{(l)}$ and dual residual $\bm{s}^{(l)}=\rho_{1}(\bm{u}^{(l)}-\bm{u}^{(l+1)})$ as the termination measure. The optimality condition [4] of ADMM shows that if both residuals are small then the objective suboptimality must be small, and thus suggests $\|\bm{r}^{(l)}\|_{2}\leq\epsilon_{1}\wedge\|\bm{s}^{(l)}\|_{2}\leq\epsilon_{2}$ as a reasonable stopping criterion. Convex subproblems in $\bm{x}$-update and $(\bm{u},\bm{\alpha})$-update need iterative methods to solve as well. The stopping criterion for these subproblems is naturally to keep iteration differences below thresholds, i.e. in $(\bm{u},\bm{\alpha})$-update, we require ${\small\|\Delta\bm{u}^{(l^{\prime})}\|_{2}\leq\epsilon^{\prime}}$ and ${\small\|\Delta\bm{\alpha}^{(l^{\prime})}\|_{2}\leq\epsilon^{\prime}}$.

Computational Complexity. Let $N_{c}$ denote the number of iterations that ADMM takes to achieve an approximate solution $\hat{\bm{x}}$ with an accuracy of $\epsilon_{x}$. Based on the convergence analysis in [14], the time complexity scales as $O(1/\epsilon_{x})$, which in our problem mostly depends on the properties of cost functions $f_{jt}$’s. Assume that all convex subproblems in $x$-update and $(u,\alpha)$-update are solved by general first-order gradient descent methods that take $N_{x}$ and $N_{\alpha}$ iterations to converge respectively, the overall complexity of the algorithm is therefore $O\big{(}N_{c}\big{(}N_{x}|\mathcal{V}|+N_{\alpha}|\mathcal{E}|+|\mathcal{E}|\big{)}\big{)}$.

Following previous work [11], we first experiment with a synthetic network in which each node attempts to solve its own support vector machine (SVM) classifier, but only has a limited amount of data to learn from. We further split the nodes into clusters where each cluster has an underlying model, i.e., nodes in the same cluster share the same underlying model. Our aim is to learn accurate models for each node by leveraging their network connections and limited observations. Note that the neighbors with different underlying models provide misleading information to each other, which potentially harms the overall performance. Therefore, we later investigate the robustness of the DANR with a varying ratio of such malicious edges.

Synthetic Network Generation. We create a network of 100 nodes, split into 5 ground-truth communities $C^{1}$, $C^{2}$, $\cdots$, $C^{5}$. Each of these communities consists of 20 nodes. We denote the community of a node $i$ with $C_{i}$. Next, we form the network connections by adding edges between nodes based on the following criterion: The probability of adding an edge between any node pair $(i,j)$ is 0.5 if $C_{i}=C_{j}$, i.e., the edge connects nodes from the same community (intra-community edge). Otherwise, we set the probability as 0.02 if $C_{i}\neq C_{j}$, i.e., the edge connects nodes from different communities (inter-community edge). The resulting synthetic network $G_{0}(\mathcal{V},\mathcal{E})$ has 100 nodes and 574 edges, with 82% of them being intra-community edges. Each ground-truth community $C^{k}$ has an underlying classifier model $\bm{X}^{k}\in R^{10}$, drawn independently from a normal distribution of zero mean and unit variance, where $k\in[1,\cdots,5]$. Next, we generate 5 random training example pairs $(\bm{w},y)$ per node, where $\bm{w}\in R^{10}$ denotes an observation and $y\in\{-1,1\}$ denotes the ground-truth value to estimate. The ground-truth value $y$ for each observation $\bm{w}$ is then computed using the underlying classifier of the community that a node belongs to:

(4.15)

$\displaystyle y=\text{sign}(\bm{w}^{T}\bm{X}^{k}+\eta)$

Optimization Problem. The objective function $f_{i}$ on each node is formulated by the soft margin SVM objective, where node $i$ estimates its model (i.e. separating hyperplane) $\bm{x}_{i}\in R^{10}$ using its five training examples $(\bm{w},y)$ as follows:

$\displaystyle\min\Big{(}\frac{1}{2}\bm{x}_{i}^{T}\bm{x}_{i}+C\sum_{l=1}^{5}\xi_{l}\Big{)}\quad\text{s.t. }y^{(l)}(\bm{x}_{i}^{T}\bm{w}^{(l)})\geqslant(1-\xi_{l}),~{}\xi_{l}\geqslant 0$

Test Results. In following subsections, we first thoroughly inspect DANR and its prototype method NL on synthetic network $G_{0}$, in order to show the benefits of introducing discrepancy-buffering variable $\bm{\alpha}$. Later in §4.1.2, we evaluate the performance and robustness of DANR and all three baseline methods under more noisy scenarios. To assess the performance, prediction accuracies are computed over a separate set of 1000 test pairs (10 per node). The test pairs are again randomly generated using Eq. (4.15).

In this section, our goal is to jointly (i) estimate the rental prices of houses by leveraging their geological location and the set of features; (ii) discover boundaries between neighborhoods. The intuition is that the houses in the same neighborhood often tend to have similar pricing models. However, such neighborhoods can have complex underlying structures (as described later in this section), which imposes additional challenges in learning accurate models.

Dataset. We experiment with two largely populated areas in Northern California: the Greater Sacramento Area (SAC) and the Bay Area (BAY). The anonymized data is provided by the property management software company Appfolio Inc, from which we further sample houses with at least one signed rental agreement during the year 2017. The resulting dataset covers 3849 houses in SAC and 1498 houses in BAY. Concerning the houses that have more than one rental agreement signed during 2017, we average the rental prices listed in all the agreements. Each house holds the information about its location (latitude/longitude), number of bedrooms, number of bathrooms, square footage, and the rental price. We regard these areas (SAC and BAY) as two separate datasets for the remainder of this section. We also randomly split 20 % of the houses in each dataset for testing.

Network Construction. After excluding test houses from both datasets, we construct two networks (one for each dataset) based on the houses’ locations. An undirected edge exists between node $i$ and $j$, if at least one of them is in the set of 10 nearest neighbors of the other. (Note that node $j$ being one of the 10 nearest neighbors of $i$ doesn’t necessarily imply that node $i$ is also in the set of nearest neighbors of $j$.) In the resulting networks, the average degree of a node is 12.16 for SAC and 12.04 for BAY. Additionally, we construct two versions of these networks, weighted and unweighted. While the weighted network has edge weights inversely proportional to the distances between houses, the unweighted network ignores the proximity between houses, and thus weights on all the edges are 1.

Optimization Problem. The model at each house estimates its rental price by solving a linear regression problem. More specifically, at each node $i$ we learn a 4-dimensional model $\bm{x}_{i}=[a_{i},b_{i},c_{i},d_{i}]$. $\bm{x}$ simply represents the coefficients of each feature, which later is used to estimate the rental price $p_{i}$ as follows:

	$\displaystyle\overline{p_{i}}=a_{i}\cdot(\#bedrooms)+b_{i}\cdot(\#bathrooms)$
	$\displaystyle+c_{i}\cdot(squarefootage)+d_{i},$

where $d_{i}$ is the bias term. The training objective is

$\displaystyle\min_{\bm{x,\alpha}}~{}\sum_{i\in\mathcal{V}}\|\overline{p_{i}}-p_{i}\|_{2}^{2}+c\cdot\|\bm{x}_{i}\|_{2}^{2}+\lambda\cdot\mathcal{R}_{S}(\bm{x},\bm{\alpha})$

where $\mathcal{R}_{S}$ represents the proposed network regularization term (see Eq. 2.3), while $c$ is the regularization weight to prevent over-fitting.

Test Results. Once training converges, we predict the rental prices on the test set. To do that, we connect each house in the test set to its 10 nearest neighbors in the training set. We then infer the new model $\bm{x}_{j}$ by taking the average of the models on its neighbors: $\bm{x}_{j}=(1/|N(j)|)\sum_{k\in N(j)}\bm{x}_{k}$. The model $\bm{x}_{j}$ is further used to estimate the rental price of the corresponding house. Alternatively, one can also infer $\bm{x}_{j}$ by solving $\min_{\bm{x}_{j}}\sum_{k\in N(j)}\|\bm{x}_{j}-\bm{x}_{k}^{*}\|_{2}$, while keeping the models on neighbors fixed. However, we empirically find that simply averaging the neighbors’ models performs better for both methods in this particular setting. We compute Mean Squared Error (MSE) over test nodes to evaluate the performance.

Table 3 displays the test results of eight different settings on both datasets. As shown, local & global estimations and rMTFL method produce high errors for both datasets. We further apply FORMULA and Network Lasso (NL) methods to both weighted and unweighted versions of the networks, while the proposed method is only applied to the unweighted version. Notice that the network weights are irrelevant for the local & global estimation settings and the rMTFL method. Intuitively, both FORMULA and NL performs better on the weighted setting compared to the unweighted setting for both datasets. As shown, the rental price estimation errors achieved by the NL are 0.4173 (weighted) vs. 0.4392 (unweighted) for BAY and 0.3022 (weighted) vs. 0.3273 (unweighted) for SAC. These results suggest that the pre-defined weights on these networks help to learn more accurate models on houses.

However, we further argue that although such pre-defined weights improve the overall clustering performance, they don’t account for more complex clustering scenarios, e.g., two nearby houses falling into different school districts or some houses having higher values compared to their neighbors due to geography, e.g., having a view of the city. That being said, DANR outperforms all the other baselines and achieves the smallest errors for both datasets; 0.4106 for BAY and 0.2978 for SAC. Especially, the DANR (unweighted) even outperforms the weighted version of the baseline approaches by a notable margin. This reveals that the DANR indeed accounts for such heterogeneities in data and provides enhanced clustering of the network.

Figure 5 shows examples of two complex scenarios that are captured by the DANR from the SAC network. We use a color code to represent the clusters in the network, where the same colored houses ($i$, $j$) are in consensus on their models ($\bm{x}_{i}=\bm{x}_{j}$). In the left subfigure, the house shown in yellow uses a different model than all of its neighbors. Interestingly, it resides at the border of three different area codes. The area code for this house is 95864, while the area code on its west is 95825 and 95826 on its south-east. Additionally, we observe similar heterogeneous behaviors in some houses that are near creeks, lakes, and rivers. As an example, Figure 5 (right) displays a creek side house (colored in blue), which again uses a different underlying model from its neighbors.

We now evaluate our method in spatio-temporal setting, where the aim is to dynamically estimate the water quality of Northeastern US lakes over years. We follow the same procedure in §4.2, but have an additional temporal regularization term in our objective that allows nodes to obtain signals from past snapshots of the network, along with their neighbors in the current snapshot.

Dataset and Network. We experiment with the geospatial and temporal dataset of lakes in 17 states in the Northeastern United States, called LAGOS [26]. The dataset holds extensive information about the physical characteristics, various nutrient concentration measurements (water quality), ecological context (surrounding land use/cover, climate etc.), and location of lakes; from which we select the following available set of features: {max depth, surface area, elevation, annual mean temperature, % of agricultural land, % of urban land, % of forest land, % of wetland}. The water quality metric to estimate is the summerly mean chlorophyll concentration of the lakes. We represent the feature vector with $\bm{w}\in\mathbb{R}^{8}$ and the water quality score with $y\in\mathbb{R}$.

During our experiments, we consider a 10-year period between 2000 and 2010. Due to sparsity in the data, we allow 2 years range between the two consecutive snapshots of the network. This results in total of 1039 lakes with all the aforementioned features and the water quality measurements available for each of the selected years. After randomly splitting 20% of the lakes for testing, we build our network by using the latitude/longitude information of lakes, where each lake (node) is connected to 10 nearest lakes.

Optimization Problem. After constructing the network, we now aim to dynamically estimate the water quality ($y_{i,t}$) of lakes by using their feature vectors ($\bm{w}_{i,t}$). More formally, for each year $t\in[2000,2002,\cdots,2010]$, we learn a model $\bm{x}_{i,t}\in\mathbb{R}^{8}$ per node by solving the following spatio-temporal problem in a streaming fashion: (4.16)

$\displaystyle\min_{\bm{x,\alpha,\beta}}~{}\sum_{i\in\mathcal{V}}f(\bm{x}_{i,t},\bm{w}_{i,t},y_{i,t})+\lambda_{1}\cdot\mathcal{R}_{S}^{t}(\bm{x},\bm{\alpha})+\lambda_{2}\cdot\mathcal{R}_{T}^{t}(\bm{x},\bm{\beta})$

Test results. For each year, we first solve the spatial problem (Eq. (4.16) without the temporal term $\mathcal{R}_{T}$) and report the Mean Squared Errors in Table 2. Note that the test nodes are again inferred by averaging the neighbors’ models in the current snapshot. Later, in order to successfully assess the improvements gained by the temporal discrepancy-buffering variables ($\bm{\beta}$), we solve the above spatio-temporal problem with three different temporal regularizers:

DANR+T-SON and DANR+T-SOS approaches simply apply two widely adopted temporal regularizers on the temporal edges (the sum-of-norms regularizer [19] and sum-of-squares reqularizer [6, 29] respectively), while ST-DANR includes the discrepancy-buffering variables on both spatial and temporal edges. Recall that we solve the Eq. (4.16) in a streaming fashion for simplicity, i.e., each snapshot of the network solves the problem while holding the models learned on the previous snapshot (if available) fixed. Potentially, the performance can further be improved by allowing updates on the previous snapshots.

As we can see from the Table 2, leveraging the temporal connections between the two consecutive snapshots significantly improves the performance. The DANR+T-SON outperforms the DANR by 8%-14%, while DANR+TSOS outperforms the DANR by 10%-15%. Moreover, the ST-DANR consistently outperforms both baselines for all the years. This confirms that the proposed formulation accounts for the heterogeneous nature of the temporal networks where some group of nodes exhibits different evolution patterns than the others.

Figure 6 further displays the models (color-coded) learned on three consecutive snapshots, corresponding to years 2000, 2002 and 2004. In 2000, the formed clusters don’t go much beyond the boundaries of the states, resulting in coarse clustering of the network. Yet some states such as Missouri and Iowa are grouped into one cluster (shown in green). This suggests that accurately estimating the chlorophyll concentration of lakes based on the selected set of features is indeed challenging. Specifically, the estimation problem depends on several other external factors that potentially affect the volume of plants and algae in lakes; cultural eutrophication, nutrient inputs from human activities to name a few [26]. However, as it can be seen from the models learned in 2002 and 2004, the clustering performance improves once we allow temporal regularization to synchronize models between two consecutive snapshots, which is analogous with the reduction in errors over years as reported in Table 2. Overall, we observe a split of clusters over time, e.g., in Wisconsin, Minnesota, Iowa and New England. Additionally, a dark brown cluster in south Missouri begins to appear in 2002 and further spreads north in 2004. This indicates that by leveraging the temporal connections, the ST-DANR learns improved models and clustering while allowing for heterogeneity in group level evolutions of nodes.