Kernel generalized least squares regression for network-structured data

TL;DR

This paper introduces new methods for predicting network-structured data using regression models, showing improved performance in time-series applications like air quality monitoring.

Contribution

The paper extends Kernel Graph Regression (KGR) to handle missing data and correlated errors, and introduces a scalable GLS-KGR algorithm.

Findings

The GLS-KGR algorithm outperforms standard techniques in time-series applications with network-structured data.

A scalable expression for marginal variance is derived using the Laplace approximation for out-of-sample prediction error.

The proposed methods are validated on both synthetic and real-world air quality data from California.

Abstract

In this paper, we study a class of non-parametric regression models for predicting graph signals {𝐲t} as a function of explanatory variables {𝐱t}. Recently, Kernel Graph Regression (KGR) and Gaussian Processes over Graph (GPoG) have emerged as promising techniques for this task. The goal of this paper is to examine several extensions to KGR/GPoG, with the aim of generalising them a wider variety of data scenarios. The first extension we consider is the case of graph signals that have only been partially recorded, meaning a subset of their elements is missing at observation time. Next, we examine the statistical effect of correlated prediction error and propose a method for Generalized Least Squares (GLS) on graphs. In particular, we examine Autoregressive AR(1) vector autoregressive processes, which are commonly found in time-series applications. Finally, we use the Laplace approximation to determine a lower bound for the out-of-sample prediction error and derive a scalable expression for the marginal variance of each prediction. These methods are tested on both real and synthetic data, with the former taken from a network of air quality monitoring stations across California. We find evidence that the generalised GLS-KGR algorithm is well-suited to such time-series applications, outperforming several standard techniques on this dataset.

Figures

Colour representation of seven low-frequency eigenvectors and two high-frequency eigenvectors of a 3D mesh graph, along with their associated eigenvalue.The mesh is constructed from the classic Utah Teapot 3D object. Here, the main body, lid and spout are all mutually disconnected, giving three zero eigenvalues.

A graphical depiction of a small Cartesian product graph.

The mean number of iterations required for convergence on synthetic data is shown over a range of values for θ and α.The upper plot shows the results of experiment (a), when ΣN=𝐈, and the lower plot shows the results for experiment (b), when ΣNij=exp((i−j)2)+0.05δij.

The Root Mean Squared Error on the test and train sets are shown as γ is varied.

The Root Mean Squared Error on the test and train sets are shown as β is varied.

The Root Mean Squared Error on the test and train sets are shown as σ is varied.

The Root Mean Squared Error on the test and train sets are shown as a function of the total percentage of node-times that were observed.

A graphical depiction of the graph connecting the 1382 available monitoring stations, with the map colour representing the elevation profile.This is the graph used in all further experiments. “Reprinted from [38] under a CC BY license, with permission from [Heriott-Watt University], original copyright [2024].” Source available at: http://hdl.handle.net/10399/5041

Left: the observed graph signal for Ozone on a particular day represented via a colour map.Right: prediction for the reconstructed latent graph signal across the entire network on the same day made by the GLS KGR method. “Reprinted from [38] under a CC BY license, with permission from [Heriott-Watt University], original copyright [2024].” Source available at: http://hdl.handle.net/10399/5041

A section of the Ozone time-series signal predicted by GLS KGR at a particular unobserved node is shown along with the ground truth.The blue shading depicts two standard deviations of prediction error arising from latent signal uncertainty as calculated by the Laplace approximation.

Equations

Full text

1 Introduction

Graphs have proven to be a useful way of describing complex datasets due to their ability to capture general relational information between entities [1]. Knowledge about pairwise relationships between data points can, for example, help enhance regularization, improve computational efficiency and robustness, or better exploit unlabelled data [2,3]. Successful applications have included social networks [4], brain imaging [5], biomolecular systems [6], sensor networks [7], image and point cloud processing [8] and chemical structure [9] to name a few.

The Graph Signal Processing (GSP) community focuses in particular on generalisations of traditional signal processing techniques on regular domains. A toolbox of analogous graph-based techniques has been built by translating concepts such as convolution and spectral decomposition into the graph domain. These have been applied to a diverse range of problems including signal reconstruction [10], graph learning [11], denoising [12] and regression [13].

Graph-based regression specifically has received some attention recently in the form of kernel regression and Gaussian processes [13–16]. In this context, we are interested in constructing a function that maps input features to an output space, which is interpreted as a graph signal. This opens interesting avenues for enhanced regularisation using GSP tools and offers novel ways of defining kernel functions.

In this paper, we examine Kernel Graph Regression (KGR) and the closely related topic of Gaussian Processes over Graph (GPoG) from the perspective of Graph Signal Reconstruction (GSR). This framework allows signal prediction at unobserved test nodes as well as test times. Next, we extend this by developing a method for Generalized Least Squares (GLS) regression on graphs. This has particular relevance for the statistical analysis of time-series data where a graph signal is sampled at regular intervals, as it allows, for example, node-level error autocorrelation to be incorporated. Such scenarios are encountered in the analysis of many common graph applications such as sensor networks and financial time series. Finally, we use the Laplace approximation to calculate a lower bound for the out-of-sample prediction uncertainty and derive a scalable expression for the marginal uncertainty at each node.

1.1 Problem overview and scope

The goal of supervised learning and regression is to estimate an optimal function that maps an input space to an output space, given a set of training examples $[eqn]$ . It will be assumed that there is a static graph $[eqn]$ with vertex set $[eqn]$ with N vertices (nodes), edge set $[eqn]$ and adjacency matrix $[eqn]$ . In this paper we are concerned specifically with real-valued vector regression inputs (time series covariates) and outputs (time series responses), that is, $[eqn]$ and $[eqn]$ , for d-dim observation vector over a subset $[eqn]$ of $[eqn]$ vertices out of a total of N vertices at time t. Clearly, it is assumed that $[eqn]$ is a partially observed graph signal, meaning the elements can be interpreted as being observed at time t on a subset $[eqn]$ of the full node-set $[eqn]$ .

In this setup, the M dimensional covariate time series $[eqn]$ may be exogenous to the graph vertices (not observed at any of the $[eqn]$ graph vertices) or they may be endogenous to the graph vertices (observed at a subset or all of the graph vertices). The methodology developed in this work admits both possible setups, but the focus in this work will be primarily on the former case of exogeneity of the regression covariate time series and this case will be considered unless otherwise stated.

The $[eqn]$ dimensional response time series vectors $[eqn]$ of the graph regression will be assumed to be observed, ideally at each time point $[eqn]$ , such that the first d coordinates $[eqn]$ for $[eqn]$ corresponds to the d-dimensional observation sub-vector of $[eqn]$ that were observed at vertex $[eqn]$ and analogously coordinates $[eqn]$ for $[eqn]$ corresponds to the d-dimensional observation sub-vector of $[eqn]$ for vertex $[eqn]$ .

