# Graph-based Virtual Sensing from Sparse and Partial Multivariate
Observations

Giovanni De Felice1 , Andrea Cini2, Daniele Zambon2, Vladimir V. Gusev1,
Cesare Alippi2,3
1University of Liverpool
2The Swiss AI Lab IDSIA & Università della Svizzera italiana
3Politecnico di Milano
Correspondence to<EMAIL_ADDRESS>
###### Abstract

Virtual sensing techniques allow for inferring signals at new unmonitored
locations by exploiting spatio-temporal measurements coming from physical
sensors at different locations. However, as the sensor coverage becomes sparse
due to costs or other constraints, physical proximity cannot be used to
support interpolation. In this paper, we overcome this challenge by leveraging
dependencies between the target variable and a set of correlated variables
(covariates) that can frequently be associated with each location of interest.
From this viewpoint, covariates provide partial observability, and the problem
consists of inferring values for unobserved channels by exploiting
observations at other locations to learn how such variables can correlate. We
introduce a novel graph-based methodology to exploit such relationships and
design a graph deep learning architecture, named GgNet, implementing the
framework. The proposed approach relies on propagating information over a
nested graph structure that is used to learn dependencies between variables as
well as locations. GgNet is extensively evaluated under different virtual
sensing scenarios, demonstrating higher reconstruction accuracy compared to
the state-of-the-art.

##  1 Introduction

Spatio-temporal data analysis plays a significant role in applied and
fundamental domains where systems evolve over both space and time. These
include, for example, environmental science, urban planning, and epidemiology
(Cressie & Wikle, 2015; Wang et al., 2020). In practice, the acquisition of
spatio-temporal data is inevitably affected by all of the typical issues
arising in real-world scenarios. Sensor and communication failures, for
instance, can result in partial or even complete data loss at certain
locations (Little & Rubin, 2019). Furthermore, the deployment of physical
sensors can incur high costs that often limit the number of monitored
locations. Nonetheless, inferring the values of observables at new target
locations holds a significant value for analysis and, eventually, decision-
making (Zhang et al., 2022; Hu et al., 2023). To address these challenges,
virtual sensing (Liu et al., 2009) (also known as spatio-temporal kriging)
(Stein, 1999) has emerged as a viable solution. Virtual sensing consists of
reconstructing complete signals at unobserved locations by utilizing data
gathered from monitored locations during the same time frame (Paepae et al.,
2021; Brunello et al., 2021). As an example, virtual sensing techniques can be
used to estimate solar irradiance at a particular location by exploiting
neighboring sensors (Jayawardene & Venayagamoorthy, 2016), allowing, for
instance, to plan for the installation of a new PV power plant. Most works
rely on spatial correlations and on interpolating observations at neighboring
locations (Appleby et al., 2020; Wu et al., 2021a; b; Zheng et al., 2023).
However, as the sensor coverage becomes sparser due to costs and/or other
(e.g., ecological) concerns, the number of neighbors from which to extract
useful information becomes limited or even null. This makes the spatial
position less informative and limits the effectiveness of existing approaches.
One such prominent example is material weathering testing (De Felice et al.,
2022), where such sparsity, together with the total lack of any historical
data regarding the variable to infer at the target location, makes virtual
sensing bound to fail. Despite being quite common in the real world, this
scenario of sparse data collection is still mostly overlooked in the
scientific literature.

In this work, we focus on the multivariate setting and show how leveraging
correlations between variables can mitigate the challenges posed by sparse
data scenarios. We consider sparsely distributed sensors, each of which is
associated with a multivariate time series (MTS). MTS can be partly
unobserved, i.e., completely lack some of the channels. In such settings, the
problem becomes that of reconstructing missing channels (target variables)
from the available ones (covariates) at the same location and from
observations at related sensors (Fig. 1).

Figure 1: Sparse multivariate virtual sensing; available data (yellow),
missing data (white), and predictions (red).

This scenario is typical of many real-world applications; in fact, variables
correlated to the reconstruction target are often available at each location,
e.g., satellite weather data can be used as covariates to infer photovoltaic
energy production or material degradation. The approach developed in this
paper follows from two observations. First, the mutual dependencies between
the target variables and the covariates can be learned from the dynamics
observed at locations where the target variable is being monitored. Such
dependencies can then be leveraged at the target location, where covariates
provide partial observability, to infer the missing variable. Second, the
target variables can be reconstructed by integrating information from sensors
that, although physically distant, can be considered close (similar) in a
certain functional latent space. Such similarities can be learned from the
data. To solve the task while taking advantage of the above considerations, we
design a novel graph-based framework allowing for modeling dependencies
between variables and propagating information across locations through graph
convolutional layers (Bronstein et al., 2021). Graph neural networks (GNNs)
(Scarselli et al., 2008; Bacciu et al., 2020) indeed allow for exploiting
dependencies as architectural biases in deep learning architectures and have
recently found widespread success in time series processing (Li et al., 2018;
Cini et al., 2023a; Jin et al., 2023). In particular, we propose a novel graph
deep learning method, named Graph-graph Network (GgNet), designed to propagate
information through a nested graph structure encompassing both sensors and
channels. Relying on such a structure, the architecture of the resulting model
is biased toward propagating representations by exploiting learned
dependencies among both variables and locations. Summarizing, our novel
contributions are as follows:

  * •

we introduce a general methodology for performing multivariate virtual
sensing, specifically designed to handle scenarios where available sensors do
not provide adequate spatial coverage (Sec. 4.1);

  * •

we propose a nested graph representation allowing for effectively modeling
dependencies within multivariate spatio-temporal data (Sec. 4.2);

  * •

we design GgNet, the first spatio-temporal graph neural network explicitly
tailored for multivariate virtual sensing from sparse observations (Sec. 4.3);

  * •

we carry out an extensive empirical evaluation by exploring different use
cases and assessing the performance of the proposed method against the state-
of-the-art (Sec. 5).

The proposed framework constitutes a first attempt at tackling a challenging
setting for deep virtual sensing methods and arguably constitutes a key
contribution to the methodological advancement of the field, as well as a
powerful tool in practical applications.

##  2 Problem formulation

This section introduces the notation and formalizes the multivariate virtual
sensing problem, with a focus on settings where sensors provide a sparse
coverage of the area of interest.

##### Notation

We indicate scalar variables and indices as lowercase $k$, constants as
uppercase $K$, vectors as bold lowercase x, matrices as bold uppercase X, and
higher order tensors and sets with calligraphic $\mathcal{X}$. We use the
notation $\vec{\textbf{x}}$ and $\vec{\textbf{X}}$ to indicate univariate and
multivariate time series, respectively. Association with the n-th location is
indicated with the notation $x[n]$. Measurements at time step $t$ are
indicated with subscript $x_{t}$; channel (variable) $d$ of a multivariate
observation is indicated with superscript $x^{d}$.

##### Multivariate spatio-temporal data

Spatio-temporal data refers to a collection of temporal observations coming
from $N$ distinct spatial locations with coordinates
$\\{\textbf{q}[n]\\}^{N}_{n=1}$ in a generic domain $\Omega$, e.g., the one
induced by the geographic placement of sensors. For each spatial location $n$,
consider a multivariate time series (MTS)
$\vec{\textbf{X}}[n]\in\mathbb{R}^{T[n]\times D}$, where $T[n]$ and $D$ are,
respectively, the number of time steps at location $n$, and the number of
channels, i.e., observable variables. If observations are synchronous and an
equal number of time steps are recorded across locations, the set of
measurements can be represented as a 3-D tensor
$\mathcal{X}\in\mathbb{R}^{N\times T\times D}$. Otherwise, the tensor
representation can usually be obtained by padding or interpolating the missing
observations (Lepot et al., 2017).

##### Univariate virtual sensing

Consider a set of $N$ sensors positioned at distinct locations in a spatial
domain $\Omega$ where a subset of $\bar{N}<N$ univariate time series
$\\{\vec{\textbf{x}}[n]\in\mathbb{R}^{T}\\}_{n=1}^{\bar{N}}$ are observed,
while the remaining are missing. In line with previous works (Cao et al.,
2018; Cini et al., 2022), we model data availability with a binary mask
$m[n]\in\\{0,1\\}$, which indicates if the $n$-th location is observed
($m[n]=1$) or missing ($m[n]=0$). The task of _univariate virtual sensing_
consists of estimating
$p(\vec{\textbf{x}}[n]\,|\,\textbf{q}[n],\,\mathcal{X})$, where time series
$\vec{\textbf{x}}[n]$ is unavailable (i.e., $m[n]=0$), and the position
$\textbf{q}[n]$ of the associated sensor is given together with the
observation set
$\mathcal{X}=\\{(\vec{\textbf{x}}[n],\,\textbf{q}[n])\,|\,m[n]=1\\}_{n=1}^{N}$.
Note that spatial coordinates are assumed to be available everywhere.

##### Multivariate virtual sensing

The problem can easily be extended to the multivariate case. Here, consider a
set of multivariate sensors $\\{\vec{\textbf{X}}[n]\in\mathbb{R}^{T\times
D}\\}_{n=1}^{N}$ where individual channels are missing at some locations. We
can denote channel availability by extending the previous binary mask to
$m^{d}[n]\in\\{0,1\\}$, which indicates if the $d$-th channel is available at
the $n$-th location. The observation set now consists of all the observed
channels at all locations
$\mathcal{X}=\\{(\vec{\textbf{x}}^{\,d}[n],\,\textbf{q}[n])\,|\,m^{d}[n]=1\\}_{n=1}^{N}$.
The task of _multivariate virtual sensing_ (MVS) then consists of modeling:

$p(\vec{\textbf{x}}^{\,d}[n]\,|\,\textbf{q}[n],\,\mathcal{X})$ (1)

for every pair ($n$, $d$) such that $m^{d}[n]=0$. Note that virtual sensing
is, a priori, a much more challenging task than the ordinary missing data
imputation problem, due to the absence of any historical data to infer both
the scale and the dynamics of the missing signal. When addressing this
problem, the sensors’ spatial proximity plays a pivotal role. If the sensors
densely cover space, it is possible, in both the univariate and multivariate
settings, to leverage the spatial proximity to available observations.
However, the problem is far more challenging in sparse settings: if the
problem is univariate, no observations are available at the target or any
neighboring location. Conversely, the multivariate settings might allow for
dealing with the virtual sensing problem by exploiting dependencies on
variables other than the target one.

##  3 Related methods

Many approaches to kriging and virtual sensing have been proposed in
geospatial analysis and machine learning literature (Cressie & Wikle, 2015).
Notably, among existing methods, some rely on Gaussian processes (Luttinen &
Ilin, 2012) and tensor decomposition (Bahadori et al., 2014). More related to
our approach, Cini et al. (2022) and Marisca et al. (2022) introduced
imputation methods based on GNNs and showed that such methods can perform
virtual sensing as well. Wu et al. (2021a) and Zheng et al. (2023) directly
tackle the virtual sensing problem exploiting inductive GNNs. Indeed, GNNs are
becoming popular in time series imputation and reconstruction (Ye et al.,
2021; Kuppannagari et al., 2021; Chen et al., 2022; Jin et al., 2023).
However, all these approaches focus on univariate settings. As a consequence,
they can only rely on the spatial proximity of the available sensors to
perform the reconstruction and/or do not target sparser settings. Besides
virtual sensing, deep learning methods have been widely applied in the context
of MTS imputation by relying on a diverse range of model architectures, e.g.,
autoregressive RNNs (Cao et al., 2018), attention-based models (Shukla &
Marlin, 2021; Du et al., 2023), generative adversarial networks (Yoon et al.,
2018a; Luo et al., 2019) and diffusion models (Tashiro et al., 2021). Although
such approaches can in principle be used to perform virtual sensing in the
multivariate case, they cannot take spatial correlations directly into
account. To the best of our knowledge, the proposed GgNet is the first graph-
based neural network architecture designed to tackle MVS in sparse scenarios
by exploiting learned relationships across both observable variables and
locations. We refer to Appendix A for a more thorough discussion of the
relevant literature for spatio-temporal kriging and MTS imputation.

##  4 Methodology

Here, we introduce a general framework for solving the MVS problems formalized
in Sec. 2, regardless of the spatial coverage offered by the available
sensors. We begin by presenting a conceptualization of the inference problem
as dealt with in our framework. Then, we introduce the details of the proposed
novel architecture, aligned with such conceptualization.

###  4.1 Conceptualization

We identify two fundamental types of dependencies characterizing the MVS task,
i.e., the relationships between channels at the same location (_intra_
location) and among observations at different locations (_inter_ location). We
model such dependencies by considering location-channel $(n,d)$ pairs and
defining the sets

$\displaystyle\mathcal{Y}^{d}[n]$
$\displaystyle=\left\\{\vec{\textbf{x}}^{\,\delta}[\nu]\ \big{|}\
\delta=d\,\land\,\nu=n\,\land\,m^{\delta}[\nu]=0\right\\}$ (target set)
$\displaystyle\mathcal{T}^{d}[\setminus n]$
$\displaystyle=\left\\{\vec{\textbf{x}}^{\,\delta}[\nu]\ \big{|}\
\delta=d\,\land\,\nu\neq n\,\land\,m^{\delta}[\nu]=1\right\\}$ (observed
target set) $\displaystyle\mathcal{C}^{\setminus d}[n]$
$\displaystyle=\left\\{\vec{\textbf{x}}^{\,\delta}[\nu]\ \big{|}\ \delta\neq
d\,\land\,\nu=n\,\land\,m^{\delta}[\nu]=1\right\\}$ (intra-location covariate
set) $\displaystyle\mathcal{C}^{\setminus d}[\setminus n]$
$\displaystyle=\left\\{\vec{\textbf{x}}^{\,\delta}[\nu]\ \big{|}\ \delta\neq
d\,\land\,\nu\neq n\,\land\,m^{\delta}[\nu]=1\right\\}$ (inter-location
covariate set)