In the case that the cardinality of vertex set $[eqn]$ is less than that of the complete graph vertex set $[eqn]$ i.e. $[eqn]$ , the formulation will accommodate missingness of observations at some vertices in the setup, allowing for a subset of the vertices, $[eqn]$ , to have complete missingness of the observation vectors per time instant. In the case that $[eqn]$ one may have full or partial observations all d-dim observation sub-vector coordinates per vertex $[eqn]$ at all vertices in $[eqn]$ .

This gives a general pattern of missingness where single observation coordinates of missingness may occur at a given vertex or an entire observation vector may be missing at a vertex at any instant of time. To capture this general pattern of missingness, a general sensing matrix $[eqn]$ will be utilised throughout, where Si,t takes value 1 if the i-th coordinate of $[eqn]$ is observed at time t and value 0 if it is missing at time t. Without loss of generality, in latter sections, we will set d = 1 and we will decompose S into two sub-matrices $[eqn]$ and $[eqn]$ .

The graph itself is assumed to be static in time with a known adjacency matrix $[eqn]$ , for a graph with N vertices. There will be an estimation challenge to estimate the graph signal reconstruction and regression models over T time points based on a potentially partially observed set of responses over N vertices at $[eqn]$ time points. The types of missingness that can be considered in this setup is quiet general and this aspect will be explained with illustrations in latter sections of the paper.

In this formulation, the input features $[eqn]$ are ‘global’ in the sense that individual elements are not necessarily associated with any particular node in the output space. For example, if the graph represents a network of corporations, and the graph signal is their quarterly revenue growth, the features could represent global factors such as inflation, GDP etc, as well as firm-specific data such as the number of employees at a particular company. This is the most general formulation and allows for maximum flexibility in problem specification. Additional constraints explicitly linking elements from one space to the other are possible but not considered in this paper. This approach is shared by other recent work such as [13,17].

Throughout this paper we take a Bayesian perspective, using graph filters to construct spectral priors that can be used to make statements about the expected profile of a graph signal. Using graph filters in this way shares aesthetic and practical features with the graph kernel framework described in [18,19] and elsewhere. Graph filters as opposed to kernels, however, provide a more direct Bayesian interpretation which is useful in this context.

As mentioned, in this work, the graph adjacency matrix $[eqn]$ is assumed to be known a priori. In many applications such as social or sensor networks, this is a reasonable assumption, as a graph may be self-evident or easy to construct. For other applications, such as financial or biological networks, the dependency structure may be more opaque. In such scenarios, we refer the reader to the large body of work on graph learning [20–23].

1.2 Contributions

This paper contributes a statistical exploration of Kernel Graph Regression and Gaussian Processes over Graphs in the context of graph signal reconstruction. A primary aim is to incorporate the possibility of error correlation, which is of practical concern in many applications. We pay particular attention to time series modelling and node-level autocorrelation and describe an algorithm for determining function and model parameters in the presence of autocorrelated noise. We also focus on deriving mathematical expressions which lend themselves to practical computation for large-scale problems and use Bayesian reasoning to provide a lower bound for the model parameter uncertainty. These contributions are clarified below.

1.2.1 Incorporating partially observed graph signals

Prior work on KGR and GPoG has assumed that the signal from all nodes under consideration is fully observed throughout training [24], [25] and [26]. We introduce a modification that incorporates ideas from graph signal reconstruction to allow a constant subset of nodes to remain unobserved at train time and use the topology of the graph to make smooth predictions at these points.

1.2.2 Incorporating correlation prediction error via Generalised Least Squares (GLS)

We relax the statistical assumption that all prediction errors should be independent and identically distributed, allowing for the possibility of cross-correlation amongst nodes, and autocorrelation over time. This creates a broader class of models with rich statistical properties. In particular, we consider the scenario where Autoregressive AR(1) autocorrelation may exist at each node which has relevance in many time-series applications.

1.2.3 Strict enforcement of worst-case O(N3  +  T3) complexity

Graph regression problems are naturally expressed in terms of matrix operations of $[eqn]$ $[eqn]$ $[eqn]$ Kronecker products i.e in terms of a matrix operations on an NT $[eqn]$ NT matrix of elements involving a high dimensional inverse. This is problematic for all but the smallest of problems as memory and computational costs escalate quickly. This work contributes an analysis of the complexity of the methods and ensures that all final mathematical expressions are given in terms of worst-case $[eqn]$ operations, ensuring medium to large problems remain tractable.

1.2.4 A Laplace approximation for parameter uncertainty

We use the Laplace approximation to estimate the full posterior and derive a tractable expression for the marginal variance. This can be used to provide a lower bound for the uncertainty over out-of-sample model predictions.

1.3 Paper organisation

Section 2 gives a brief introduction to some fundamental GSP concepts such as the Graph Fourier Transform (GFT) and graph filters. In section 3 we revisit the problem of Kernel Graph Regression (KGR) from a Bayesian perspective and incorporate the concept of partially observed graph signals. Section 4 builds on this, introducing the GLS model and outlining an iterative algorithm for estimating the prediction error, focusing particularly on an autoregressive model. Here, we also run several experiements on synthetic data to provide intuition for the hyperparameters this algorithm relies on. Finally, section 5 uses the Laplace approximation to derive a tractable lower bound for the marginal prediction uncertainty. Finally, section 6 analyses these methods on a real dataset concerning the prediction of pollutant levels across a network of monitoring stations in California.

1.4 Notation

Effort is made to adhere to the following variable naming conventions throughout this paper, see Table 1.

Table 1: The notational conventions used in this paper.

2. Preliminaries

2.1 Graphs and the GFT

Consider a weighted, undirected graph $[eqn]$ made up of N vertices $[eqn]$ and edge set $[eqn]$ connecting vertex pairs. $[eqn]$ is the weighted adjacency matrix, where $[eqn]$ represents the strength of interaction between nodes i and j. $[eqn]$ implies $[eqn]$ and $[eqn]$ implies $[eqn]$ . For graph $[eqn]$ one may define the graph Laplacian matrix as $[eqn]$ − $[eqn]$ where $[eqn]$ is the diagonal degree matrix; $[eqn]$ being the ith row sum (or column sum) of $[eqn]$ . One key property of the graph Laplacian is that it can be used to measure the smoothness of a signal $[eqn]$ with respect to the graph.

$$[eqn]$$

As visible, the quadratic product $[eqn]$ measures the square difference in signal value between connected nodes, weighted by the respective adjacency matrix entry. In this way, it can be considered a measure of ‘roughness’ of the signal f, with smaller values indicating a smoother signal, achieving a minimum of zero when the function is constant on each disconnected subgraph [27]. It can be useful to view the Laplacian in terms of its eigendecomposition, $[eqn]$ , where the columns of the orthogonal matrix $[eqn]$ are the normalised eigenvectors of $[eqn]$ , forming an orthonormal spanning set. (The eigenvalues are typically ordered such that $[eqn]$ ). Any graph signal can be expressed as a linear combination of the eigenvectors of the Laplacian as $[eqn]$ , where $[eqn]$ is the i-th column of the matrix $[eqn]$ and $[eqn]$ is known as the Graph Fourier Transform of f [28]. This gives an alternative way to express the quadratic form of Eq (1).