where $\setminus n$ indicates all indices except $n$. Sets
$\mathcal{Y}^{d}[n],\mathcal{T}^{d}[\setminus n],\mathcal{C}^{\setminus d}[n]$
and $\mathcal{C}^{\setminus d}[\setminus n]$ encompass, respectively, the
target variable at the target location, the target variables at all other
locations except $n$, the observed variables at the target location, and the
observed variables at all other locations except $n$. By considering the above
sets, we distinguish between the types of dependencies that allow for
reconstructing a missing channel $\mathbf{x}^{d}[n]$ in the target set.

T$\rightarrow$Y:

Relations between observed targets at other locations
($\mathcal{T}^{d}[\setminus n]$) and the target $\mathbf{x}^{d}[n]$.

C$\rightarrow$Y:

Relations between observed covariates at the target location
($\mathcal{C}^{\setminus d}[n]$) and the target $\mathbf{x}^{d}[n]$.

The subtask of learning from T$\rightarrow$Y dependencies is similar to
learning how to exploit observations at related locations typical of
univariate virtual sensing. To weigh the importance of different observations
in the observed target set w.r.t. the target, a notion of similarity between
locations is required. Kriging methods address it by relying on spatial
proximity (dense scenario), e.g., as derived from the sensor positions
$\\{\mathbf{q}[n]\\}$, to weigh the reciprocal importance of observations at
different locations. In contrast, with the sparse settings in mind, we model
T$\rightarrow$Y dependencies in a data-driven fashion, learning a similarity
score for each pair of locations. In particular, as detailed in Sec. 4.2, we
associate a representation (node embedding) (Cini et al., 2023a) to each
location and learn a score function taking as input pairs of such
representations. At the same time, these are used as additional learned
covariates to tailor C$\to$Y in modeling location-specific dynamics. Such a
framework allows for the extension of virtual sensing to generic MTS datasets,
as it does not rely on pre-defined spatial relationships such as physical
proximity. In the following, we design a model that directly accounts for both
types of dependencies with architectural inductive biases that align the
processing with the modeling of intra and inter-location dependencies.

###  4.2 Nested graph structure

Figure 2: Dataset as a collection of multivariate time series (left) and its
nested graph representation (right). Each location is represented as a node in
the inter-location graph (G), capturing generic relations between locations.
Within each node in G, a smaller intra-location graph (g) models dependencies
between the channels.

To take advantage of both T$\to$Y and C$\to$Y relations, we learn a nested
graph structure composed of an inter-location graph and an intra-location
graph (Fig. 2). The inter-location graph connects different spatial locations,
while the intra-location one explicits dependencies among channels.

Inter-location graph $G$:

We associate each sensor location to a node of a graph $G$, which we call
_inter-location_ graph, and model the association between different sensors as
edges of $G$, thus accounting for the T$\to$Y modeling.

As physical proximity cannot be exploited in sparse settings, we learn the
graph topology ($N\times N$ adjacency matrix $\mathbf{A}_{G}$) from the data
to account for dependencies among sensors that can be far apart in space. In
particular, following previous works on local representations (Cini et al.,
2023a), we associate each $n$-th node in $G$ with a learnable (static)
embedding $\mathbf{e}_{G}[n]\in\mathbb{R}^{h_{e,G}}$. These are learned as
parameters, jointly with the network weights, and are used to estimate
relationships relevant to the downstream MVS task. The weighted adjacency
matrix is then obtained from the similarities between embeddings. In
particular, we model $\textbf{A}_{G}$ as

$\textbf{A}_{G}=Softmax(\,\textnormal{MLP}_{1}(\textbf{E}_{G})\,\cdot\,\textnormal{MLP}_{2}(\textbf{E}_{G})^{T})$
(2)

where the matrix $\textbf{E}_{G}\in\mathbb{R}^{N\times h_{e,G}}$ contains the
learnable node embeddings for all nodes in $G$. Convolutions on such a graph
directly account for T$\rightarrow$Y dependencies by propagating information
across locations. Similarly, we learn an adjacency matrix for intra-location
dependencies.

Intra-location graph $g$:

We identify each of the $D$ channels, at every $n$-th sensor, with a node of a
second graph, which we name intra-location graph $g$, accounting for the
relationships between channels (C$\to$Y).

Similarly to inter-location case, we associate each $d$-th node in $g$ with a
learnable embedding $\mathbf{e}_{g}[d]$. We learn adjacency matrix
$\textbf{A}_{g}\in\mathbb{R}^{D\times D}$ of $g$ as a function of a $D\times
D$ matrix $\Phi$ of free parameters, treated as edge scores. As the physical
quantities considered at each location are the same, we assume $g$ to be the
same in all locations.

Note that, in this paper, for simplicity, we do not consider graphs $G$ and
$g$ that change over time and do not constrain $G$ and $g$ to be discrete
binary objects. That being said, more advanced graph learning methods (e.g.,
Kipf et al. 2018; Niculae et al. 2023; Cini et al. 2023b) can be incorporated
into the framework and constitute a possible extension for future works. Note
that the node embeddings $\mathbf{e}_{G}[n]$ need to be learned for every
location, which makes the model intrinsically transductive. Despite that, the
fine-tuning of new node embeddings can be done efficiently (Cini et al.,
2023a).

###  4.3 Graph-graph Network

In this section, we present the Graph-graph Network (GgNet), our proposed
architecture relying on the nested graph representations provided by $G$ and
$g$. The model is composed of several blocks, as detailed below and depicted
in Fig. 3.

##### Input encoder

To account for specifics along the temporal, spatial and channel dimensions,
we obtain a hidden representation for each time step $t$, location $n$, and
channel $d$ given the input data $\mathcal{X}$ and the learned node embeddings
$\mathbf{E}_{G}$, $\mathbf{E}_{g}$ similarly to Marisca et al. (2022)

$\mathbf{h}_{t}^{d}[n]=m_{t}^{d}[n]\cdot\text{MLP}_{\text{enc,}1}(\mathbf{x}_{t}^{d}[n])+\text{MLP}_{\text{enc,}2}(\mathbf{e}_{G}[n],\mathbf{e}_{g}[d])$
(3)

where $\mathbf{h}_{t}^{d}[n]\in\mathbb{R}^{H}$ are the resulting encoded
representations, while MLP${}_{\text{enc,}1}$ and MLP${}_{\text{enc,}2}$ are
two distinct multilayer perceptions. Note that representations are obtained by
exploiting observations wherever available and relying exclusively on node
embeddings for the targets.

##### Temporal processing

Figure 3: Overview of GgNet. Temporal convolutions (blue) encode temporal
patterns, $G$-convolutions (green) propagate information across the inter
location graph and model dependencies between locations; $g$-convolutions
(purple) propagate information across the intra-location graph and model
dependencies between channels. Starred modules refer to channel-wise
operations.

Stacks of Temporal Convolutions (TC) (Bai et al., 2018) are used to encode
temporal information for each location and channel as

$\vec{\mathbf{H}}^{d}[n]=\textnormal{TC}^{d}\left(\vec{\mathbf{H}}^{d}[n],\vec{\mathbf{m}}^{d}[n],\textbf{e}_{G}[n]\right)$
(4)

where $\vec{\mathbf{H}}^{d}[n]\in\mathbb{R}^{T\times H}$. As channels may
present largely heterogeneous temporal features, we use a different set of
weights (filters) for processing data in each channel. Temporal convolutions
are highly parallelizable making them preferable, in this setting, over, e.g.,
recurrent neural networks. Furthermore, we symmetrically pad the input
sequence to avoid duplicating the architecture to encode forward and backward
dynamics and exploit exponentially increasing dilation rates to capture both
short and long-range temporal dependencies.

##### Inter-location (spatial) processing

Spatial information is propagated across locations (T$\to$Y) by means of Graph
Convolutions (GC) (Bronstein et al., 2021) over the inter-location graph
($G$-convolution in Fig. 3):

$\mathbf{H}_{t}^{d}=\text{GC}_{(G)}\left(\mathbf{H}_{t}^{d},\mathbf{m}_{t}^{d},\mathbf{E}_{G},\mathbf{E}_{g};\mathbf{A}_{G}\right)$
(5)

where $\mathbf{H}_{t}^{d}\in\mathbb{R}^{N\times H}$ and the same convolution
is performed independently on each channel and each time step. Specifics of
each channel are accounted for via the conditioning on the embeddings in
$\mathbf{E}_{g}$. Note that G-convolution is a synchronous operation that
assumes that the time frames of each time series are aligned.

##### Intra-location (channel-wise) processing

As the next step, information is propagated across channels (C$\to$Y) by
performing GCs over the intra-location graph ($g$-convolution in Fig. 3):

$\mathbf{H}_{t}[n]=\text{GC}_{(g)}\left(\mathbf{H}_{t}[n],\mathbf{m}_{t}[n],\mathbf{E}_{G},\mathbf{E}_{g};\mathbf{A}_{g}\right)$
(6)

for all time steps and locations, and where
$\mathbf{H}_{t}[n]\in\mathbb{R}^{D\times H}$. $g$-convolutions allow for
inferring observations at unavailable channels by modeling dependencies among
targets and covariates. Note that weights are shared across locations, as
local dynamics are accounted for by conditioning on node embeddings
$\mathbf{E}_{G}$.

Multiple $T$-, $G$-, and $g$-convolution blocks can be stacked, allowing for
designing deep architectures. We refer to Appendix E.2 for an extensive
ablation study assessing the impact of each component.

##### Readout

Finally, a readout composed of $d$ MLPs maps representations to predictions at
all target locations and time steps as

$\hat{\mathbf{x}}_{t}^{d}[n]=\text{MLP}_{\text{dec}}^{\,d}(\mathbf{h}_{t}^{d}[n],\,m_{t}^{d}[n],\,\mathbf{e}_{G}[n],\,\mathbf{e}_{g}[d]).$
(7)

Note that the GgNet implementation considered here is just one of many
possible instantiations of the framework. The model can be extended, e.g., by
including more sophisticated graph processing and learning modules.

###  4.4 Model training

Figure 4: Training and evaluation splits for MVS.

For a generic MVS task, the reconstruction error can be defined as

$\mathcal{L}(\hat{\mathcal{X}},\mathcal{X},\mathcal{M})=\frac{\sum_{n}\sum_{d}\bar{m}^{d}[n]\cdot\ell(\hat{\vec{\textbf{x}}}^{\,d}[n],\vec{\textbf{x}}^{\,d}[n])}{\sum_{n}\sum_{d}\bar{m}^{d}[n]}$
(8)

where $\hat{\vec{\textbf{x}}}^{\,d}[n]$ and $\vec{\textbf{x}}^{\,d}[n]$ are,
respectively, the predicted and true values, and $\ell(\cdot,\cdot)$ is a loss
function between time series, e.g., the mean absolute error or the mean
squared error. The model is trained and evaluated by simulating the presence
of missing channels partitioning the available $N\times D$ channels into
training, validation, and test channels (Fig. 4), simulating a transductive
learning setting. To bias the training toward the virtual sensing task, we
mask out a fraction of the available channels for each training batch and
train the model to reconstruct the input data, giving a higher weight to the
masked out observations. We have found this strategy to be crucial for models
like GgNet that do not rely on forecasting as a surrogate objective. Finally,
we mask out an additional small fraction of data points, at random, in the
training data to enforce robustness to random missing values (Cini et al.,
2022; Du et al., 2023).

##  5 Experimental evaluation

In this section, we evaluate the proposed method in several MVS tasks.
Experiments are conducted across three different datasets, allowing for the
investigation of settings characterized by (i) different temporal resolutions
and (ii) different degrees of sparsity, as well as (iii) assessing the
limitations of the different approaches. Alongside GgNet, we provide results
for the following baselines: KNN, i.e., a simple baseline averaging the
nearest geographical neighbors, BRITS (Cao et al., 2018) as representative of
autoregressive RNN global models111With the term “global” we refer to models
whose parameters are shared in all locations, in contrast to “local” ones that
have location-specific parameters (Montero-Manso & Hyndman, 2021; Cini et al.,
2023a)., SAITS (Du et al., 2023) as representative of transformer-based global
models unaware of spatial relations, GRIN (Cini et al., 2022) as
representative of graph-based imputation models, leveraging spatial
correlations. Moreover, we consider a variation of the standard GRIN model,
named GRINm, modified to better deal with multivariate data at each location;
specifically, we add self-loops in GRIN’s Spatial Decoder for conditioning on
covariates at the same time step of the reconstruction target. To assess the
impact of some of the introduced design choices, we also compare $4$
progressively more advanced recurrent architectures: 1) RNN, a standard global
recurrent neural network; 2) $\textbf{RNN}_{\textbf{BiD.}}$, its bidirectional
extension; 3) $\textbf{RNN}_{\textbf{Emb.}}$, which, besides being
bidirectional, incorporates a local component by concatenating learnable node
embeddings to the input observations; 4) $\textbf{RNN}_{\textbf{G}}$, which
adds a convolution operation over an inter-location graph on top of the
previous architecture. For each model, we consider different hyperparameter
configurations, with increasing complexity; reported results correspond to the
best-performing model, given the reconstruction accuracy on the validation
set. Regarding the preprocessing steps, we follow the common practice of
standardizing the data, ensuring that, across locations and time, each of the
$D$ channels has zero mean and unit variance. Missing values, if present, are
masked out. Details about the baselines can be found in Appendix B, while the
complete experimental setup is described in Appendix C.