$$[eqn]$$

These two facts provide the interpretation of the eigenvectors of the Laplacian as having varying degrees of smoothness, ordered by the size of the associated eigenvalue. The number of zero eigenvalues is the number of disconnected subgraphs. Beyond that, larger eigenvalues, analogous to frequency in classical signal processing, are associated with increasingly ‘rough’ eigenvectors. Fig 1 shows an example of this behaviour for a 3D mesh graph.

Colour representation of seven low-frequency eigenvectors and two high-frequency eigenvectors of a 3D mesh graph, along with their associated eigenvalue.The mesh is constructed from the classic Utah Teapot 3D object. Here, the main body, lid and spout are all mutually disconnected, giving three zero eigenvalues.

2.2 Graph filters

The concept of a graph filter can be constructed by extending the analogy of the GFT [29]. A low-pass filter, $[eqn]$ , which operates on a graph signal f, attenuates the high-frequency graph Fourier components and can be used to smooth or denoise a signal [30]. It can be constructed by defining a decreasing function $[eqn]$ , parametrised by $[eqn]$ , though we will suppress the parameter when using the function unless explicitly required, for which $[eqn]$ and $[eqn]$ . The filter, $[eqn]$ , is then defined by

$$[eqn]$$

where $[eqn]$ denotes the application of $[eqn]$ element wise to the diagonal. The filtered signal $[eqn]$ is then calculated as $[eqn]$ . Note that the operator $[eqn]$ is necessarily symmetric and semi-positive definite. Numerous possibilities exist for $[eqn]$ , some of which are outlined in Table 2.

Table 2: Possible graph filter functions.

Graph filters constructed in this way bear a close resemblance to graph kernels as described in [18], where instead the authors define an increasing function $[eqn]$ . In practical terms, a graph filter can be considered the inverse of a graph kernel as defined in this paper or, more accurately, $[eqn]$ . “Graph kernel" in this context should not be confused with another concept of the same name, which concerns distance metrics between pairs of graphs [31].

2.3 Graph spectral priors

Graph filters can also be used to construct covariance matrices, for example in [14,32]. Consider a Gaussian white noise graph signal $[eqn]$ with precision $[eqn]$ . If a filter $[eqn]$ is applied to give $[eqn]$ , the resultant signal will be smoother with respect to the graph and will be drawn from a distribution

$$[eqn]$$

In a Bayesian setting, $[eqn]$ can act as a covariance matrix which encodes the assumption that signals observed over a graph are likely to be smooth with respect to the underlying topology.

The selection of an appropriate graph filter in this context can be guided by prior knowledge about the expected power spectrum. Adaptive techniques have also been proposed for jointly learning a graph filter such as [15,33] but are not considered in this paper.

3. Kernel graph regression

3.1 Data and Definitions

This problem runs over a set of T discrete “times" $[eqn]$ . (These need not be regularly-sampled time intervals but in many practical problems will be). At each time, there exists a fixed graph $[eqn]$ with N nodes as specified earlier by vertex set $[eqn]$ and related by a static adjacency matrix $[eqn]$ .

At a distinct subset of times $[eqn]$ where $[eqn]$ we observe supervised learning pairs $[eqn]$ . Here, $[eqn]$ is a vector of explanatory variables, and $[eqn]$ is a partially observed graph signal. This signal is measured on a subset of the nodes $[eqn]$ where $[eqn]$ which remains constant over all $[eqn]$ . In a sensor network, this may be the case if only a subset of the sensors are functioning correctly, or we wish to estimate the reading at locations where no physical sensor is placed. The goal is to predict the graph signal at the unobserved nodes for $[eqn]$ and at all nodes for $[eqn]$ .

It is useful to define two binary selection matrices $[eqn]$ and $[eqn]$ which help in mapping between observed and unobserved node sets. $[eqn]$ has $[eqn]$ -th entry equal to one for $[eqn]$ to $[eqn]$ and $[eqn]$ , with the rest set to zero. Similarly $[eqn]$ has entries $[eqn]$ equal to one for $[eqn]$ to $[eqn]$ and $[eqn]$ , with the rest set to zero. For simplicity, and without loss of generality, we can assume that

$$[eqn]$$

which implies that the first $[eqn]$ nodes are observed, and the first $[eqn]$ times contain the labelled examples. Algebraically, this is always possible under a reordering of nodes/times. However, statistically, depending on the regression assumptions of the model, one should note that such reordering may influence the regression results, depending on what assumptions on the autocorrelation structures and conditional dependence assumptions of the response variables, given the covariates, are made for the observations over time in the graph regression. If the observations are assumed conditionally uncorrelated (or independent if under a Gaussian error assumption) in time, given the covariates, then reordering in this manner will not influence the regression results, in the case that the covariates time series $[eqn]$ does not display autocorrelation or cross autocorrelation in any coordinates (of course assuming covariates in the regression are also reordered accordingly). If however, the covariates do have a temporal cross correlation then mixing assumptions may need to be considered when deciding about the suitability of reordering, such that reordering takes place only between time intervals with lags that are sufficiently uncorrelated in the design space covariate time series. If additionally, there is deemed to be a conditional autocorrelation in the response time series, even after conditioning on regressors, in the graph regression, then such reordering should be considered carefully as it may influence the results of the regression including the smoothness of signals learnt. In this instance, one may perform a decorrelating transformation on the observation covariates prior to performing the reordering, of course, this would assume knowledge of the conditional autocorrelation matrix of the response in the regression. When this is not known, it leads to a generalised iterative Graph regression that is the extension of Generalised Least Squares (see [34]) into the graph regression setting. Such an example is illustrated in the AR(1) error process in Section 4.3 below. Alternatively, if the assumptions required to perform re-ordering are not easily satisfied for a given regression application, then there is no problem with working with matrices $[eqn]$ and $[eqn]$ without any reordering, it is just less convenient algebraically.

3.2 KGR via signal reconstruction

We will now derive a Kernel Graph Regression model from the perspective of graph signal smoothing and reconstruction. Let us assume that there exists an underlying noiseless graph signal $[eqn]$ at each time that we wish to estimate. The observation $[eqn]$ is modelled as a partial noisy observation of $[eqn]$ . This statement is summarised by the following model.

$$[eqn]$$

Here $[eqn]$ is a horizontal stacking of each partial graph observation such that the t-th column is the partial graph signal observed at time t, and the n-th row is the time series observed at node n. $[eqn]$ has a similar structure, but represents the true underlying function at the full, length-N, node-set at all times. $[eqn]$ is an error matrix of standard normal i.i.d. noise. The probability distribution of $[eqn]$ can therefore be expressed as

$$[eqn]$$

where $[eqn]$ has the usual meaning of a vertical stacking of matrix columns, and $[eqn]$ is the Kronecker product. In order to estimate the latent signal $[eqn]$ , we must specify a prior distribution indicating our belief at the outset about the likelihood of different signals. This provides the necessary Tikhonov regularization term and avoids under-specification. In this case, an appropriate prior for $[eqn]$ , which has also been used in [14], is the following

$$[eqn]$$

where matrix $[eqn]$ is a filtered graph Laplacian matrix, as defined in Equation (3), for a user specified graph filter $[eqn]$ , such as one of those specified in Table 2 and $[eqn]$ is a precision parameter which controls the regularisation strength.

Here, $[eqn]$ is a kernel (or Gram) matrix defined by the relation $[eqn]$ , where $[eqn]$ is a Mercer kernel (see [35]) with parameter(s) $[eqn]$ . In a kernel matrix $[eqn]$ , the entries represent the inner product between pairs of (potentially infinite-dimensional) basis function representations of the features. That is, $[eqn]$ . An expample of a typical kernel used in machine learning is the Gaussian kernel:

$$[eqn]$$

A natural question to ask is why we use the prior given in Eq (8), with a covariance matrix given by the Kronecker product between a kernel matrix and a squared graph filter. Intuitively, it makes sense that the correlation between two node-times $[eqn]$ and $[eqn]$ is expected to be high when the nodes are closely connected and the explanetory variables are similar, i.e. $[eqn]$ − $[eqn]$ is small. Eq (8) directly maps onto this intuition, since the prior covariance between two signal elements $[eqn]$ and $[eqn]$ is given by

$$[eqn]$$

In the supplementary materials, we provide a more formal justification for why this prior is appropriate. This is derived by considering the “weight-space" view of multivariate Gaussian processes [36]. This derivation begins with a multivariate linear regression model with a spherical prior over the regression coefficients. Then, the so-called ‘kernel trick’ in used to translate this into a non-parametric model with a kernel matrix $[eqn]$ . The only substantial modification to this original derivation that we make is to place a graph spectral prior over the regression coefficients, rather than the usual spherical prior. This encourages predictions which are smooth with respect to the graph topology.

Applying Bayes’ rule to Eqs (7) and (8) results in the following Maximum a Posteriori (MAP) optimisation problem for $[eqn]$ .

$$[eqn]$$

where we use the notation $[eqn]$ to represent a probability density function. Using the definition of the density function for a matrix normal distribution, $[eqn]$ can be written as

$$[eqn]$$

Taking the derivative of this expression with respect to $[eqn]$ and setting the result to zero yields a solution for $[eqn]$ . We refer to this as the KGR solution, but can also equivalently be considered the mean of a GPoG solution.

$$[eqn]$$

Applying a well-known matrix identity ([37], Eq 162) allows the dimension of the inverse to be reduced from $[eqn]$ to $[eqn]$ .

$$[eqn]$$

where $[eqn]$ and $[eqn]$ . However, in this form, computation is impractical for large problems (T or N or both large) due to the high complexity involved in inverting the Kronecker-structured matrix given by the form $[eqn]$ . Thankfully this can be overcome by performing eigendecomposition on $[eqn]$ and $[eqn]$ separately and leveraging properties of the Kronecker product. For a detailed derivation, consult the supplementary materials. Using the notation $[eqn]$ to represent the Hadamard product, the result is

$$[eqn]$$

where

$$[eqn]$$

and $[eqn]$ has the elements given by

$$[eqn]$$

A key benefit of the solution in this form is that the necessary dense eigendecompositions, which are typically the most computationally taxing step, are performed on $[eqn]$ and $[eqn]$ , of size $[eqn]$ and $[eqn]$ respectively, whereas a naive implementation would require this to be $[eqn]$ and $[eqn]$ . This can be a significant speed-up for large problems, especially when a meaningful portion of the data is unlabelled. (Decomposition of $[eqn]$ is also required, but this is typically sparse so can often be accelerated using sparsity-specific linear algebra tools, see [38] and references therein).

A full outline of all the steps for graph kernel regression is highlighted in algorithm 1.

$$[eqn]$$

3.3 KGR and Cartesian product graphs

In this short section, we highlight the connection between Kernel Graph Regression and Graph Signal Reconstruction (GSR). In particular, the mathematical formalism of KGR described so far in this paper bears a strong resemblance to Bayesian signal reconstruction as applied to a Cartesian product graph. The reconstruction of signals defined over Cartesian product graphs is an area that has received increasing attention in recent years [39,40], particularly in the context of Time-Vertex (T-V) problems [41]. By adapting the KGR algorithm slightly, we effectively get a signal reconstruction algorithm with little effort. This also highlights the connection between KGR and GSR which can provide insight into both areas.

A Cartesian product graph is described by two adjacency matrices $[eqn]$ and $[eqn]$ which, in the following, we refer to as the ‘time-like’ graph and the ‘space-like’ graph respectively. Fig 2 gives a graphical depiction of a small Cartesian product graph. The resultant graph has NT nodes, with an adjacency matrix $[eqn]$ , given by the Kronecker sum of the two individual adjacency matrices, specified as

A graphical depiction of a small Cartesian product graph.

$$[eqn]$$

Similarly, the Laplacian matrix of the Cartesian product, $[eqn]$ , is given by the Kronecker sum of the individual Laplacian matrices, that is

$$[eqn]$$

The Laplacian can be eigendecomposed as follows.

$$[eqn]$$

where $[eqn]$ and $[eqn]$ . Therefore, a general graph filter $[eqn]$ associated with the product graph can be understood as applying a decreasing function $[eqn]$ to the diagonal eigenvalue matrix $[eqn]$ .

$$[eqn]$$

For a certain class of graph filter functions, for example, the exponential filter, it is the case that

$$[eqn]$$

We refer to graph filters of this type as separable. If this condition holds, it implies that the total graph filter can also be expressed as a Kronecker product.

$$[eqn]$$

where $[eqn]$ and $[eqn]$ . Now consider a graph signal reconstruction problem, where the task is to estimate a smooth underlying signal $[eqn]$ , given a signal $[eqn]$ observed over a constant subset of the space-like and time-like nodes. Once again, we assume that $[eqn]$ is a noisy partial observation of $[eqn]$ given by Eq (6). The only key difference from the KGR formulation is that the prior distribution for $[eqn]$ is no longer given by Eq (8), but instead is given by

$$[eqn]$$

Note that the only material difference is that the kernel matrix $[eqn]$ has been replaced by another filter matrix $[eqn]$ . Following the Bayesian logic of the previous section immediately leads to a graph signal reconstruction method, given explicitly in algorithm 2.

$$[eqn]$$

4 GLS kernel graph regression

In many real-world applications of regression modelling, the assumption that the error terms of Eq (6) are independent and identically distributed is unlikely to hold. In this section, we consider the situation where the differences between the observed and underlying signals are instead correlated according to a matrix normal distribution. Under this model, the elements of $[eqn]$ are distributed as

$$[eqn]$$