###  5.1 Missing-channel reconstruction: climatic data

In the first experiment, we address reconstructing missing-at-random (MaR)
channels in multivariate spatio-temporal datasets. We build two datasets with
different temporal resolutions of climatic variables from the NASA Langley
Research Center POWER Project. More information about the data and APIs is
available in Appendix D. To emulate a reduced spatial coverage, we collect
data in correspondence with the $235$ national capitals. For each location, we
collect multiple correlated variables with daily and hourly resolutions.
Within this setting, we aim to reconstruct MaR climatic variables at different
locations from the available data. Among all $N\times D$ channels, we use
$70\%$ for training, $10\%$ for validation, and the remaining $20\%$ for
testing.

##### Daily climatic data

Considering a daily resolution, we collect $10$ climatic variables (listed in
Appendix D) encompassing $30$ years, which results in $N=235,T=10958,D=10$.
Tab. 1 reports the Mean Relative Error (MRE) for each test channel, grouped by
variable and averaged across locations. All results are consistent across
different metrics as shown in Appendix E.8.

Table 1: Climatic dataset: channel-wise MRE (%) for daily data and average
accuracy for daily (D) and hourly (H) data. Results are averaged across
locations and 5 random seeds. The best-performing method is in bold, the
second-best is underlined.

Ch. | Temp. | Temp. | Temp. | Wind | Rel. | Prec. | Temp. | Clouds | Irr. | Irr. | Avg | Avg
---|---|---|---|---|---|---|---|---|---|---|---|---
| mean | range | max | speed | hum. |  | dew |  | short | long | (D) | (H)
KNN | $15.6_{\pm 2.5}$ | $42.1_{\pm 6.3}$ | $12.4_{\pm 1.5}$ | $33.7_{\pm 4.5}$ | $12.3_{\pm 0.4}$ | $102.8_{\pm 4.7}$ | $23.9_{\pm 2.9}$ | $34.0_{\pm 0.6}$ | $18.0_{\pm 0.5}$ | $6.0_{\pm 1.0}$ | $30.1_{\pm 1.0}$ | $36.2_{\pm 1.9}$
BRITS | $4.0_{\pm 0.6}$ | $17.5_{\pm 2.3}$ | $3.3_{\pm 0.5}$ | $32.1_{\pm 1.5}$ | $5.5_{\pm 0.8}$ | $80.5_{\pm 2.4}$ | $8.4_{\pm 1.1}$ | $29.9_{\pm 0.9}$ | $20.9_{\pm 1.0}$ | $3.4_{\pm 0.4}$ | $20.5_{\pm 0.5}$ | $26.3_{\pm 0.8}$
GRIN | $4.2_{\pm 0.7}$ | $18.6_{\pm 1.9}$ | $4.3_{\pm 0.5}$ | $30.7_{\pm 2.3}$ | $5.6_{\pm 0.6}$ | $77.3_{\pm 3.8}$ | $8.7_{\pm 1.0}$ | $29.1_{\pm 0.3}$ | $15.3_{\pm 0.7}$ | $3.9_{\pm 0.3}$ | $19.8_{\pm 0.6}$ | $23.3_{\pm 0.7}$
GRINm | $2.9_{\pm 0.9}$ | $12.9_{\pm 1.3}$ | $2.8_{\pm 0.6}$ | $30.7_{\pm 3.9}$ | $3.6_{\pm 0.9}$ | $69.9_{\pm 1.6}$ | $6.4_{\pm 1.1}$ | $\underline{20.5}_{\pm 0.9}$ | $\underline{11.7}_{\pm 0.6}$ | $3.6_{\pm 0.4}$ | $16.5_{\pm 0.6}$ | $22.7_{\pm 0.4}$
SAITS | $\underline{2.4}_{\pm 0.4}$ | $\underline{11.2}_{\pm 1.3}$ | $\underline{2.3}_{\pm 0.4}$ | $\underline{26.9}_{\pm 1.0}$ | $\underline{3.3}_{\pm 1.0}$ | $\underline{66.0}_{\pm 2.5}$ | $\underline{5.0}_{\pm 1.0}$ | $20.6_{\pm 0.7}$ | $14.2_{\pm 0.9}$ | $\textbf{2.8}_{\pm 0.3}$ | $\underline{15.5}_{\pm 0.5}$ | $\underline{22.2}_{\pm 0.5}$
GgNet | $\textbf{2.1}_{\pm 0.4}$ | $\textbf{9.6}_{\pm 0.7}$ | $\textbf{2.0}_{\pm 0.2}$ | $\textbf{23.9}_{\pm 2.2}$ | $\textbf{2.7}_{\pm 0.7}$ | $\textbf{60.6}_{\pm 2.1}$ | $\textbf{4.2}_{\pm 0.8}$ | $\textbf{16.5}_{\pm 0.5}$ | $\textbf{9.2}_{\pm 0.8}$ | $\underline{2.9}_{\pm 0.3}$ | $\textbf{13.4}_{\pm 0.2}$ | $\textbf{20.4}_{\pm 0.6}$
% Imp. | 12.5% | 14.3% | 13.0% | 11.1% | 18.1% | 8.1% | 16.0% | 19.9% | 21.4% | -3.6% | 13.5% | 8.1%
RNN | $5.6_{\pm 0.7}$ | $21.9_{\pm 0.7}$ | $5.8_{\pm 0.5}$ | $32.7_{\pm 1.2}$ | $6.9_{\pm 0.7}$ | $83.8_{\pm 1.8}$ | $10.8_{\pm 1.1}$ | $35.5_{\pm 1.0}$ | $21.8_{\pm 0.3}$ | $3.8_{\pm 0.4}$ | $22.9_{\pm 0.2}$ | $28.8_{\pm 0.8}$
\+ BiD. | $4.2_{\pm 0.6}$ | $19.1_{\pm 1.0}$ | $4.7_{\pm 0.4}$ | $31.1_{\pm 0.8}$ | $5.9_{\pm 0.6}$ | $78.0_{\pm 2.6}$ | $8.6_{\pm 1.0}$ | $31.4_{\pm 0.6}$ | $19.0_{\pm 0.4}$ | $3.6_{\pm 0.4}$ | $20.6_{\pm 0.3}$ | $25.7_{\pm 0.5}$
\+ Emb. | $4.1_{\pm 0.5}$ | $17.5_{\pm 1.5}$ | $4.3_{\pm 0.4}$ | $30.2_{\pm 2.0}$ | $5.5_{\pm 0.6}$ | $77.1_{\pm 1.6}$ | $8.5_{\pm 1.0}$ | $30.2_{\pm 0.4}$ | $17.5_{\pm 0.7}$ | $3.6_{\pm 0.3}$ | $19.8_{\pm 0.5}$ | $26.2_{\pm 0.9}$
\+ G | $3.6_{\pm 0.3}$ | $16.5_{\pm 1.3}$ | $3.9_{\pm 0.4}$ | $27.6_{\pm 1.6}$ | $5.2_{\pm 0.4}$ | $74.7_{\pm 1.8}$ | $7.3_{\pm 0.9}$ | $27.7_{\pm 0.5}$ | $14.0_{\pm 0.8}$ | $3.6_{\pm 0.4}$ | $18.4_{\pm 0.5}$ | $23.2_{\pm 1.1}$

From the results, we first observe that GgNet outperforms models relying
exclusively on pre-defined spatial relations (GRIN), demonstrating superior
performance in handling sparse settings. Notably, GRIN underperforms
especially in correspondence with remote locations, where no neighboring
sensor is available (Appendix E.3). The GRINm variant improves over the
results of GRIN, but nonetheless underperforms w.r.t. GgNet.

Figure 5: (Left) Reconstruction of wind speed (poorly correlated with the
other channels) in Road Town (which can exploit observations at nearby
Caribbean capitals). (Right) t-SNE representation of the node embeddings
learned by GgNet; colours refer to the Köppen-Geiger climate classification.

SAITS performs reasonably well on this dataset, ranking second in most of the
setups. In particular, SAITS can account for intra-location covariates, but
cannot take advantage of inter-location relations. This is particularly
evident for missing channels that are hard to infer from the other variables
at the same location (e.g., wind speed, precipitations, clouds and shortwave
irradiation), but can be reconstructed from spatially close observations, as
in the example shown in Fig. 5 (left). Tab. 1’s lower section reports results
for the RNN model variants. The notable performance improvement obtained by
the bidirectional variant demonstrates the importance of accounting for both
past and future values. Node embeddings yield modest yet consistent
improvements. Additionally, the introduction of the inter-location graph G
results in a substantial performance gain, supporting the importance of inter-
location relations. Alongside MVS results, we provide in Fig. 5 (right) a
t-SNE visualization of the node embeddings (EG, Sec. 4.2) as learned by GgNet
w.r.t. each location. Each embedding is coloured according to the Köppen-
Geiger climate classification (Peel et al., 2007) that divides the globe into
different climates, with similar colours corresponding to similar climates.
Note that the model does not have access, at any point, to either these
climate labels or to the geographical coordinates of any location. The
correspondence between such labels and the formed clusters serves as a
qualitative assessment of the capability of GgNet to capture similarities in
the observed dynamics reflecting the climate at each location. Appendix E
reports several additional analyses.

##### Hourly climatic data

Then, we address the same MaR task in the case of $7$ climatic variables
sampled with hourly resolution. Data cover one year, which results in
$N=235,T=8760,D=7$. Differently from daily data, hourly data are more
difficult to predict, being characterized by high-frequency components both
temporally and spatially. As a consequence, exploiting inter-location
dependencies can help in providing better conditioning and hence lead to more
accurate predictions. Average MRE is reported in the last column of Tab. 1,
while channel-wise results are shown in Appendix E.8. Graph-based models
effectively exploit the additional information to obtain a relatively good
reconstruction accuracy. In GgNet, accounting for the interplay between
spatial aggregation (T$\rightarrow$Y) and global modeling of the relations
across channels (C$\rightarrow$Y) arguably results in superior average
performance.

###  5.2 Extending spatial coverage: Photovoltaic energy production

We consider the problem of predicting a specific target variable (whose direct
observation is expensive) at new target locations by exploiting a set of
available covariates. This constitutes an interesting use case for our
multivariate virtual sensing task. For example, renewable energy production or
material degradation can be inferred from widely available satellite climatic
data (De Felice et al., 2022). In particular, we consider as target variable
the daily photovoltaic power (PV) production simulated over continental North
America by Hu et al. (2022). To simulate a sparse scenario, we randomly sample
a few locations from the provided dense grid. As covariates, we use the same
$10$ climatic variables from Sec. 5.1 collected in correspondence with the
sampled locations. For both target and covariates, we consider $1$ full year
of daily observations, which results in $T=365,D=11$. By varying the number of
locations between $10$ and $200$, we aim to infer the PV production at new
target locations from the climate at such locations. We use 70% of the sensors
for training, 10% for validation, and 20% for testing. Tab. 2 reports the
reconstruction accuracy in terms of the Mean Absolute Error (MAE).

Table 2: Photovoltaic dataset: MAE (W*h). For each N, results are averaged
across 5 sampling of the locations and 5 different runs for each set (25 runs
in total). The best-performing method is in bold, the second-best is
underlined.

N${}_{\text{locations}}$ | 10 | 20 | 30 | 50 | 70 | 100 | 150 | 200
---|---|---|---|---|---|---|---|---
GRIN | $581\pm 170$ | $536\pm 52$ | $533\pm 49$ | $514\pm 36$ | $531\pm 34$ | $524\pm 28$ | $523\pm 40$ | $528\pm 43$
GRINm | $497\pm 117$ | $458\pm 49$ | $\underline{446}\pm 42$ | $\underline{424}\pm 23$ | $\underline{428}\pm 21$ | $423\pm 13$ | $\underline{417}\pm 13$ | $418\pm 12$
SAITS | $\underline{466}\pm 68$ | $\textbf{451}\pm 46$ | $448\pm 37$ | $427\pm 29$ | $436\pm 24$ | $\underline{420}\pm 16$ | $\underline{417}\pm 14$ | $\underline{416}\pm 14$
GgNet | $\textbf{452}\pm 67$ | $\underline{455}\pm 39$ | $\textbf{439}\pm 32$ | $\textbf{415}\pm 26$ | $\textbf{410}\pm 26$ | $\textbf{391}\pm 18$ | $\textbf{367}\pm 15$ | $\textbf{357}\pm 16$
% Imp. | 3.0% | -0.9% | 1.6% | 2.1% | 4.2% | 6.9% | 12.0% | 14.2%

As more locations are considered, results show a progressive increment of the
relative improvement of GgNet over SAITS. This could be related to the
superior capability of GgNet to exploit spatial information.

##  6 Conclusion