where $[eqn]$ and $[eqn]$ respectively represent the time correlation and vertex covariance matrices. We choose to parameterise the covariance matrix in this way for several reasons. Firstly, some method for reducing the dimensionality of the problem is certainly necessary for the estimation of the matrix. Since the available data to make the estimate is a single $[eqn]$ observation $[eqn]$ , estimating a full $[eqn]$ matrix would not be possible without other regularisation assumptions. Secondly, expressing the covariance matrix as a Kronecker product of two smaller matrices enables tractable estimators. This is primarily due to the property of Kronecker products that $[eqn]$ . In addition, holding a dense $[eqn]$ matrix in memory may be out of the question. Thirdly, the Kronecker product assumption has an intuitive interpretation in terms of independent factor-specific contributions to the overall covariance. Put simply, the covariance between node-times $[eqn]$ and $[eqn]$ will be $[eqn]$ $[eqn]$ $[eqn]$ . If there is no correlation between two times t1 and t2, then there should be no correlation between node-times $[eqn]$ and $[eqn]$ , no matter what the value of $[eqn]$ . Similarly, if the two nodes are totally uncorrelated, then there should be no correlation between node-times $[eqn]$ and $[eqn]$ , no matter what the value of $[eqn]$ . The actual value, $[eqn]$ $[eqn]$ $[eqn]$ can be seen as counting the contribution from both the time correlation and the node-level correlation multiplicatively. Finally, this model has also been used in many applications relevant to network problems such as geospatial, econometric, and EEG models [42–44]. For other recent approaches, see [45,46].

4.1 KGR with Gauss–Markov estimator

To begin, we derive a Gauss–Markov estimator for $[eqn]$ , that is, the Best Linear Unbiased Estimator (BLUE) assuming that the covariance matrix $[eqn]$ $[eqn]$ $[eqn]$ is known. (This restriction is relaxed in the subsequent section). In this case, the log-likelihood of making an observation $[eqn]$ is altered such that the cost function of Eq (10) becomes

$$[eqn]$$

The BLUE estimator for a given $[eqn]$ and $[eqn]$ is defined to be the value of $[eqn]$ which minimises this expression and can be found by differentiating it with respect to $[eqn]$ and setting the result equal to zero. By again following the same steps as in section 3, this procedure results in

$$[eqn]$$

Once again, in this form, the solution remains prohibitively expensive to compute for large $[eqn]$ and $[eqn]$ . However, it can be expressed alternatively with $[eqn]$ complexity by making the following definitions. First eigendecompose $[eqn]$ and $[eqn]$ .

$$[eqn]$$

Then perform eigendecomposition on the following matrices:

$$[eqn]$$

Note that both these matrices are guaranteed to be symmetric with positive, real eigenvalues [47]. Finally, make the following definitions:

$$[eqn]$$

The GLS estimator for $[eqn]$ is given by

$$[eqn]$$

While a simpler O(N^3^  +  T^3^) solution is possible, the above method ensures that the eigendecompositions are performed on reduced size matrices, effectively scaling as $[eqn]$  +  $[eqn]$ . For a detailed derivation with intermediate steps, we refer the reader to the supplementary materials.

$$[eqn]$$

4.2 Estimating ΣNSigN and ΣTSigT given 𝐅F

In this section, we construct estimators for the covariance matrices $[eqn]$ and $[eqn]$ , given that the latent signal $[eqn]$ , and therefore the prediction error $[eqn]$ , is known. It follows from Eq (23), that $[eqn]$ has a matrix normal distribution of the following form.

$$[eqn]$$

Past literature on parameter estimation in matrix-normal models has mostly focused on the task of estimating the covariance given multiple realisations of the error matrix $[eqn]$ [48–51]. In order to estimate both $[eqn]$ and $[eqn]$ in full, the number of observations of $[eqn]$ must be greater than [max $[eqn]$  +  1] [52]. However, in the present case, only a single observation comprising $[eqn]$ values is available. It is therefore necessary to assume a simplified, parametrised structure for at least one of these matrices to reduce the number of degrees of freedom.

For this reason, we assume that the matrix $[eqn]$ describes the autocorrelation between different time measurements and is a function of a parameter or set of parameters $[eqn]$ . The elements of $[eqn]$ are constrained to be dimensionless Pearson correlation coefficients, with 1 along the diagonal. Conversely, $[eqn]$ is assumed to be an arbitrary covariance matrix that can take on any value satisfying semi-positive definiteness. By explicitly parametrising one as a correlation matrix, parameterised by a function of a parameter vector $[eqn]$ in much lower dimension than T, while letting the other act as a fully parameterised covariance matrix, the problem of the solution being under-determined (since $[eqn]$ gives the same matrix for any a) is avoided.

The likelihood of observing an error matrix $[eqn]$ , given covariance matrices $[eqn]$ and $[eqn]$ , is

$$[eqn]$$

By specifying two independent prior distributions $[eqn]$ and $[eqn]$ , the posterior density over $[eqn]$ and $[eqn]$ can be expressed via Bayes’ rule as

$$[eqn]$$

Depending on the size of the problem, this could be used along with a Monte Carlo algorithm to sample from the whole posterior or, more likely, be used to calculate single MAP estimates for $[eqn]$ and $[eqn]$ . In either case, the negative log posterior, up to a scaling and additive constant, is equal to

$$[eqn]$$

A MAP estimate can be found by minimising this quantity jointly with respect to the parameters. There are several possible strategies for solving this, such as simple gradient descent. However, since an analytical solution for $[eqn]$ is often possible given $[eqn]$ and vice versa, an efficient algorithm using temporary approximate solutions can often be found. This so-called ‘flip-flop’ strategy is widely used in the literature on Kronecker-product covariance estimation [52]. Crucially, the number of iterations required for convergence does not depend greatly on the size of the system, meaning $[eqn]$ complexity is retained.

4.3 Example: AR(1) process

The previous section stated the general formulation of the problem without any particular parametrisation of the matrix $[eqn]$ . In this section, we give a concrete example by considering the widely used AR(1) model, which assumes simple, serial correlation between prediction errors. Denoting the t-th column of $[eqn]$ as $[eqn]$ , the general AR(1) model assumes that

$$[eqn]$$

where $[eqn]$ is a matrix of autoregression coefficients [53]. In the following, we assume a simplified version, where the n-th element of $[eqn]$ is assumed only to be serially correlated with the n-th element of $[eqn]$ , in a way that is both stationary and uniform across all n. In effect, this requires that $[eqn]$ . This is, in essence, the multivariate extension of the AR(1) model considered in the seminal paper by Cochrane and Orcutt [54]. In the limit as $[eqn]$ , the correlation matrix $[eqn]$ , and its inverse $[eqn]$ , have a well-known form, which can be truncated approximately for large enough T to the form given by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where $[eqn]$ is a matrix with 0 on the diagonal and with 1 on the first upper and lower diagonal, and $[eqn]$ is the identity matrix with zero on the first and last diagonal entries [55]. Note also that the determinant of $[eqn]$ is given by $[eqn]$ [56]. To be a valid stationary process, the parameter $[eqn]$ must lie on the interval (–1,1). An appropriate prior distribution for $[eqn]$ could therefore be

$$[eqn]$$