We presented a novel framework for virtual sensing from sparse multivariate
spatio-temporal observations. In doing so, we introduced a novel methodology
to exploit dependencies between covariates and relations across locations. In
this context, we proposed GgNet, a graph deep learning architecture leveraging
a nested graph structure to account for such relations. The relational
structure is learned end-to-end to maximize reconstruction accuracy. Compared
to state-of-the-art, GgNet can achieve superior performance in settings with
poor sensor coverage, where other methods fail. As future directions, we
believe it would be interesting to explore the application of a similar
methodology to collections of asynchronous time series, typically found in
healthcare applications. Finally, GgNet is limited, in its current form, to
transductive learning settings and its extension to inductive learning would
broaden the applicability of the proposed method.

#### Reproducibility statement

Python code to reproduce the experiments is available online at
https://github.com/gdefe/ggnet-virtual-sensing. A script to download climatic
data from the database in correspondence with the world capitals is included
in the provided code. PV datasets at different spatial densities, together
with the corresponding climatic variables, are provided alongside the code.
Complete climatic and photovoltaic datasets are publicly available online
through the references provided in Appendix D. The sampling of locations
(photovoltaic experiment) and simulation of missing channels (both climatic
and photovoltaic experiment) are controlled by a fixed seed for
reproducibility purposes. Hyperparameters for all models are provided in
Appendix C.

#### Acknowledgments

G.D.F. acknowledges the Beckers Group for funding this research. V.V.G. thanks
Leverhulme Trust for support via the Leverhulme Research Centre for Functional
Materials Design and EPSRC under grant number EP/V026887. This research was
partly funded by the Swiss National Science Foundation under grant 204061:
_High-Order Relations and Dynamics in Graph Neural Networks_. Climate data was
obtained from the National Aeronautics and Space Administration (NASA) Langley
Research Center (LaRC) Prediction of Worldwide Energy Resource (POWER) Project
funded through the NASA Earth Science/Applied Science Program. The authors
wish to thank the anonymous reviewers whose comments considerably improved the
final version of the paper.

## References

  * Appleby et al. (2020) Gabriel Appleby, Linfeng Liu, and Li-Ping Liu.  Kriging convolutional networks.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34(04):3187–3194, 2020.
  * Bacciu et al. (2020) Davide Bacciu, Federico Errica, Alessio Micheli, and Marco Podda.  A gentle introduction to deep learning for graphs.  _Neural Networks_ , 129:203–221, 2020.
  * Bahadori et al. (2014) Mohammad Taha Bahadori, Qi Rose Yu, and Yan Liu.  Fast multivariate spatio-temporal analysis via low rank tensor learning.  _Advances in neural information processing systems_ , 27, 2014.
  * Bai et al. (2020) Lei Bai, Lina Yao, Can Li, Xianzhi Wang, and Can Wang.  Adaptive graph convolutional recurrent network for traffic forecasting.  _Advances in neural information processing systems_ , 33:17804–17815, 2020.
  * Bai et al. (2018) Shaojie Bai, J Zico Kolter, and Vladlen Koltun.  An empirical evaluation of generic convolutional and recurrent networks for sequence modeling.  _arXiv preprint arXiv:1803.01271_ , 2018.
  * Biewald (2020) Lukas Biewald.  Experiment tracking with weights and biases, 2020.  URL https://www.wandb.com/.  Software available from wandb.com.
  * Bonilla et al. (2007) Edwin V Bonilla, Kian Chai, and Christopher Williams.  Multi-task gaussian process prediction.  _Advances in neural information processing systems_ , 20, 2007.
  * Bronstein et al. (2021) Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković.  Geometric deep learning: grids, groups, graphs, geodesics, and gauges.  _arXiv preprint arXiv:2104.13478_ , 2021.
  * Brunello et al. (2021) Andrea Brunello, Andrea Urgolo, Federico Pittino, András Montvay, and Angelo Montanari.  Virtual sensing and sensors selection for efficient temperature monitoring in indoor environments.  _Sensors_ , 21(8):2728, 2021.
  * Cao et al. (2018) Wei Cao, Dong Wang, Jian Li, Hao Zhou, Lei Li, and Yitan Li.  BRITS: Bidirectional recurrent imputation for time series.  _Advances in neural information processing systems_ , 31, 2018.
  * Che et al. (2018) Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu.  Recurrent neural networks for multivariate time series with missing values.  _Scientific reports_ , 8(1):6085, 2018.
  * Chen et al. (2022) Yakun Chen, Zihao Li, Chao Yang, Xianzhi Wang, Guodong Long, and Guandong Xu.  Adaptive graph recurrent network for multivariate time series imputation.  In _International Conference on Neural Information Processing_ , pp. 64–73. Springer, 2022.
  * Cini & Marisca (2022) Andrea Cini and Ivan Marisca.  Torch Spatiotemporal, 3 2022.  URL https://github.com/TorchSpatiotemporal/tsl.
  * Cini et al. (2022) Andrea Cini, Ivan Marisca, and Cesare Alippi.  Filling the g_ap_s: multivariate time series imputation by graph neural networks.  In _International Conference on Learning Representations_ , 2022.
  * Cini et al. (2023a) Andrea Cini, Ivan Marisca, Daniele Zambon, and Cesare Alippi.  Taming local effects in graph-based spatiotemporal forecasting.  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023a.
  * Cini et al. (2023b) Andrea Cini, Daniele Zambon, and Cesare Alippi.  Sparse graph learning from spatiotemporal time series.  _Journal of Machine Learning Research_ , 24(242):1–36, 2023b.
  * Cressie & Wikle (2015) Noel Cressie and Christopher K Wikle.  _Statistics for spatio-temporal data_.  John Wiley & Sons, 2015.
  * De Felice et al. (2022) Giovanni De Felice, Vladimir Gusev, John Y Goulermas, Michael Gaultois, Matthew Rosseinsky, and Catherine Vincent Gauvin.  Spatio-temporal weathering predictions in the sparse data regime with gaussian processes.  In _NeurIPS 2022 AI for Science: Progress and Promises_ , 2022.
  * De Felice et al. (2023) Giovanni De Felice, John Y Goulermas, and Vladimir Gusev.  Time series kernels based on nonlinear vector autoregressive delay embeddings.  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023.
  * Du et al. (2023) Wenjie Du, David Côté, and Yan Liu.  SAITS: self-attention-based imputation for time series.  _Expert Systems with Applications_ , 219:119619, 2023.
  * Gao & Ribeiro (2022) Jianfei Gao and Bruno Ribeiro.  On the equivalence between temporal and static equivariant graph representations.  In _International Conference on Machine Learning_ , pp. 7052–7076. PMLR, 2022.
  * Goovaerts (1998) Pierre Goovaerts.  Ordinary cokriging revisited.  _Mathematical Geology_ , 30:21–42, 1998.
  * Hu et al. (2023) Junfeng Hu, Yuxuan Liang, Zhencheng Fan, Hongyang Chen, Yu Zheng, and Roger Zimmermann.  Graph neural processes for spatio-temporal extrapolation.  In _Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , pp. 752–763, 2023.
  * Hu et al. (2022) Weiming Hu, Guido Cervone, Andre Merzky, Matteo Turilli, and Shantenu Jha.  A new hourly dataset for photovoltaic energy production for the continental usa.  _Data in Brief_ , 40:107824, 2022.
  * Jayawardene & Venayagamoorthy (2016) Iroshani Jayawardene and Ganesh K Venayagamoorthy.  Spatial predictions of solar irradiance for photovoltaic plants.  In _2016 IEEE 43rd Photovoltaic Specialists Conference (PVSC)_ , pp. 0267–0272. IEEE, 2016.
  * Jin et al. (2023) Ming Jin, Huan Yee Koh, Qingsong Wen, Daniele Zambon, Cesare Alippi, Geoffrey I Webb, Irwin King, and Shirui Pan.  A survey on graph neural networks for time series: forecasting, classification, imputation, and anomaly detection.  _arXiv preprint arXiv:2307.03759_ , 2023.
  * Kingma & Ba (2015) Diederik P Kingma and Jimmy Ba.  Adam: a method for stochastic optimization.  In _International Conference on Learning Representations_ , 2015.
  * Kipf et al. (2018) Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel.  Neural relational inference for interacting systems.  In _International conference on machine learning_ , pp. 2688–2697. PMLR, 2018.
  * Koenker & Bassett Jr (1978) Roger Koenker and Gilbert Bassett Jr.  Regression quantiles.  _Econometrica: journal of the Econometric Society_ , pp. 33–50, 1978.
  * Kong et al. (2023) Xiangjie Kong, Wenfeng Zhou, Guojiang Shen, Wenyi Zhang, Nali Liu, and Yao Yang.  Dynamic graph convolutional recurrent imputation network for spatiotemporal traffic missing data.  _Knowledge-Based Systems_ , 261:110188, 2023.
  * Kuppannagari et al. (2021) Sanmukh R Kuppannagari, Yao Fu, Chung Ming Chueng, and Viktor K Prasanna.  Spatio-temporal missing data imputation for smart power grids.  In _Proceedings of the Twelfth ACM International Conference on Future Energy Systems_ , pp. 458–465, 2021.
  * Lei et al. (2022) Mengying Lei, Aurelie Labbe, Yuankai Wu, and Lijun Sun.  Bayesian kernelized matrix factorization for spatiotemporal traffic data imputation and kriging.  _IEEE Transactions on Intelligent Transportation Systems_ , 23(10):18962–18974, 2022.
  * Lepot et al. (2017) Mathieu Lepot, Jean-Baptiste Aubin, and François HLR Clemens.  Interpolation in time series: an introductive overview of existing methods, their performance criteria and uncertainty assessment.  _Water_ , 9(10):796, 2017.
  * Li et al. (2018) Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu.  Diffusion convolutional recurrent neural network: data-driven traffic forecasting.  In _International Conference on Learning Representations_ , 2018.
  * Little & Rubin (2019) Roderick JA Little and Donald B Rubin.  _Statistical analysis with missing data_ , volume 793.  John Wiley & Sons, 2019.
  * Liu et al. (2009) Lichuan Liu, Sen M Kuo, and MengChu Zhou.  Virtual sensing techniques and their applications.  In _2009 International Conference on Networking, Sensing and Control_ , pp. 31–36. IEEE, 2009.
  * Luo et al. (2019) Yonghong Luo, Ying Zhang, Xiangrui Cai, and Xiaojie Yuan.  E2GAN: end-to-end generative adversarial network for multivariate time series imputation.  In _Proceedings of the 28th international joint conference on artificial intelligence_ , pp. 3094–3100, 2019.
  * Luttinen & Ilin (2012) Jaakko Luttinen and Alexander Ilin.  Efficient gaussian process inference for short-scale spatio-temporal modeling.  In _Artificial Intelligence and Statistics_ , pp. 741–750. PMLR, 2012.
  * Ma et al. (2019) Jiawei Ma, Zheng Shou, Alireza Zareian, Hassan Mansour, Anthony Vetro, and Shih-Fu Chang.  CDSA: cross-dimensional self-attention for multivariate, geo-tagged time series imputation.  _arXiv preprint arXiv:1905.09904_ , 2019.
  * Marisca et al. (2022) Ivan Marisca, Andrea Cini, and Cesare Alippi.  Learning to reconstruct missing data from spatiotemporal graphs with sparse observations.  _Advances in Neural Information Processing Systems_ , 35:32069–32082, 2022.
  * Mikalsen et al. (2018) Karl Øyvind Mikalsen, Filippo Maria Bianchi, Cristina Soguero-Ruiz, and Robert Jenssen.  Time series cluster kernel for learning similarities between multivariate time series with missing data.  _Pattern Recognition_ , 76:569–581, 2018.
  * Montero-Manso & Hyndman (2021) Pablo Montero-Manso and Rob J Hyndman.  Principles and algorithms for forecasting groups of time series: Locality and globality.  _International Journal of Forecasting_ , 37(4):1632–1653, 2021.
  * Niculae et al. (2023) Vlad Niculae, Caio F Corro, Nikita Nangia, Tsvetomila Mihaylova, and André FT Martins.  Discrete latent structure in neural networks.  _arXiv preprint arXiv:2301.07473_ , 2023.
  * Paepae et al. (2021) Thulane Paepae, Pitshou N Bokoro, and Kyandoghere Kyamakya.  From fully physical to virtual sensing for water quality assessment: a comprehensive review of the relevant state-of-the-art.  _Sensors_ , 21(21):6971, 2021.
  * Peel et al. (2007) Murray C Peel, Brian L Finlayson, and Thomas A McMahon.  Updated world map of the köppen-geiger climate classification.  _Hydrology and earth system sciences_ , 11(5):1633–1644, 2007.
  * Roth & Liebig (2022) Andreas Roth and Thomas Liebig.  Forecasting unobserved node states with spatio-temporal graph neural networks.  In _2022 IEEE International Conference on Data Mining Workshops (ICDMW)_ , pp. 740–747. IEEE, 2022.
  * Scarselli et al. (2008) Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.  The graph neural network model.  _IEEE transactions on neural networks_ , 20(1):61–80, 2008.
  * Shukla & Marlin (2021) Satya Narayan Shukla and Benjamin Marlin.  Multi-time attention networks for irregularly sampled time series.  In _International Conference on Learning Representations_ , 2021.
  * Shuman et al. (2013) David I Shuman, Sunil K Narang, Pascal Frossard, Antonio Ortega, and Pierre Vandergheynst.  The emerging field of signal processing on graphs: extending high-dimensional data analysis to networks and other irregular domains.  _IEEE signal processing magazine_ , 30(3):83–98, 2013.
  * Stein (1999) Michael L Stein.  _Interpolation of spatial data: some theory for kriging_.  Springer Science & Business Media, 1999.
  * Takeuchi et al. (2017) Koh Takeuchi, Hisashi Kashima, and Naonori Ueda.  Autoregressive tensor factorization for spatio-temporal predictions.  In _2017 IEEE international conference on data mining (ICDM)_ , pp. 1105–1110. IEEE, 2017.
  * Tashiro et al. (2021) Yusuke Tashiro, Jiaming Song, Yang Song, and Stefano Ermon.  CSDI: conditional score-based diffusion models for probabilistic time series imputation.  _Advances in Neural Information Processing Systems_ , 34:24804–24816, 2021.
  * Van Rossum et al. (1995) Guido Van Rossum, Fred L Drake, et al.  _Python reference manual_ , volume 111.  Centrum voor Wiskunde en Informatica Amsterdam, 1995.
  * Wang et al. (2020) Senzhang Wang, Jiannong Cao, and S Yu Philip.  Deep learning for spatio-temporal data mining: a survey.  _IEEE transactions on knowledge and data engineering_ , 34(8):3681–3700, 2020.
  * Wu et al. (2022) Xuesong Wu, Mengyun Xu, Jie Fang, and Xiongwei Wu.  A multi-attention tensor completion network for spatiotemporal traffic data imputation.  _IEEE Internet of Things Journal_ , 9(20):20203–20213, 2022.
  * Wu et al. (2021a) Yuankai Wu, Dingyi Zhuang, Aurelie Labbe, and Lijun Sun.  Inductive graph neural networks for spatiotemporal kriging.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 35(5):4478–4485, 2021a.
  * Wu et al. (2021b) Yuankai Wu, Dingyi Zhuang, Mengying Lei, Aurelie Labbe, and Lijun Sun.  Spatial aggregation and temporal convolution networks for real-time kriging.  _arXiv preprint arXiv:2109.12144_ , 2021b.
  * Ye et al. (2021) Yongchao Ye, Shiyao Zhang, and James JQ Yu.  Spatial-temporal traffic data imputation via graph attention convolutional network.  In _International Conference on Artificial Neural Networks_ , pp. 241–252. Springer, 2021.
  * Yoon et al. (2018a) Jinsung Yoon, James Jordon, and Mihaela Schaar.  GAIN: Missing data imputation using generative adversarial nets.  In _International conference on machine learning_ , pp. 5689–5698. PMLR, 2018a.
  * Yoon et al. (2018b) Jinsung Yoon, William R Zame, and Mihaela van der Schaar.  Estimating missing data in temporal data streams using multi-directional recurrent neural networks.  _IEEE Transactions on Biomedical Engineering_ , 66(5):1477–1490, 2018b.
  * Zhang et al. (2022) Xiyuan Zhang, Ranak Roy Chowdhury, Jingbo Shang, Rajesh Gupta, and Dezhi Hong.  ESC-GAN: extending spatial coverage of physical sensors.  In _Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining_ , pp. 1347–1356, 2022.
  * Zheng et al. (2023) Chuanpan Zheng, Xiaoliang Fan, Cheng Wang, Jianzhong Qi, Chaochao Chen, and Longbiao Chen.  INCREASE: inductive graph representation learning for spatio-temporal kriging.  In _Proceedings of the ACM Web Conference 2023_ , pp. 673–683, 2023\.

## Appendix

##  Appendix A Related methods

Methods that can perform virtual sensing can be broadly categorized into
spatio-temporal kriging methods and generic MTS imputation approaches. Kriging
methods exploit the density of the network and infer the missing sensors by
interpolating neighboring observations. Among those, ordinary kriging (Stein,
1999) is a linear interpolation procedure that assigns optimal weights for
neighboring data points based on a variogram model. Similarly, cokriging
(Goovaerts, 1998; Cressie & Wikle, 2015) also exploits correlation with
auxiliary covariate information (such as elevation, soil type, or temperature)
and estimates them simultaneously with the variable of interest. However,
these are linear spatial interpolation methods, therefore, they fail to
capture any temporal patterns, as well as nonlinear dependencies within the
data. Multi-task Gaussian Processes leverage flexible kernel structures to
capture spatio-temporal correlations (Bonilla et al., 2007; Luttinen & Ilin,
2012). In addressing scalability issues when dealing with extensive datasets,
matrix/tensor completion methods capture global consistency in the data by
imposing a low-rank assumption, as well as local consistency by diverse key
regularization structures (Bahadori et al., 2014; Takeuchi et al., 2017; Lei
et al., 2022; Wu et al., 2022). More related to our approach, imputation
methods based on GNNs have demonstrated great effectiveness in modeling
spatial dependencies between time series (Ye et al., 2021; Kuppannagari et
al., 2021; Roth & Liebig, 2022; Kong et al., 2023; Jin et al., 2023). Notably,
GRIN (Cini et al., 2022) accounts for both temporal and spatial correlation
within the data by combining autoregressive modeling with message-passing
operations; SPIN (Marisca et al., 2022) introduces a spatiotemporal attention
mechanism to weigh discrete points in time and space; IGNNK (Wu et al., 2021a)
tackle inductive virtual sensing by specializing spatio-temporal graph neural
networks, later extended in SATCN (Wu et al., 2021b) by allowing multiple
transformations of the spatial information to contribute to the adjacency
matrix; INCREASE (Zheng et al., 2023) also considers different heterogeneous
spatial relations. However, they all assume that the spatial information
characterizing the placement of each sample is available to the model; either
as a given graph structure or as geographical information to compute one. In
sparse settings, or when such information is missing, the main limitation of
these methods is deriving an appropriate adjacency matrix. This problem is
even exacerbated in the multivariate case. As evidence, MTS similarity is a
challenging problem on its own (Mikalsen et al., 2018; De Felice et al.,
2023). Adaptive methods (Bai et al., 2020), such as AGRN (Chen et al., 2022),
address this by learning the adjacency matrix together with the network
weights. Despite this, existing GNN-based approaches only target dense and
univariate sensor networks and do not explicitly model the relations between
channels in individual locations.

On the other hand, recent advancements in generic MTS imputation have led to
models that can be adapted to perform virtual sensing. Among these, the most
widely adopted are deep autoregressive methods based on recurrent neural
networks (Che et al., 2018; Yoon et al., 2018b), e.g. BRITS (Cao et al., 2018)
is a bidirectional model that also takes advantage of the relationships
between dimensions. Moreover, imputation can be performed by modeling the
underlying data distribution, generally employing generative adversarial
neural networks (Yoon et al., 2018a; Luo et al., 2019). In recent times,
imputation techniques have been proposed that rely on attention mechanisms (Ma
et al., 2019; Shukla & Marlin, 2021; Du et al., 2023) and diffusion models
(Tashiro et al., 2021). However, these models cannot take spatial correlations
directly into account. Indeed, these methods learn a single model shared
across locations and take into account a single MTS at a time. One could still
exploit such approaches by considering all the available data as a single MTS,
but the resulting model would not be appropriate for virtual sensing, as the
same variables at different locations would be modeled as entirely different
channels.

##  Appendix B Baseline details

##### KNN:

imputes missing test channels by averaging the corresponding train channels in
the $k$-nearest geographical neighbors. The number of neighbors
$k\in\\{1,2,3,5,10\\}$ is selected based on the best performance on the
validation set. This baseline has access to the geographical coordinates, with
which neighbors are identified.

##### RNN:

global model where data at different locations
$\textbf{X}[n]\in\mathbb{R}^{T\times D}$ are considered as different training
instances to train a shared model. Specifically, this model imputes missing
values from a linear transformation of the hidden state of a global Recurrent
Neural Network (RNN). Formally, the following operations are sequentially
iterated for every time step $t$ in the considered temporal window:

$\hat{\textbf{x}}_{t}=\textbf{W}_{out}\,\textbf{h}_{t}+\textbf{b}$ (9)
$\textbf{x}_{t}=\textbf{m}_{t}\odot\textbf{x}_{t}+\bar{\textbf{m}}_{t}\odot\hat{\textbf{x}}_{t}$
(10)
$\textbf{h}_{t}=\textnormal{GRU}(\,[\,\textbf{x}_{t}\,||\,\textbf{m}_{t}\,],\,\textbf{h}_{t-1}\,)$
(11)

where the notation $[n]$ is implicit for all inputs, masks and hidden states.

#####  ${\textbf{RNN}}_{\textbf{BiD}}$:

this model replicates the previous RNN model both in the forward direction and
in the backward direction. Imputations performed by Eq. 9 are discarded and
the final imputation is obtained with a separate MLP acting on the aggregated
hidden representations:

$\hat{\textbf{x}}_{t}=\textnormal{MLP}(\,[\,\textbf{h}^{fwd}_{t}\,||\,\textbf{h}^{bwd}_{t}\,]\,)$
(12)

#####  ${\textbf{RNN}}_{\textbf{emb}}$:

as global models, it is not possible for the two previous models to specialize
the modeling on the specifics of each location. In this extension, we
incorporate local components into RNN${}_{\text{BiD}}$ by concatenating the
node embeddings e[n] to all input time steps before the state update and
readout. Formally, the model sequentially iterates, on both directions, for
every time step $t$, Eq. 9, Eq. 10 and:

$\textbf{h}_{t}[n]=\textnormal{GRU}(\,[\,\textbf{x}_{t}[n]\,||\,\textbf{m}_{t}[n]\,||\,\textbf{e}[n]\,],\,\textbf{h}_{t-1}[n]\,)$
(13)

where we have here explicit the notation $[n]$. Finally, as before, we
aggregate hidden representations to output the final result:

$\hat{\textbf{x}}_{t}[n]=\textnormal{MLP}(\,[\,\textbf{h}^{fwd}_{t}[n]\,||\,\textbf{h}^{bwd}_{t}[n]\,||\,\textbf{e}[n]\,]\,)$
(14)

#####  ${\textbf{RNN}}_{\textbf{G}}$:

we further build here on the previous model by introducing synchronous
message-passing operations between locations. The resulting architecture
mimics time-then-space models that are typical to the spatio-temporal graph
processing literature (Gao & Ribeiro, 2022). First, inputs are processed by
the RNN${}_{\text{emb}}$, as above, up to Eq. 13; then, a GC layer performs a
convolution operation over the inter-location graph $G$:

$\textbf{H}^{\prime}_{t}=\textbf{A}_{G}\,\textbf{H}_{t}\,\Theta+\textbf{H}_{t}\,\Theta_{skip}+\textbf{b}$
(15)

where $\textbf{H}_{t}=[\,\textbf{H}^{fwd}_{t}\,||\,\textbf{H}^{bwd}_{t}\,]$
and $\textbf{A}_{G}$ is obtained from Eq. 2. Finally, RNN and GNN states are
concatenated and fed to an MLP readout to produce the final imputation:

$\hat{\textbf{x}}_{t}[n]=\textnormal{MLP}(\,[\,\textbf{h}^{fwd}_{t}[n]\,||\,\textbf{h}^{bwd}_{t}[n]\,||\,\textbf{h}^{\prime}_{t}[n]\,||\,\textbf{e}[n]\,]\,)$
(16)

##### BRITS:

global bidirectional RNN model for generic MTS imputation from Cao et al.
(2018) 222https://github.com/caow13/BRITS. The 2D form is recovered by
processing one MTS at a time, i.e., in practice, flattening the location
dimension along the batch dimension.

##### SAITS:

global self-attention-based MTS imputation model from Du et al. (2023)
333https://github.com/WenjieDu/SAITS. The 2D form is recovered by processing
one MTS at a time, i.e., in practice, flattening the location dimension along
the batch dimension.

##### GRIN:

graph-based model for spatio-temporal data imputation from Cini et al. (2022)
444https://github.com/Graph-Machine-Learning-Group/grin. This baseline has
access to the geographical coordinates, from which we compute the adjacency
matrix with a thresholded Gaussian kernel (Shuman et al., 2013):

$A_{ij}=A_{ji}=\begin{cases}\exp{\bigg{(}-\dfrac{dist(i,j)^{2}}{\sigma^{2}}\bigg{)}}&\text{if}\quad
dist(i,j)<\delta\\\ 0&\text{otherwise}\end{cases}$ (17)

with $dist(\cdot,\cdot)$ being the Haversine distance, $\sigma^{2}$ the
distances variance, and $\delta$ a custom threshold. For all datasets with a
number of nodes $N\geq 70$, we set $\delta$ to the required value to obtain an
average of 10 edges per node, e.g., $\delta(\text{climatic datasets})=3500$
km. For smaller datasets, we fix $\delta$ to the corresponding value at
$N=70$. The GRINm variant adds self-loops to the original Spatial Decoder of
GRIN.

##  Appendix C Detailed experimental setting

##### Shared settings

For GgNet and all baselines, we adopt the Adam optimizer (Kingma & Ba, 2015)
with a learning rate $lr=0.001$ paired with a cosine annealing learning rate
scheduler. All models are trained for a maximum of 500 epochs, with a 30
epochs patience for early stopping. All batches have a size set to 32 and
consider temporal windows of $t_{w}=24$ time steps. As for the settings that
are specific to virtual sensing: we set, for all methods, the probability of
randomly masking training points to $p_{withen-points}=0.05$, the probability
of randomly masking entire channels to $p_{withen-channels}=0.3$ and the
weight of masked points during loss calculation to $w_{whiten}=5$. Different
values $w_{whiten}\in\\{3,5,10\\}$ led to a negligible difference in
performance. A brief study on the robustness under changes in the $p_{withen-
channels}$ mask can be found in Appendix E.6. Experiments are conducted within
the Python (Van Rossum et al., 1995) library Torch Spatiotemporal (Cini &
Marisca, 2022), which is also used to implement all methods. Experiments are
tracked using Weight and Biases (Biewald, 2020). Code to reproduce the
experiments is available online.555https://github.com/gdefe/ggnet-virtual-
sensing