for some scalar parameter $[eqn]$ . Depending on the value of $[eqn]$ , this prior is roughly uniform across the majority of the interval, whilst rapidly decreasing in likelihood towards –1 and 1. This effectively encodes a belief that the time series is stationary across all nodes, with $[eqn]$ controlling the strength of that stationarity assumption. The resultant $[eqn]$ -dependant part of the cost function of Eq (29) therefore becomes

$$[eqn]$$

Whilst this is not a quadratic function of $[eqn]$ , the unique maximum likelihood estimator $[eqn]$ can be found by differentiating this expression and setting the result to zero. This results in a MAP estimator which is the real root of a cubic polynomial.

$$[eqn]$$

where

$$[eqn]$$

A derivation of this expression can be found in the supplementary materials. Note that this estimator is a function of $[eqn]$ .

In terms of the cross-covariance matrix $[eqn]$ , we choose to implement a modified form of the estimator proposed in [57] and further developed in [58]. In essence, these papers considered a weighted combination of the high-variance and often ill-posed maximum likelihood estimate, $[eqn]$ , and the heavily-biased but well-conditioned estimate, $[eqn]$ . The only modification presented here is the introduction of $[eqn]$ into the ML estimate, that is, $[eqn]$ , which is necessary to account for the effect of autocorrelation.

$$[eqn]$$

[57] and [58] also investigate the optimal setting of the shrinkage coefficient $[eqn]$ , which is a parameter between 0 and 1. Here, we choose to implement the Rao-Blackwell Ledoit-Wolf (RBLW) estimator described in [58]. This is given by

$$[eqn]$$ $$[eqn]$$

[Eq (34)]

4.4 GLS kernel graph regression

To complete the GLS kernel regression algorithm, estimation of both $[eqn]$ and $[eqn]$ $[eqn]$ $[eqn]$ must be performed simultaneously. Since both rely on each other to be estimated, it is necessary to implement an alternating iterative algorithm, where first $[eqn]$ is estimated for some reasonable initial guess of $[eqn]$ and $[eqn]$ , and then $[eqn]$ and $[eqn]$ are estimated for this value of $[eqn]$ . This process then continues until some convergence criterion is met. These steps are outlined in algorithm 5.

Estimation of the hyperparameters of the kernel is performed using a grid search procedure as discussed further in Section 7.1 below. Additional details of this hyperparameter estimation approach are available in the github repository for the code and detailed worked examples on synthetic and real data, including those contained in this work and additional real data examples, see https://github.com/nickelnine37 and further discussions in [38]. In addition, in the following Section 4.5 explores the model sensitivity to various hyperparameters.

4.5 Investigating model sensitivity to hyperparameter selection

In this subsection, we study the effect that the various hyperparameters have on the properties of the GLS-KGR algorithm, using synthetic data. The aim is to provide some intuition for how these variables should be set or learned on real data. The key variables of interest are γ, which acts as a global regularisation parameter, β, which dictates how smooth the graph signals are expected to be, α, which controls the strength of the stationarity assumption, and σ, which affects the local variance scale of the Gaussian kernel. We are interested in two key effects that these hyperparameters have on the algorithm, namely, their effect on the prediction accuracy, and their effect on the convergence rate.

$$[eqn]$$

[algorithm 3]

[algorithm 4]

Each experiment was set up in the following way. First, a smooth underlying graph signal $[eqn]$ was generated by applying a chain-graph diffusion filter $[eqn]$ to a matrix of Gaussian noise. 20% of the times and nodes were chosen uniformly at random to be hidden, giving two selection matrices $[eqn]$ and $[eqn]$ . A model error matrix $[eqn]$ was then generated from a matrix normal distribution with known covariance matrices $[eqn]$ and $[eqn]$ , and added to $[eqn]$ to create the observed signal $[eqn]$ . Finally, a matrix of Gaussian noise with 12 columns $[eqn]$ was created to act as the explanatory variables. This model was run under a variety of hyperparameter conditions, with the number of steps required for convergence, and the Root Mean Squared Error (RMSE) on both the seen and unseen node-times recorded on each run.

First, we analyse the convergence properties. The experimental variables $[eqn]$ , $[eqn]$ , and $[eqn]$ were found to have a negligible effect on the convergence rate of the GLS Kernel Graph Regression algorithm. That is, when $[eqn]$ , $[eqn]$ and $[eqn]$ remained fixed, varying any of the aforementioned experimental parameters did not affect the number of iterations required to reach a minimum level of precision across the model outputs $[eqn]$ and $[eqn]$ .

On the other hand, we found that $[eqn]$ , $[eqn]$ and $[eqn]$ all had a significant effect on convergence. In order to analyse this systematically, we fixed N = 100, T = 120, $[eqn]$ , $[eqn]$ and $[eqn]$ . We then varied $[eqn]$ from $[eqn]$ to $[eqn]$ in 9 increments of 0.2, and varied $[eqn]$ from 10^−3^ to 10^1^ in 50 logarithmic-spaced increments. For each unique $[eqn]$ pair we generated 50 unique sets of input data $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ randomly as specified at the start of this section, and ran the GLS KGR model to give $[eqn]$ unique trials. For each trial, we counted the total number of iterations required for convergence, including both the inner loop detailed in algorithm 4 and the outer loop of algorithm 5. This experiment was performed twice: In experiment (a) we set $[eqn]$ , and in experiment (b) we set $[eqn]$ , where $[eqn]$ is the Kronecker-delta symbol such that $[eqn]$ iff i = j and $[eqn]$ otherwise.

The mean number of iterations required for convergence for each $[eqn]$ pair is shown in Fig 3.

The mean number of iterations required for convergence on synthetic data is shown over a range of values for θ and α.The upper plot shows the results of experiment (a), when ΣN=𝐈, and the lower plot shows the results for experiment (b), when ΣNij=exp((i−j)2)+0.05δij.

In general, we can see that $[eqn]$ , θ and α have a complex interaction with the number of iterations. In both experiments (a) and (b) we see that convergence is roughly symmetric for $[eqn]$ , with this symmetry slightly broken in experiment (b). In both cases, a spike in the number of iterations is observed for a certain value of α, although this is not observed in the experimental α range for $[eqn]$ in experiment (a). Further investigation into this interaction from a theoretical standpoint would be valuable, however, we leave this to a future study.

Next, we studied the effect that the hyperparameters had on prediction accuracy across the test and train sets. Here, by test set, we mean all unobserved node-times, and by train set, we mean all observed node-times. Since 80% of the nodes and 80% of the times were observed, this means that $[eqn]$ of all node-times were in the train set, and 36% were in the test set. In this case, the most important parameters were γ, β and σ, with $[eqn]$ having a negligible impact on the prediction accuracy.

First, we used gradient descent to find the optimal parameters, with an objective function that evaluated the mean test set RMSE across 50 random realisations of the input data. The optimal values, in this case, were found to be $[eqn]$ , $[eqn]$ and $[eqn]$ . Next, we fixed two of these parameters and varied the third across a range of values, measuring the mean test and training set RMSE across 50 random realisations of the input data. The results are shown in Figs 4, 5, 6, 7. As is visible, the hyperparameter with the largest effect, and therefore the most important to set appropriately, is γ. This is expected, since it acts as a global regularisation parameter. β and σ also have a small but noticeable effect.