##### Hyperparameter settings on climatic dataset

In the climatic data, all models are evaluated under different hyperparameter
configurations, with increasing complexity. Separately for the daily and
hourly datasets, the best-performing configuration of parameters on the
validation set across the training epochs is used at testing. For GgNet, this
resulted in the same set of parameters for both datasets: hidden size
$d_{h}=128$ for all convolutions and MLPs; location node embedding size
$h_{e,G}=16$, channel node embedding size $h_{e,g}=8$; 2 blocks, each of which
is composed of 3 layers of temporal convolution (with filter length
$k_{\text{temp. filter}}=3$ centered at each time step $t$), 1 layer of
G-convolution and 1 layer of g-convolution; “elu” activation functions for the
encoding MLPs, all graph convolutions and decoding MLPs; “Tanh” activation
functions for the MLPs transforming the embeddings (Eq. 2). Hyperparameter
choices for all baselines are available within the configurations file in the
provided code. RNN variants are set to the same hidden size $h=64$ for
comparability.

##### Hyperparameter settings on photovoltaic dataset

In the photovoltaic dataset, to assess performance scaling with the number of
locations, we set a common hidden size of $h=64$ to all models and leave all
other parameters to the best-performing setting in the climatic daily
experiment. For all models, we also set a dropout value of $dr=0.1$ in all
experiments up to $N=70$, and $dr=0$ afterward.

###  C.1 Loss functions

##### GgNet loss function

As our loss function $\ell(\cdot,\cdot)$, we employ the sum of three quantile
losses, (also referred to as ’pinball loss’) (Koenker & Bassett Jr, 1978),
with quantile levels selected as follows: $q_{-\sigma}=0.159$, $q_{\mu}=0.5$,
and $q_{+\sigma}=0.841$. The median quantile aligns with the Mean Absolute
Error (MAE), while $q_{-\sigma}$ and $q_{+\sigma}$ allow for the concurrent
estimation of uncertainties. We believe this is particularly relevant to
virtual sensing, as not all channels in all locations can be equally
reconstructed from the available data. As in (Cini et al., 2023a), further
embedding regularizations can be considered in addition to the training
objective.

##### Baseline loss functions

RNN and ${\text{RNN}}_{\text{{BiD}}}$ use MAE loss;
${\text{RNN}}_{{\text{emb}}}$ and ${\text{RNN}}_{{G}}$ use the same three-
quantiles loss as GgNet; BRITS, SAITS and GRIN use same loss functions as the
authors provide in the respective code implementations.

##  Appendix D Datasets details

##### Climatic dataset:

We obtain daily and hourly climatic data from the POWER Project’s Daily and
Hourly 2.3.5 version on 2023/02/26. This online database contains data about
surface solar energy fluxes and other meteorological quantities obtained
through satellite systems and further reanalysis. For radiation data, spatial
resolution is 1.0° (latitude) by 1.0° (longitude); for meteorological data,
this is ${1/2}$° (latitude) by ${5/8}$° (longitude). For each point on this
global-scale grid, data are provided in time-series format with customizable
temporal resolution. Further information, together with data and API, is
available at the project website 666https://power.larc.nasa.gov/. For daily
data, we select the following 10 variables: mean temperature ($C$),
temperature range ($C$), maximum temperature ($C$), wind speed ($m/s$),
relative humidity ($\%$), precipitation ($mm/\text{day}$), dew/frost point
($C$), cloud amount ($\%$), all-sky surface shortwave irradiance ($W/m^{2}$)
and all-sky surface longwave irradiance ($W/m^{2}$). Daily data extend for 30
years (1991-01-01 to 2022-12-31). For hourly data, we reduce the
dimensionality to 7 variables, as maximum temperature, temperature range and
cloud amount are not available at this resolution. Hourly data extend for 1
year (2022-01-01 to 2022-12-31).

##### Photovoltaic dataset:

We collect photovoltaic data from Hu et al. (2022), which simulate
photovoltaic power production densely over the Continental North America.
Simulations extend over the entire year 2019, with a temporal resolution of
one hour, which we average across days to obtain daily data. The spatial
resolution of the mesh grid is 12 km, resulting in 56,776 points, from which
we randomly sample a few locations to recreate a sparse scenario. From the 13
simulated photovoltaic panels, we select module 00. Finally, we take the
average across the 21 forecast members of the provided ensemble.

##  Appendix E Additional analysis

###  E.1 Scalability and computation time

The asymptotic computational complexity of GgNet originates from the sum of
the complexities of its components. The time and space complexity are both
$\mathcal{O}(N^{2}TD)+\mathcal{O}(NTD^{2})$. As it is often safe to assume
that $D^{2}\ll N^{2}$, the overall complexity reduces to
$\mathcal{O}(N^{2}TD)$. The main bottleneck in GgNet is due to the inter-
locations graph which makes message-passing operations scale quadratically
w.r.t. the number of locations. If scalability to large $N$ is a concern,
existing sparse graph learning methods can be considered (Niculae et al.,
2023; Cini et al., 2023b). However, note that scenarios that involve
significantly more sensors than the selected datasets would not likely be
sparse and a different set of techniques should be used.

We also report here the computational times of different methods on the daily
climatic dataset ($N=235,T=10958,D=10$) under the best hyperparameter setting.
It takes approximately (on average) 2h for BRITS to complete one training, 7h
for SAITS, 16h for GgNet, 14h for GRIN. Approximate time per epoch are: 40s
for BRITS, 4m 40s for SAITS, 6m for GgNet, 8m 30s for GRIN. The temporal
modeling in GgNet is based on temporal convolutions, which makes it faster
than the popular graph-based representative GRIN, which uses a recurrent model
instead. As for BRITS and SAITS, these are faster than GgNet as they do not
account for spatial dependencies. Timings are taken on a machine equipped with
an NVIDIA A100 GPU.

###  E.2 Ablation study

The GgNet architecture is composed of the following principal constituents: 3
layers of temporal convolutions, which we indicate with $3$T, and one pairs of
alternate spatial and channel convolutions: G - $g$. The sequence is iterated
twice to allow for deeper processing, resulting in the following layers’
pattern: $2$($3$T - G - $g$). We provide in Tab. 3 an ablation study by
removing individual components; model’s hidden size is $h=64$, without
residual connections and variable embeddings (Eg).

Table 3: Ablation study. Results are averaged across 3 runs. $2$($3$T - G -
$g$) represents the architecture for which results are presented in the main
text of the paper.

Ch. | Temp. | Temp. | Temp. | Wind | Rel. | Prec. | Temp. | Clouds | Irr. | Irr. | Avg
---|---|---|---|---|---|---|---|---|---|---|---
| mean | range | max | speed | hum. |  | dew |  | short | long | (D)
2(3T - G) | ${9.2}_{\pm 0.5}$ | $25.9_{\pm 3.0}$ | $8.4_{\pm 1.2}$ | $33.1_{\pm 1.1}$ | $8.5_{\pm 0.4}$ | $83.0_{\pm 3.7}$ | $14.8_{\pm 1.6}$ | $33.1_{\pm 0.6}$ | $17.7_{\pm 1.5}$ | $4.5_{\pm 0.1}$ | $23.8_{\pm 0.7}$
2(G - $g$) | ${4.8}_{\pm 0.8}$ | $18.5_{\pm 2.4}$ | $3.9_{\pm 1.0}$ | $28.0_{\pm 2.4}$ | $5.6_{\pm 0.8}$ | $69.1_{\pm 0.7}$ | $7.6_{\pm 1.3}$ | $20.6_{\pm 1.9}$ | $11.3_{\pm 1.2}$ | $3.2_{\pm 0.2}$ | $17.3_{\pm 1.1}$
2(3T - $g$) | ${3.1}_{\pm 0.2}$ | $11.9_{\pm 0.7}$ | $2.5_{\pm 0.3}$ | $29.1_{\pm 0.9}$ | $2.8_{\pm 0.5}$ | $68.0_{\pm 3.0}$ | $4.9_{\pm 0.2}$ | $20.8_{\pm 0.4}$ | $14.1_{\pm 1.0}$ | $3.4_{\pm 0.3}$ | $16.1_{\pm 0.4}$
2(3T - G - $g$) | ${2.6}_{\pm 0.2}$ | $10.9_{\pm 1.7}$ | $2.6_{\pm 0.3}$ | $25.0_{\pm 1.4}$ | $2.8_{\pm 0.3}$ | $61.7_{\pm 1.4}$ | $4.4_{\pm 0.5}$ | $17.3_{\pm 0.4}$ | $9.2_{\pm 0.8}$ | $3.2_{\pm 0.2}$ | $14.0_{\pm 0.1}$

Considerably worse performances are obtained by removing the $g$-convolution,
i.e., $2$(3T - G), supporting the importance of modeling the relations with
the covariates in reconstructing the missing channel in a sparse setting.
Removing the temporal convolution, i.e., $2$(G - $g$), forces the model to
exchange only instantaneous information across locations and channels, missing
the temporal context. Removing the G-convolution, i.e., $2$($3$T - $g$),
negates the model from exploiting information from other locations at
inference. This variant has access to the same information as other models
that do not consider spatial information, e.g., SAITS.

###  E.3 Geographically isolated locations

In this section, we qualitatively discuss the disadvantages of graph-based
approaches in reconstructing virtual sensors in geographically sparse
locations. This is done, in Fig. 6 (left), by an illustrative example on the
daily climatic dataset.

Figure 6: (Left) Reconstruction of maximum temperature (well correlated with
other available variables) in Port-aux-Français, which cannot exploit
observations at nearby locations as geographically isolated in the south
Indian Ocean. (Right) colourmap visualization of the GgNet’s learned intra-
location graph ($g$), representation of the mutual dependencies between
climatic channels.

The plot compares the reconstruction of the maximum temperature channel, with
different models, in Port-aux-Français, a geographically isolated location in
the south Indian Ocean. For such a remote location, no information can be
inferred from nearby points, making T $\rightarrow$ Y dependencies ineffective
in performing virtual sensing. As a consequence, graph-based methods, e.g.,
GRIN, heavily underperform, as they rely on geographical proximity to build
the graph and model spatial dependencies. On the contrary, GgNet or global
recurrent models, e.g., SAITS, can effectively infer the missing channel by
leveraging, at the target location, the learned dependencies between the
target and the available covariates (C $\rightarrow$ Y dependencies). With the
intent of enhancing this effect, we deliberately choose maximum temperature as
target variable, which is highly correlated with other variables in the
dataset.

###  E.4 Learned dependencies between variables

In GgNet, relationships between channels are accounted for by means of the
intra-location graph $g$. In particular, we learn the edge scores of its
adjacency matrix from a D $\times$ D free parameter matrix, normalized with a
softmax function. In Fig. 6 (right) we provide a colourmap visualization of
the learned weighted (dense and bidirected) adjacency matrix. Interestingly,
after processing the daily climatic dataset, several meaningful correlations
have been captured: e.g., temperatures correlate with each other; temperature
range affects precipitation and wind speed; relative humidity correlates with
temperatures, precipitation and wind speed; precipitation correlate with
clouds; longwave irradiation correlates with clouds and shortwave irradiation.

###  E.5 Location-wise analysis

Tables throughout the paper present results in a channel-wise fashion, i.e.,
for each channel $d$, performance metrics are averaged across test locations.
In this additional analysis, we present, for the daily climatic dataset,
results averaged across variables, i.e., location-wise results.

Figure 7: (left) location-wise analysis of GgNet performance over SAITS for
the daily climatic dataset. (right) Robustness of GgNet to the fraction of
masked training channels on the daily climatic dataset.

With this analysis, we aim to demonstrate that the improvement of GgNet over
SAITS is not restricted to a few critical locations but is actually spread
across the considered locations. As evidenced by Fig. 7 (left), GgNet improves
over SAITS for $\sim 2/3$ of the considered target locations.

###  E.6 Robustness tests

With reference to Sec. 4.4, the training of GgNet and the adopted baselines is
guided by means of a specific mask. In particular, we mask out a fraction of
the available training channels (entirely for all timestamps in a batch) and
train the models to reconstruct such masked portions of data. This is a core
aspect, as it pushes the training toward learning the virtual-sensing task. In
this regard, we investigate the accuracy of GgNet on the daily climatic
dataset for different values of the fraction of masked training channels. The
study is reported in Fig. 7 (right); the model’s layer organization is 4T -
2(G-g), model’s hidden size is $h=64$, without residual connections and
variable embeddings (Eg). Results show a rapid improvement as soon as a small
fraction is used, followed by robust performances.

###  E.7 Uncertainty estimates

Virtual sensing models are often employed to guide decision-making at
unsampled locations. In this regard, it is important to pair the estimate for
the virtual channel with a prediction of the corresponding uncertainty.

Figure 8: Daily climatic dataset: GgNet predictions for different channels at
different locations, uncertainty estimates appear consistent with the
difficulty of the reconstruction.

In fact, channels at some locations may be, in principle, more challenging to
reconstruct than others. In GgNet, uncertainty estimates are obtained by
employing the sum of three quantile losses (Appendix C.1). Fig. 8 explores
four different scenarios where the reconstruction is gradually more
challenging. The extent of the uncertainty estimates is consistent with the
dominant factors in assessing the quality of the predictions, i.e., the
correlation between the target and the available covariates and the presence
of similar latent representations in the training set.