The Root Mean Squared Error on the test and train sets are shown as γ is varied.

The Root Mean Squared Error on the test and train sets are shown as β is varied.

The Root Mean Squared Error on the test and train sets are shown as σ is varied.

The Root Mean Squared Error on the test and train sets are shown as a function of the total percentage of node-times that were observed.

As a final experiment, we also measured the accuracy across both test and train sets as a function of the percentage of nodes and times that were observed. Here, we ran the experiment as before, with the optimal hyperparameters, but increased the total percentage of node-times that were observed from p = 0.05 to p = 0.99 in 50 equally spaced increments. On each run, we set $[eqn]$ and $[eqn]$ , rounded to the nearest integer.

5 Latent signal uncertainty via the Laplace approximation

In this final methodology section, we outline the steps for estimating the uncertainty over the latent signal $[eqn]$ via the Laplace approximation. In the case of KGR, this is in fact not an approximation, but the exact Gaussian process prior. On the other hand, in the GLS case, it is indeed an approximation as point estimates for the uncertain covariance matrices $[eqn]$ and $[eqn]$ are used. Here, we derive the approximate posterior for the GLS case only but note that the exact posterior for simple KGR can be restored by setting $[eqn]$ and $[eqn]$ equal to the identity matrix of an appropriate size.

Since the latent signal $[eqn]$ has size $[eqn]$ , the true covariance matrix specifying its uncertainty has size NT $[eqn]$ NT. Deriving a mathematical expression for this matrix is simple, however, it is generally impractical to compute and store in memory for all but the smallest of problems. For this reason, we take steps to derive an efficient expression for the $[eqn]$ matrix representing the marginal variance of each element of $[eqn]$ , which is generally the most useful information and allows scaling to relatively large problems.

Consider the negative log-likelihood given in Eq (24). Applying Laplace’s approximation gives the covariance matrix for the uncertainty over the vectorized latent signal.

$$[eqn]$$ $$[eqn]$$

Holding this dense matrix in memory is often impractical, and inverting it is generally out of the question. Instead, we introduce the matrix $[eqn]$ defined by $[eqn]$ . That is, $[eqn]$ is the diagonal of $[eqn]$ of length NT stacked into an $[eqn]$ matrix of marginal variances, such that element (n,t) is the square uncertainty for the n-th node at time t. Here we present an expression for $[eqn]$ followed by a short derivation.

$$[eqn]$$

where $[eqn]$ is the Hadamard product,

$$[eqn]$$

and

$$[eqn]$$

Note that the eigendecompositions are no longer necessarily performed on symmetric matrices. This means that while $[eqn]$ and $[eqn]$ cannot be assumed to be unitary, the eigenvalues are still guaranteed to be positive and real [47]. To prove this expression, first note that $[eqn]$ can be factorized into the following:

$$[eqn]$$

From here, substitute in the eigendecomposition definitions.

$$[eqn]$$

This can be split into the following sum.

$$[eqn]$$

where $[eqn]$ is defined to be a $[eqn]$ matrix of zeros with element (t,t) equal to one. From this, it can be seen that the re-stacked diagonal $[eqn]$ is given by the following outer product sum.

$$[eqn]$$

where

$$[eqn]$$

Or more compactly

$$[eqn]$$

Further detail concerning this derivation can be found in the supplementary materials.

6 Experimental results

6.1 Spatio-temporal pollutant analysis

In this section, we consider the problem of predicting the concentration of various airborne pollutants, measured across a network of $[eqn]$ air quality monitoring stations in and around California. Each pollutant type is measured at a unique set of locations and is therefore treated as an independent graph regression task. The goal is to make accurate predictions using weather and environmental features from the previous day. The methods of Kernel Graph Regression and AR(1) GLS Kernel Graph Regression are analysed and compared to some standard baseline algorithms.

Pollutant concentration data was taken from the US Environmental Protection Agency’s air quality monitoring program [59]. Specifically, daily measurements of Ozone, Carbon Monoxide (CO), Nitrogen Dioxide (NO_2_), PM_2.5_ and PM_10_ were taken from January 2017 to April 2021, giving a total of T = 1570 days. This data set also contains daily measurements of humidity, pressure, wind speed and temperature at various locations, which we use as additional explanatory variables. In addition, data concerning historical wildfires in California was sourced from the Department of Forestry and Fire Protection in California [60].

6.2 Graph construction

A key decision that can significantly impact the effectiveness of a graph signal processing method is how to construct the underlying graph. In certain applications, such as social networks, a sparse graph may be self-evident or relatively simple to construct. However, for a large portion of practical problems, including the current case of geographically placed sensors, it is required to either propose a sensible construction method or learn a graph from available signal data. In this paper we opt for the former, and omit a full discussion of the available techniques, which can be found in, for example, [61].

The graph construction method we use makes use of both pairwise geodesic distances between monitors and the intermediate elevation profile. This is important as environmental processes can be strongly influenced by topography, especially in mountainous regions. Elevation data is sourced from the GLOBE30 project [62]. This dataset gives the approximate height above sea level over a 30 arc-second latitude/longitude grid. While more complex models incorporating land use, prevailing wind direction etc. are possible, we opt for a simpler model for the sake of brevity.

Our first step is to create an $[eqn]$ symmetric distance matrix $[eqn]$ . We define the “distance" between two monitors $[eqn]$ to be a weighted combination of their geodesic distance and the vertical relief along the intermediate path. This introduces a hyperparameter defining the relative importance of each component which we learn later via cross-validation. We then use the perturbed minimum spanning tree algorithm outlined in [63] to construct a sparse, fully connected graph. A representation can be found in Fig 8.

A graphical depiction of the graph connecting the 1382 available monitoring stations, with the map colour representing the elevation profile.This is the graph used in all further experiments. “Reprinted from [38] under a CC BY license, with permission from [Heriott-Watt University], original copyright [2024].” Source available at: http://hdl.handle.net/10399/5041

6.3 Data pre-processing

Several steps were taken to pre-process the raw sensor data before feeding it into the models. Firstly, any stations with more than a total of 15 null readings due to equipment failure, over the three-year period, were discarded. Any remaining missing values were interpolated linearly between the available readings. Secondly, meter readings for each of the pollutant variables, which generally have units of mass or particle count per unit volume, were transformed via $[eqn]$ . This rectified the left skew of the readings, generally producing more symmetric aggregate histograms. The same transformation was also applied to wind speed, however pressure, humidity and temperature, which already had largely symmetric data distributions, were not transformed in this way.

To construct a daily set of features for wildfires, a historical list of all recorded wildfire events in California was used. This dataset includes information about the start and end date of each fire, along with the central location and the total ground area burned. For each fire we assumed that the burn rate rose and fell with a Gaussian shape, peaking at the midpoint of the burn window, such that 95% of the area burns within the specified time range. From this we estimate the total area burning on any given day in eleven different regions in California.