###  E.8 Other metrics

Results are provided in the main text in terms of the most relevant metric for
the task at hand. This is MRE for the climatic dataset, as different channels
are expressed in different units, and MAE for the photovoltaic dataset, as
power output is the only variable being reconstructed. Here, we also define a
custom metric, Variance Rescaled Error (VRE), tailored to this task:

$\text{VRE}(\hat{\textbf{x}}^{d}[n],\textbf{x}^{d}[n])=\dfrac{1}{T}\sum_{t=1}^{T}\dfrac{|\,\hat{x}^{d}_{t}[n]-x^{d}_{t}[n]\,|}{\sigma^{d}[n]}$
(18)

where $\hat{\textbf{x}}^{d}[n]$ and $\textbf{x}^{d}[n]$ are the predicted and
true $d$-th channel in $n$-th location, $\sigma^{d}[n]$ is the standard
deviation along the temporal axis of the true signal. Results are consistent
across different metrics. We support this statement by reporting below the
testing accuracy, for all experiments, in terms of the other metrics.
Specifically, Tab. 6 contains the channel-wise results for the hourly climatic
dataset that were omitted from Tab. 1 to respect the space constraints.

Table 4: Daily climatic dataset: channel-wise MAE. Results are averaged across
locations and 5 random seeds. The best-performing method is in bold, the
second-best is underlined.

Ch. | Temp. | Temp. | Temp. | Wind | Rel. | Prec. | Temp. | Clouds | Irr. | Irr.
---|---|---|---|---|---|---|---|---|---|---
| mean | range | max | speed | hum. |  | dew |  | short | long
KNN | $3.21_{\pm 0.39}$ | $2.84_{\pm 0.26}$ | $2.97_{\pm 0.35}$ | $1.47_{\pm 0.19}$ | $9.00_{\pm 0.30}$ | $2.93_{\pm 0.27}$ | $3.49_{\pm 0.39}$ | $19.53_{\pm 0.33}$ | $36.46_{\pm 1.02}$ | $21.87_{\pm 3.23}$
BRITS | $0.82_{\pm 0.09}$ | $1.19_{\pm 0.16}$ | $0.78_{\pm 0.12}$ | $1.40_{\pm 0.04}$ | $4.04_{\pm 0.64}$ | $2.30_{\pm 0.21}$ | $1.22_{\pm 0.15}$ | $17.16_{\pm 0.65}$ | $42.41_{\pm 2.18}$ | $12.49_{\pm 1.05}$
GRIN | $0.86_{\pm 0.10}$ | $1.26_{\pm 0.14}$ | $1.04_{\pm 0.11}$ | $1.33_{\pm 0.06}$ | $4.15_{\pm 0.48}$ | $2.20_{\pm 0.16}$ | $1.26_{\pm 0.12}$ | $16.70_{\pm 0.32}$ | $31.14_{\pm 1.22}$ | $14.30_{\pm 0.92}$
GRINm | $0.59_{\pm 0.16}$ | $0.88_{\pm 0.09}$ | $0.67_{\pm 0.12}$ | $1.33_{\pm 0.12}$ | $2.64_{\pm 0.65}$ | $2.00_{\pm 0.19}$ | $0.93_{\pm 0.14}$ | $\underline{11.79}_{\pm 0.39}$ | $\underline{23.70}_{\pm 0.94}$ | $13.13_{\pm 1.22}$
SAITS | $\underline{0.49}_{\pm 0.07}$ | $\underline{0.76}_{\pm 0.10}$ | $\underline{0.56}_{\pm 0.10}$ | $\underline{1.17}_{\pm 0.05}$ | $\underline{2.42}_{\pm 0.75}$ | $\underline{1.88}_{\pm 0.20}$ | $\underline{0.73}_{\pm 0.14}$ | $11.86_{\pm 0.38}$ | $28.80_{\pm 1.58}$ | $\textbf{10.07}_{\pm 0.78}$
GgNet | $\textbf{0.42}_{\pm 0.07}$ | $\textbf{0.65}_{\pm 0.02}$ | $\textbf{0.48}_{\pm 0.06}$ | $\textbf{1.04}_{\pm 0.08}$ | $\textbf{1.98}_{\pm 0.55}$ | $\textbf{1.73}_{\pm 0.13}$ | $\textbf{0.62}_{\pm 0.12}$ | $\textbf{9.48}_{\pm 0.35}$ | $\textbf{18.73}_{\pm 1.42}$ | $\underline{10.39}_{\pm 0.81}$
RNN | $1.16_{\pm 0.10}$ | $1.49_{\pm 0.13}$ | $1.40_{\pm 0.12}$ | $1.42_{\pm 0.05}$ | $5.08_{\pm 0.50}$ | $2.39_{\pm 0.20}$ | $1.57_{\pm 0.14}$ | $20.38_{\pm 0.24}$ | $44.15_{\pm 1.02}$ | $13.97_{\pm 1.17}$
RNNbi | $0.87_{\pm 0.09}$ | $1.30_{\pm 0.14}$ | $1.12_{\pm 0.10}$ | $1.35_{\pm 0.06}$ | $4.34_{\pm 0.47}$ | $2.22_{\pm 0.16}$ | $1.25_{\pm 0.13}$ | $18.06_{\pm 0.19}$ | $38.59_{\pm 0.90}$ | $13.13_{\pm 1.21}$
RNNemb | $0.85_{\pm 0.07}$ | $1.19_{\pm 0.12}$ | $1.03_{\pm 0.08}$ | $1.31_{\pm 0.07}$ | $4.01_{\pm 0.46}$ | $2.20_{\pm 0.18}$ | $1.23_{\pm 0.13}$ | $17.36_{\pm 0.26}$ | $35.48_{\pm 1.18}$ | $13.18_{\pm 0.71}$
RNNG | $0.74_{\pm 0.04}$ | $1.12_{\pm 0.12}$ | $0.94_{\pm 0.07}$ | $1.20_{\pm 0.06}$ | $3.80_{\pm 0.30}$ | $2.13_{\pm 0.20}$ | $1.06_{\pm 0.11}$ | $15.90_{\pm 0.33}$ | $28.47_{\pm 1.09}$ | $13.04_{\pm 1.13}$

Table 5: Daily climatic dataset: channel-wise VRE (%). Results are averaged
across locations and 5 random seeds. The best-performing method is in bold,
the second-best is underlined.

Ch. | Temp. | Temp. | Temp. | Wind | Rel. | Prec. | Temp. | Clouds | Irr. | Irr. | Avg
---|---|---|---|---|---|---|---|---|---|---|---
| mean | range | max | speed | hum. |  | dew |  | short | long | (D)
KNN | $105.8_{\pm 15.3}$ | $213.3_{\pm 55.6}$ | $90.9_{\pm 7.4}$ | $116.4_{\pm 26.3}$ | $102.3_{\pm 3.4}$ | $56.8_{\pm 4.2}$ | $100.8_{\pm 10.0}$ | $72.1_{\pm 2.7}$ | $59.7_{\pm 3.7}$ | $95.9_{\pm 14.8}$ | $101.4_{\pm 9.3}$
BRITS | $27.5_{\pm 3.1}$ | $75.1_{\pm 11.2}$ | $23.4_{\pm 3.4}$ | $104.0_{\pm 13.2}$ | $45.2_{\pm 6.6}$ | $41.6_{\pm 1.2}$ | $30.9_{\pm 4.1}$ | $63.7_{\pm 3.4}$ | $67.2_{\pm 4.4}$ | $53.9_{\pm 5.2}$ | $53.2_{\pm 1.9}$
GRIN | $29.7_{\pm 3.6}$ | $75.8_{\pm 12.2}$ | $31.0_{\pm 2.5}$ | $99.4_{\pm 13.7}$ | $48.2_{\pm 4.1}$ | $39.9_{\pm 1.4}$ | $36.1_{\pm 2.1}$ | $61.5_{\pm 1.6}$ | $50.4_{\pm 3.0}$ | $64.3_{\pm 9.5}$ | $53.6_{\pm 3.2}$
GRINm | $23.3_{\pm 5.1}$ | $55.6_{\pm 7.3}$ | $21.9_{\pm 5.2}$ | $101.8_{\pm 18.5}$ | $30.1_{\pm 6.6}$ | $37.0_{\pm 2.8}$ | $27.4_{\pm 3.1}$ | $\underline{43.6}_{\pm 1.8}$ | $\underline{39.6}_{\pm 3.3}$ | $64.8_{\pm 10.2}$ | $44.5_{\pm 2.7}$
SAITS | $\underline{16.8}_{\pm 1.4}$ | $\underline{44.8}_{\pm 5.5}$ | $\underline{16.5}_{\pm 2.9}$ | $\underline{88.8}_{\pm 8.9}$ | $\underline{26.0}_{\pm 7.4}$ | $\underline{33.7}_{\pm 1.6}$ | $\underline{20.1}_{\pm 3.2}$ | $44.0_{\pm 2.1}$ | $47.5_{\pm 5.4}$ | $\textbf{45.2}_{\pm 5.7}$ | $38.3_{\pm 1.5}$
GgNet | $\textbf{14.6}_{\pm 3.2}$ | $\textbf{35.8}_{\pm 3.0}$ | $\textbf{14.4}_{\pm 1.8}$ | $\textbf{80.5}_{\pm 14.5}$ | $\textbf{22.1}_{\pm 5.3}$ | $\textbf{30.9}_{\pm 1.3}$ | $\textbf{17.3}_{\pm 2.2}$ | $\textbf{35.2}_{\pm 1.7}$ | $\textbf{31.3}_{\pm 3.8}$ | $\underline{48.9}_{\pm 6.3}$ | $\textbf{33.1}_{\pm 1.9}$
RNN | $37.2_{\pm 4.1}$ | $89.6_{\pm 5.4}$ | $39.3_{\pm 3.6}$ | $105.2_{\pm 12.0}$ | $56.8_{\pm 4.3}$ | $44.0_{\pm 2.2}$ | $42.0_{\pm 2.6}$ | $74.7_{\pm 1.1}$ | $69.1_{\pm 3.3}$ | $60.0_{\pm 7.0}$ | $61.8_{\pm 1.1}$
RNNbi | $28.9_{\pm 3.4}$ | $74.8_{\pm 9.8}$ | $32.7_{\pm 3.3}$ | $101.0_{\pm 11.6}$ | $48.7_{\pm 4.2}$ | $40.9_{\pm 2.2}$ | $33.3_{\pm 2.0}$ | $66.4_{\pm 1.8}$ | $61.4_{\pm 3.9}$ | $57.3_{\pm 7.5}$ | $54.5_{\pm 0.8}$
RNNemb | $26.9_{\pm 4.1}$ | $63.2_{\pm 3.8}$ | $29.7_{\pm 3.0}$ | $99.8_{\pm 12.7}$ | $45.8_{\pm 3.9}$ | $40.0_{\pm 1.4}$ | $33.3_{\pm 3.0}$ | $63.8_{\pm 1.4}$ | $57.3_{\pm 4.6}$ | $58.4_{\pm 7.8}$ | $51.8_{\pm 2.3}$
RNNG | $24.5_{\pm 3.2}$ | $60.2_{\pm 2.5}$ | $28.0_{\pm 2.3}$ | $93.1_{\pm 11.8}$ | $43.6_{\pm 2.4}$ | $38.6_{\pm 2.0}$ | $29.0_{\pm 2.1}$ | $58.5_{\pm 1.5}$ | $46.2_{\pm 3.3}$ | $58.6_{\pm 10.5}$ | $48.0_{\pm 2.0}$

Table 6: Hourly climatic dataset: channel-wise MRE (%). Results are averaged
across locations and 5 random seeds. The best-performing method is in bold,
the second-best is underlined.