The final step before estimation was to transform all variables so that they were scaled and translated to achieve a unit marginal variance and zero mean. PCA dimensionality reduction was then performed on each of the feature groups individually keeping enough dimensions to preserve 90% of the total variance for that group. The resultant input data for $[eqn]$ and $[eqn]$ had 1570 columns representing the dates from 2017-01-02 to 2021-04-20. Each of the N rows of $[eqn]$ held a time series of readings for a unique monitoring sensor with a well-defined location. The columns of $[eqn]$ contain first the PCA components of the weather conditions, second the approximate number of acres burning in each of the eleven regions. Regression was then performed for each of the five different pollutants. Table 3 gives information about the total number of features that were abvailable after perfroming these processing steps.

Table 3: Information about the number of available features before and after PCA.

6.4 Results

To test the performance of the KGR and GLS KGR methods, we compared them against several other regression algorithms, namely Ordinary Least Squares (OLS) regression, Ridge regression, Lasso regression and Elastic Net regression.

For each model tested, all relevant hyperparameters were tuned using cross-validation. First, the input data was split such that the first 80% of the days served as a training/validation set and the final 20% of the days as a test set. The train/validation set was then split into four folds of 219 days each uniformly at random. Hyperparameters were set by attempting to minimise the mean squared error averaged across each of these four folds, using three to train and one to validate accordingly. Final results were then reported on the test set. Fig 9 shows the predicted Ozone levels across the full graph (right) on a particular day given the partially observed signal (left). As expected, the output is fairly smooth across the network, indicating that the model has successfully utilised information from the graph structure. Fig 10 shows the output from the GLS KGR algorithm at a particular unobserved node. As is visible, the estimated reading approximates the ground truth well, especially during train times. In addition, the uncertainty about the estimate notably increased during the test times, further indicating a healthy model output.

Left: the observed graph signal for Ozone on a particular day represented via a colour map.Right: prediction for the reconstructed latent graph signal across the entire network on the same day made by the GLS KGR method. “Reprinted from [38] under a CC BY license, with permission from [Heriott-Watt University], original copyright [2024].” Source available at: http://hdl.handle.net/10399/5041

A section of the Ozone time-series signal predicted by GLS KGR at a particular unobserved node is shown along with the ground truth.The blue shading depicts two standard deviations of prediction error arising from latent signal uncertainty as calculated by the Laplace approximation.

Table 4 shows the mean squared error as reported on the test set after hyperparameter tuning. As is visible, either GLS KGR or KGR performs the best by this metric across all pollutants. Table 5 shows the total compute time for each method on a 4-core i7 intel CPU.

Table 4: Mean squared Error on all unseen data.

Table 5: The run time for each algorithm in seconds.

7 Discussion

7.1 Hyperparameter tuning

One weakness with KGR and GLS KGR in practice is that several hyperparameters must be tuned to produce an accurate model, namely $[eqn]$ and $[eqn]$ , as well as any used in the graph construction period. Each is important and can have a significant effect on the output produced. For example, insufficient regularisation via $[eqn]$ can cause severe over-fitting and result in unduly small entries within the estimated cross-covariance matrix $[eqn]$ . Similarly, setting $[eqn]$ too high can effectively drive the estimated value of $[eqn]$ to zero, while setting it too low can result in singular matrices if the data is at all non-stationary. This creates a large, non-convex search space that can be both costly to explore and full of local minima. The degree to which this is an issue depends largely on the size of the problem at hand. For relatively small problems, such as the one considered here with $[eqn]$ nodes and $[eqn]$ time steps, well-known scalar minimization algorithms such as Nelder-Mead can help automate this process. However, with significantly larger problems hyperparameters would have to be carefully selected based on domain knowledge.

7.2 The normalised graph Laplacian

One model choice not discussed so far is the question of whether to use the normalised or un-normalised version of the graph Laplacian. Until now, we have assumed the regular Laplacian, $[eqn]$ , is being used, however, the normalised version, defined by $[eqn]$ , has properties that may make it preferable in some circumstances. For example, graph penalties constructed using $[eqn]$ give equal weight to all nodes, whereas penalties constructed using $[eqn]$ implicitly favour high-degree nodes because they appear more frequently in the sum of Eq (1) [27]. This may or may not be desirable depending on the problem at hand. Another potential benefit of $[eqn]$ is that its eigenvalues are guaranteed to fall in the interval [0,2] which makes the graph-regularisation parameter $[eqn]$ easier to set and interpret, as well as making a sensible comparison across problems possible. In the case of environmental monitoring networks, our initial experiments indicate that the normalized Laplacian may result in slightly better performance, however, further investigation is necessary.

7.3 Options for increasing scalability

In the previous analysis, we assumed that, while the graph Laplacian $[eqn]$ is often sparse, the filter and kernel matrices $[eqn]$ and $[eqn]$ were dense. This typically means that the primary bottleneck for KGR and GLS KGR is the eigendecomposition of these matrices, even in their down-sampled form. For large problems, it may be sufficient to only calculate the first k eigenvectors and eigenvalues, which can be performed efficiently for sparse matrices. Therefore finding sparse representations for $[eqn]$ and $[eqn]$ can be highly beneficial for large applications. Numerous methods exist in the literature on kernel regression for improving computational robustness. One, in particular, is the Nystrom method, which uses a low-rank approximation to the kernel matrix [64]. This can reduce the computational complexity of solving the linear system to $[eqn]$ , and the memory requirements to O(M^2^ + MN), where M is the chosen number of data examples.

One way to create a sparse graph filter is to simply use a low-degree polynomial for the filter function $[eqn]$ . In this way, $[eqn]$ can be calculated efficiently as $[eqn]$ and decomposed faster. When it comes to the kernel matrix, a large literature also exists on compactly supported kernels, for example, the Wendland kernel [65]. Here, an integral operator is applied recursively to Askey’s truncated power functions to create a set of compactly supported radial basis functions, which are guaranteed to be smooth, continuous and positive definite. Fast algorithms for the computation of these basis functions have been outlined in [66].

8 Conclusions and future work

This paper has contributed a statistical analysis of regression methods for signals defined over networks. Drawing upon recent work on Kernel Graph Regression, Gaussian Processes over Graphs and graph signal reconstruction, we proposed the method of GLS KGR and demonstrated its effectiveness on a relevant task. In particular, we addressed the situation where a partial observation of an N-node graph signal is made at a subset of T time-points, deriving the steps for an AR(1) autocorrelation regression model, with general cross-correlation. By assuming a matrix-normal error distribution, an algorithm was designed with O(N^3^  +  T^3^) complexity at each iteration. Finally, the Laplace approximation was used to derive a lower bound for the marginal prediction error arising from latent signal uncertainty. This was tested and shown to be effective on real data taken from a network of pollutant monitoring stations.

One assumption in this work was that the unobserved nodes remain so over the lifetime of the problem. This is somewhat unrealistic in applications such as sensor networks where equipment may temporarily fail. In the future, it would be valuable to work on adapting these algorithms for arbitrary missing values. Another worthwhile investigation would be a systematic study into the effect of different graph construction methods within the context of graph regression.