Ch. | Temp. | Wind | Rel. | Prec. | Temp. | Irr. | Irr. | Avg
---|---|---|---|---|---|---|---|---
| mean | speed | hum. |  | dew | short | long | (H)
KNN | $15.9_{\pm 2.4}$ | $38.0_{\pm 5.2}$ | $14.6_{\pm 1.8}$ | $133.6_{\pm 7.8}$ | $23.1_{\pm 3.4}$ | $21.4_{\pm 1.9}$ | $6.7_{\pm 0.4}$ | $36.2_{\pm 1.9}$
BRITS | $7.4_{\pm 2.5}$ | $34.6_{\pm 1.4}$ | $6.5_{\pm 1.2}$ | $85.4_{\pm 1.1}$ | $9.1_{\pm 1.5}$ | $36.5_{\pm 2.0}$ | $4.6_{\pm 0.4}$ | $26.3_{\pm 0.8}$
GRIN | $6.5_{\pm 2.3}$ | $35.5_{\pm 1.0}$ | $5.3_{\pm 1.1}$ | $86.9_{\pm 2.3}$ | $8.6_{\pm 1.2}$ | $\textbf{16.2}_{\pm 1.0}$ | $4.4_{\pm 0.3}$ | $23.3_{\pm 0.7}$
GRINm | $6.0_{\pm 2.2}$ | $34.5_{\pm 2.1}$ | $4.8_{\pm 0.9}$ | $85.5_{\pm 1.1}$ | $8.0_{\pm 1.4}$ | $\textbf{16.1}_{\pm 0.8}$ | $4.3_{\pm 0.3}$ | $22.7_{\pm 0.4}$
SAITS | $\underline{4.7}_{\pm 1.8}$ | $\textbf{29.0}_{\pm 0.8}$ | $\textbf{4.3}_{\pm 0.8}$ | $\underline{82.8}_{\pm 1.2}$ | $\underline{6.6}_{\pm 1.2}$ | $24.4_{\pm 1.6}$ | $\underline{3.7}_{\pm 0.2}$ | $\underline{22.2}_{\pm 0.5}$
GgNet | $\textbf{4.5}_{\pm 1.9}$ | $\underline{29.2}_{\pm 2.5}$ | $\textbf{4.3}_{\pm 0.7}$ | $\textbf{78.4}_{\pm 1.1}$ | $\textbf{5.7}_{\pm 0.9}$ | $\underline{17.4}_{\pm 3.8}$ | $\textbf{3.6}_{\pm 0.2}$ | $\textbf{20.4}_{\pm 0.6}$
RNN | $8.0_{\pm 2.7}$ | $38.9_{\pm 1.3}$ | $6.6_{\pm 0.8}$ | $90.8_{\pm 1.3}$ | $9.7_{\pm 1.5}$ | $43.2_{\pm 3.3}$ | $4.5_{\pm 0.3}$ | $28.8_{\pm 0.8}$
RNNbi | $5.9_{\pm 2.0}$ | $36.5_{\pm 1.5}$ | $5.3_{\pm 0.9}$ | $87.1_{\pm 1.3}$ | $7.7_{\pm 1.2}$ | $33.3_{\pm 2.0}$ | $4.1_{\pm 0.3}$ | $25.7_{\pm 0.5}$
RNNemb | $6.6_{\pm 2.8}$ | $36.4_{\pm 0.9}$ | $5.6_{\pm 1.1}$ | $86.8_{\pm 1.1}$ | $8.2_{\pm 1.5}$ | $35.6_{\pm 3.2}$ | $4.5_{\pm 0.4}$ | $26.2_{\pm 0.9}$
RNNG | $5.7_{\pm 2.4}$ | $32.2_{\pm 2.7}$ | $5.1_{\pm 0.9}$ | $83.0_{\pm 3.0}$ | $7.1_{\pm 1.0}$ | $25.1_{\pm 2.7}$ | $4.0_{\pm 0.4}$ | $23.2_{\pm 1.1}$

Table 7: Hourly climatic dataset: channel-wise MAE. Results are averaged
across locations and 5 random seeds. The best-performing method is in bold,
the second-best is underlined.

Ch. | Temp. | Wind | Rel. | Prec. | Temp. | Irr. | Irr.
---|---|---|---|---|---|---|---
| mean | speed | hum. |  | dew | short | long
KNN | $3.24_{\pm 0.29}$ | $1.64_{\pm 0.25}$ | $10.59_{\pm 1.05}$ | $0.16_{\pm 0.01}$ | $3.54_{\pm 0.49}$ | $43.52_{\pm 4.19}$ | $24.54_{\pm 1.55}$
BRITS | $1.49_{\pm 0.44}$ | $1.49_{\pm 0.13}$ | $4.68_{\pm 0.76}$ | $\underline{0.10}_{\pm 0.01}$ | $1.39_{\pm 0.24}$ | $74.13_{\pm 4.00}$ | $16.84_{\pm 1.39}$
GRIN | $1.32_{\pm 0.42}$ | $1.53_{\pm 0.09}$ | $3.84_{\pm 0.74}$ | $\underline{0.10}_{\pm 0.01}$ | $1.31_{\pm 0.20}$ | $\textbf{32.91}_{\pm 2.22}$ | $15.99_{\pm 1.11}$
GRINm | $1.23_{\pm 0.40}$ | $1.49_{\pm 0.16}$ | $3.47_{\pm 0.60}$ | $\underline{0.10}_{\pm 0.00}$ | $1.22_{\pm 0.23}$ | $\textbf{32.74}_{\pm 1.75}$ | $15.60_{\pm 0.99}$
SAITS | $\underline{0.96}_{\pm 0.32}$ | $\textbf{1.25}_{\pm 0.08}$ | $\textbf{3.11}_{\pm 0.50}$ | $\underline{0.10}_{\pm 0.01}$ | $\underline{1.00}_{\pm 0.18}$ | $49.67_{\pm 3.34}$ | $\underline{13.31}_{\pm 0.77}$
GgNet | $\textbf{0.91}_{\pm 0.35}$ | $\underline{1.26}_{\pm 0.14}$ | $\textbf{3.11}_{\pm 0.47}$ | $\textbf{0.09}_{\pm 0.01}$ | $\textbf{0.87}_{\pm 0.15}$ | $\underline{35.28}_{\pm 7.33}$ | $\textbf{13.14}_{\pm 0.49}$
RNN | $1.61_{\pm 0.48}$ | $1.68_{\pm 0.12}$ | $4.76_{\pm 0.50}$ | $0.11_{\pm 0.00}$ | $1.49_{\pm 0.23}$ | $87.86_{\pm 6.70}$ | $16.21_{\pm 1.09}$
RNNbi | $1.19_{\pm 0.37}$ | $1.57_{\pm 0.13}$ | $3.82_{\pm 0.59}$ | $0.10_{\pm 0.01}$ | $1.18_{\pm 0.18}$ | $67.75_{\pm 4.15}$ | $15.00_{\pm 1.08}$
RNNemb | $1.34_{\pm 0.51}$ | $1.57_{\pm 0.09}$ | $4.09_{\pm 0.71}$ | $0.10_{\pm 0.01}$ | $1.25_{\pm 0.24}$ | $72.50_{\pm 7.11}$ | $16.29_{\pm 1.30}$
RNNG | $1.14_{\pm 0.43}$ | $1.39_{\pm 0.14}$ | $3.72_{\pm 0.62}$ | $0.10_{\pm 0.00}$ | $1.09_{\pm 0.17}$ | $50.99_{\pm 5.59}$ | $14.71_{\pm 1.28}$

Table 8: Hourly climatic dataset: channel-wise VRE (%). Results are averaged
across locations and 5 random seeds. The best-performing method is in bold,
the second-best is underlined.

Ch. | Temp. | Wind | Rel. | Prec. | Temp. | Irr. | Irr. | Avg
---|---|---|---|---|---|---|---|---
| mean | speed | hum. |  | dew | short | long | (H)
KNN | $77.4_{\pm 5.2}$ | $96.9_{\pm 13.9}$ | $83.9_{\pm 11.9}$ | $49.1_{\pm 6.6}$ | $106.9_{\pm 21.2}$ | $16.3_{\pm 1.4}$ | $83.5_{\pm 8.3}$ | $73.4_{\pm 4.1}$
BRITS | $35.9_{\pm 7.0}$ | $83.3_{\pm 4.4}$ | $36.1_{\pm 5.8}$ | $29.0_{\pm 0.8}$ | $38.4_{\pm 5.0}$ | $26.8_{\pm 1.5}$ | $59.2_{\pm 6.5}$ | $44.1_{\pm 2.9}$
GRIN | $37.8_{\pm 7.1}$ | $85.0_{\pm 2.2}$ | $32.6_{\pm 8.2}$ | $29.2_{\pm 1.1}$ | $42.7_{\pm 8.6}$ | $\textbf{12.2}_{\pm 0.8}$ | $59.1_{\pm 4.3}$ | $42.7_{\pm 1.1}$
GRINm | $40.1_{\pm 8.3}$ | $82.8_{\pm 6.3}$ | $30.0_{\pm 5.5}$ | $29.2_{\pm 0.4}$ | $40.0_{\pm 8.3}$ | $\textbf{12.1}_{\pm 0.6}$ | $58.2_{\pm 7.4}$ | $41.8_{\pm 2.7}$
SAITS | $\underline{24.3}_{\pm 4.3}$ | $\textbf{69.9}_{\pm 4.0}$ | $\textbf{25.5}_{\pm 3.5}$ | $\underline{28.3}_{\pm 0.6}$ | $\underline{28.1}_{\pm 3.7}$ | $18.3_{\pm 1.3}$ | $\textbf{46.5}_{\pm 1.1}$ | $\underline{34.4}_{\pm 1.3}$
GgNet | $\textbf{24.0}_{\pm 6.0}$ | $\underline{70.6}_{\pm 2.8}$ | $\underline{26.3}_{\pm 5.2}$ | $\textbf{27.3}_{\pm 1.4}$ | $\textbf{25.3}_{\pm 3.1}$ | $\underline{12.9}_{\pm 2.5}$ | $\underline{46.9}_{\pm 4.4}$ | $\textbf{33.3}_{\pm 0.9}$
RNN | $40.5_{\pm 7.2}$ | $93.5_{\pm 4.7}$ | $36.5_{\pm 4.2}$ | $31.2_{\pm 0.9}$ | $44.3_{\pm 4.5}$ | $31.8_{\pm 2.3}$ | $56.0_{\pm 2.0}$ | $47.7_{\pm 2.0}$
RNNbi | $29.9_{\pm 5.9}$ | $87.4_{\pm 4.3}$ | $30.5_{\pm 5.4}$ | $29.5_{\pm 1.0}$ | $33.3_{\pm 3.1}$ | $24.8_{\pm 1.4}$ | $51.7_{\pm 2.2}$ | $41.0_{\pm 1.3}$
RNNemb | $35.6_{\pm 8.4}$ | $88.4_{\pm 3.2}$ | $33.8_{\pm 6.9}$ | $29.1_{\pm 1.0}$ | $33.9_{\pm 2.9}$ | $26.2_{\pm 2.5}$ | $57.2_{\pm 1.8}$ | $43.4_{\pm 2.2}$
RNNG | $32.9_{\pm 11.1}$ | $78.0_{\pm 4.6}$ | $31.0_{\pm 6.6}$ | $27.8_{\pm 1.6}$ | $33.7_{\pm 5.0}$ | $18.5_{\pm 1.8}$ | $52.9_{\pm 2.2}$ | $39.3_{\pm 2.7}$

Table 9: Photovoltaic dataset: MRE (%). For each N, results are averaged
across 5 sampling of the locations and 5 different runs for each set (25 runs
in total). The best-performing method is in bold, the second-best is
underlined.

N${}_{\text{locations}}$ | 10 | 20 | 30 | 50 | 70 | 100 | 150 | 200
---|---|---|---|---|---|---|---|---
GRIN | $20.1\pm 6.3$ | $18.8\pm 2.7$ | $18.6\pm 2.2$ | $18.2\pm 2.2$ | $18.8\pm 1.8$ | $18.4\pm 1.4$ | $18.3\pm 1.7$ | $18.5\pm 1.6$
GRINm | $17.2\pm 4.7$ | $16.0\pm 2.4$ | $\underline{15.6}\pm 2.1$ | $\underline{15.0}\pm 1.5$ | $\underline{15.1}\pm 1.2$ | $14.8\pm 0.8$ | $\underline{14.6}\pm 0.7$ | $\underline{14.6}\pm 0.6$
SAITS | $\underline{16.0}\pm 3.0$ | $\textbf{15.8}\pm 2.3$ | $\underline{15.6}\pm 1.8$ | $15.1\pm 1.4$ | $15.4\pm 1.3$ | $\underline{14.7}\pm 1.0$ | $\underline{14.6}\pm 0.7$ | $\underline{14.6}\pm 0.7$
GgNet | $\textbf{15.6}\pm 3.5$ | $\underline{15.9}\pm 2.0$ | $\textbf{15.3}\pm 1.9$ | $\textbf{14.7}\pm 1.5$ | $\textbf{14.5}\pm 1.3$ | $\textbf{13.7}\pm 1.0$ | $\textbf{12.8}\pm 0.7$ | $\textbf{12.5}\pm 0.7$

Table 10: Photovoltaic dataset: VRE (%). For each N, results are averaged
across 5 sampling of the locations and 5 different runs for each set (25 runs
in total). The best-performing method is in bold, the second-best is
underlined.

N${}_{\text{locations}}$ | 10 | 20 | 30 | 50 | 70 | 100 | 150 | 200
---|---|---|---|---|---|---|---|---
GRIN | $40.7\pm 11.9$ | $38.2\pm 4.3$ | $37.8\pm 3.6$ | $36.0\pm 2.5$ | $37.4\pm 2.7$ | $37.0\pm 2.1$ | $36.9\pm 2.8$ | $37.4\pm 3.0$
GRINm | $34.8\pm 7.1$ | $32.8\pm 4.1$ | $\underline{31.7}\pm 3.3$ | $\underline{29.8}\pm 1.8$ | $\underline{30.3}\pm 1.7$ | $30.0\pm 1.2$ | $\underline{29.5}\pm 1.2$ | $29.8\pm 1.2$
SAITS | $\underline{32.8}\pm 4.6$ | $\textbf{32.4}\pm 3.8$ | $\underline{31.7}\pm 2.6$ | $30.0\pm 2.0$ | $30.8\pm 1.9$ | $\underline{29.7}\pm 1.3$ | $\underline{29.5}\pm 1.2$ | $\underline{29.6}\pm 1.3$
GgNet | $\textbf{31.8}\pm 4.4$ | $\underline{32.6}\pm 3.4$ | $\textbf{31.2}\pm 2.6$ | $\textbf{29.2}\pm 2.0$ | $\textbf{29.0}\pm 2.0$ | $\textbf{27.8}\pm 1.7$ | $\textbf{26.0}\pm 1.4$ | $\textbf{25.5}\pm 1.4$