A Gaussian Process Method with Uncertainty Quantification for Air Quality Monitoring

Abstract

The monitoring and forecasting of particulate matter (e.g., ${\mathrm{PM}}_{2.5}$) and gaseous pollutants (e.g., NO, ${\mathrm{NO}}_2$, and ${\mathrm{SO}}_2$) is of significant importance, as they have adverse impacts on human health. However, model performance can easily degrade due to data noises, environmental and other factors. This paper proposes a general solution to analyse how the noise level of measurements and hyperparameters of a Gaussian process model affect the prediction accuracy and uncertainty, with a comparative case study of atmospheric pollutant concentrations prediction in Sheffield, UK, and Peshawar, Pakistan. The Neumann series is exploited to approximate the matrix inverse involved in the Gaussian process approach. This enables us to derive a theoretical relationship between any independent variable (e.g., measurement noise level, hyperparameters of Gaussian process methods), and the uncertainty and accuracy prediction. In addition, it helps us to discover insights on how these independent variables affect the algorithm evidence lower bound. The theoretical results are verified by applying a Gaussian processes approach and its sparse variants to air quality data forecasting.

1. Introduction

It is generally believed that urban areas provide better opportunities in terms of economic, political, and social facilities compared to rural areas. As a result, more and more people are migrating to urban areas. At present, more than fifty percent of people worldwide live in urban areas, and this percentage is increasing with time. This has led to several environmental issues in large cities, such as air pollution [1].

Landrigan reported that air pollution caused $6.4$ million deaths worldwide in 2015 [2]. According to World Health Organization (WHO) statistical data, three million premature deaths were caused by air pollution worldwide in 2012 [3]. Air pollution has a strong link with dementia, causing 850,000 people to suffer from dementia in the UK [4]. Children growing up in residential houses near busy roads and junctions have a much higher risk of developing various respiratory diseases, including asthma, due to high levels of air pollution [5]. Polluted air, especially air with high levels of NO, ${\mathrm{NO}}_2$, and ${\mathrm{SO}}_2$ and particulate matter (PM${}_{2.5}$), is considered the most serious environmental risk to public health in urban areas [6]. Therefore, many national and international organisations are actively working on understanding the behaviour of various air pollutants [7]. This eventually leads to the development of air quality forecasting models so that people can be alerted in time [8].

Essentially, being like a time series, air quality data can be easily processed by models that are capable of time series data processing. For instance, Shen applies an autoregressive moving average (ARMA) model in PM${}_{2.5}$ concentration prediction in a few Chinese cities [9]. Filtering techniques like Kalman filter are also applied to adjust data biases to improve air quality prediction accuracy [10]. These methods, though with good results reported, are limited by the requirement of a prior model before data processing. Machine learning methods, on the other hand, can learn a model from the data directly. This has enabled them to attract wide attention in recent decades in the field of air quality forecasting. For instance, Lin et al. propose the support vector regression with logarithm preprocessing procedure and immune algorithms (SVRLIA) method, which outperforms general regression neural networks (GRNN) [11] and BackPropagation neural networks (BPNN) [12] in Taiwan air quality forecasting [13].

Recently, inspired by the fact that large scale data are accumulated, deep learning models have been applied in air quality prediction [14]. Some work has added these deep learning models with the ability to quantify uncertainties introduced by inputs. For instance, Garriga-Alonso et al. endow a deep convolutional network with uncertainty quantification, by taking it as an equivalent of a Gaussian processes (GPs) model [15]. This is because GPs predictions are accompanied by confidence intervals, which are usually taken as a metric to measure prediction uncertainties. Applications of GPs in air quality forecasting can be found in [16,17]. However, the involvement of matrix inversion in GPs limits their application in large-scale datasets [18]. This has inspired research on improving the efficiency of GP models, and a series of efficient GP models have been published [19]. We also proposed an efficient GP model with application in air quality forecasting [17]. Despite the rich number of GP models published, there lacks work that investigates how noise level, hyperparameters, etc. affect the performance of GP models. It is necessary because air quality data vary due to seasonal variations and sensor degradations. A well-trained GP model may not work when fed with new data, simply due to measurement noise level change. By knowing how the variation of GPs performance can be attributed to noise level and hyperparameters, etc., we will still be able perform analysis when noise level or hyperparameters vary.

Aiming at this, a general solution is proposed in this paper. It provides insights on how a GP model’s performance is related to measurement noise level and hyperparameters, etc. The main contribution of this work includes (1) a general method for analysing how noise level and hyperparameters of a GP model affect the prediction performance. The variation of the evidence lower bound (ELBO) and the upper bound of the marginal likelihood (UBML) with respect to the noise level and hyperparameters are also given. (2) Neumann series is exploited to approximate the matrix inversion involved in GPs. This helps construct an analytical relation between noise level, hyperparameters, etc., and model performance. (3) A comparative air quality forecasting study between Sheffield, UK, and Pershawar, Pakistan is given, demonstrating that the proposed solution is able to capture how noise level and hyperparameters affect GPs performance.

The remaining part of this paper is as follows. Section 2 provides the theoretical fundamentals involved in this paper; Section 3 elaborates the proposed uncertainty quantification solution. In Section 4, we provide a comparative study of air quality prediction in the same period between the British city Sheffield and Pakistani city Pershawar, and the paper is concluded in Section 5. Appendix A describes the data collection process in Peshawar, Pakistan, and in Sheffield, United Kingdom, and presents maps of the considered areas of these cities. Appendix B gives the World Health Organisation (WHO) criteria for air pollutants. Appendix C gives the approximate derivatives of the GP kernel.

2. Background Knowledge

2.1. Gaussian Processes

Given a set of training data $\mathcal{D}=\{({\mathbf{x}}_i,y_i),i=1,\cdots ,n\}$ where ${\mathbf{x}}_i\in \mathcal{X}$ is the input and $y_i\in \mathbb{R}$ is the observation, we can determine a GP model $f(·)$ to predict $y_{*}$ for a new input ${\mathbf{x}}_{*}$. For instance, when the output is one-dimensional, the GP model is formulated as

$$f\sim \mathcal{GP}\left(\overline{f}\left(\mathbf{x}\right),k(\mathbf{x},{\mathbf{x}}^{\prime })\right),y=f\left(\mathbf{x}\right)+\epsilon ,\epsilon \sim \mathcal{N}(0,{\sigma }^2),$$

(1)

where $\overline{f}:\mathcal{X}\to \mathbb{R}$ is the mean function defined as

$$\overline{f}\left(\mathbf{x}\right)=\mathbb{E}\left[f\left(\mathbf{x}\right)\right],$$

(2)

and $k:\mathcal{X}\times \mathcal{X}\to \mathbb{R}$ is the kernel function [18] defined as

$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\mathbb{E}\left[(f\left(\mathbf{x}\right)-\overline{f}\left(\mathbf{x}\right))(f\left({\mathbf{x}}^{\prime }\right)-\overline{f}\left({\mathbf{x}}^{\prime }\right))\right],$$

(3)

where $\epsilon$ is the additive, independent, identically distributed Gaussian measurement noise with variance ${\sigma }^2\ne 0$, and $\mathbb{E}$ denotes the mathematical expectation operation.

Given ${\mathbf{x}}_i$ a $D\times 1$ vector, the n inputs can be aggregated into a matrix ${\mathbf{X}}_{D\times n}$, or briefly $\mathbf{X}$ with the corresponding output vector ${\mathbf{y}}_{n\times 1}$, or $\mathbf{y}$. Similarly, the function values at the test inputs ${\mathbf{X}}_{*}$ with dimensions of $D\times N$ can be denoted as ${\mathbf{f}}_{*}$, and we next write the joint distribution of $\mathbf{y}$ and ${\mathbf{f}}_{*}$ as

$$\left[\begin{array}{c}\mathbf{y}\\ {\mathbf{f}}_{*}\end{array}\right]\sim \mathcal{N}\left(\mathbf{0},\left[\begin{array}{cc}{\mathbf{K}}_{nn}+{\sigma }^2\mathbf{I}& {\mathbf{K}}_{nN}\\ {\mathbf{K}}_{Nn}& {\mathbf{K}}_{NN}\end{array}\right]\right),$$

(4)

where $\mathbf{I}$ represents the identity matrix. ${\mathbf{K}}_{nn}+{\sigma }^2\mathbf{I}$ is the $n\times n$ prior covariance matrix of $\mathbf{y}$ with entry ${\mathbf{K}}_{ij}=k({\mathbf{x}}_i,{\mathbf{x}}_j)+{\sigma }^2{\delta }_{ij}$, where ${\delta }_{ij}$ is one iff $i=j$ and zero otherwise, and ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ are column vectors from $\mathbf{X}$. The matrix ${\mathbf{K}}_{NN}$ denotes the $N\times N$ prior covariance matrix of ${\mathbf{f}}_{*}$ with entry ${\mathbf{K}}_{ij}=k({\mathbf{x}}_i,{\mathbf{x}}_j)$, where ${\mathbf{x}}_i$ and ${\mathbf{x}}_j$ are column vectors from ${\mathbf{X}}_{*}$. The matrices ${\mathbf{K}}_{Nn}$ and ${\mathbf{K}}_{nN}$ satisfy ${\mathbf{K}}_{Nn}={\mathbf{K}}_{nN}^{\mathrm{T}}$, and the entry of the $N\times n$ prior covariance matrix of ${\mathbf{f}}_{*}$ and $\mathbf{y}$ is ${\mathbf{K}}_{ij}=k({\mathbf{x}}_i,{\mathbf{x}}_j)$, where ${\mathbf{x}}_i$ is a column vector from ${\mathbf{X}}_{*}$ and ${\mathbf{x}}_j$ is a column vector from $\mathbf{X}$.

By deriving the conditional distribution of ${\mathbf{f}}_{*}$ from (4), where the prior mean is set to be zero for simplicity [20], we have the predictive posterior at new inputs ${\mathbf{X}}_{*}$ as

$${\mathbf{f}}_{*}|\mathbf{X},\mathbf{y},{\mathbf{X}}_{*}\sim \mathcal{N}\left({\overline{\mathbf{f}}}_{*},\mathrm{cov}\left({\mathbf{f}}_{*}\right)\right),$$

(5)

where

$${\overline{\mathbf{f}}}_{*}\triangleq \mathbb{E}\left[{\mathbf{f}}_{*}|\mathbf{X},\mathbf{y},{\mathbf{X}}_{*}\right]={\mathbf{K}}_{Nn}{\left[{\mathbf{K}}_{nn}+{\sigma }^2\mathbf{I}\right]}^{-1}\mathbf{y},$$

(6)

is the prediction at ${\mathbf{X}}_{*}$, and

$$\mathrm{cov}\left({\mathbf{f}}_{*}\right)={\mathbf{K}}_{NN}-{\mathbf{K}}_{Nn}{\left[{\mathbf{K}}_{nn}+{\sigma }^2\mathbf{I}\right]}^{-1}{\mathbf{K}}_{Nn}^{\mathrm{T}},$$

(7)

denotes the covariance of ${\mathbf{f}}_{*}$.

The hyperparameter $\mathbf{\theta }$ incorporated in the mean and covariance functions underpin the predictive performance of GP models, and they are usually estimated by maximising the logarithm of the marginal likelihood

$$logp\left(\mathbf{y}\right|\mathbf{X})=-\frac{1}{2}{\mathbf{y}}^{\mathrm{T}}{\left({\mathbf{K}}_{nn}+{\sigma }^2\mathbf{I}\right)}^{-1}\mathbf{y}-\frac{1}{2}log|{\mathbf{K}}_{nn}+{\sigma }^2\mathbf{I}|-\frac{n}{2}log2\pi .$$

(8)

2.2. Neumann Series Approximation

Given a matrix inverse ${\mathbf{A}}^{-1}$, it can be expanded as the following Neumann series [21]

$${\mathbf{A}}^{-1}=\sum _{n=0}^{\infty }{\left({\mathbf{X}}^{-1}(\mathbf{X}-\mathbf{A})\right)}^n{\mathbf{X}}^{-1},$$

(9)

which holds if ${lim}_{n\to \infty }{\left(\mathbf{I}-{\mathbf{X}}^{-1}\mathbf{A}\right)}^n=\mathbf{0}$ is satisfied. In our case, suppose

$$\mathbf{A}=\mathbf{K}+{\sigma }_n^2\mathbf{I}\triangleq {\mathbf{D}}_A+{\mathbf{E}}_A,$$

(10)

where ${\mathbf{D}}_A$ is the main diagonal of $\mathbf{A}$ and ${\mathbf{E}}_A$ is the hollow. If we substitute $\mathbf{X}$ in Equation (9) by ${\mathbf{D}}_A$, we get

$${\mathbf{A}}^{-1}=\sum _{n=0}^{\infty }{\left(-{\mathbf{D}}^{-1}{\mathbf{E}}_A\right)}^n{\mathbf{D}}_A^{-1},$$

(11)

which is guaranteed to converge when ${lim}_{n\to \infty }{\left(-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A\right)}^n=\mathbf{0}$. We investigated the convergence condition in [17], where we proved that if $\mathbf{A}$ is diagonally dominant, then Neumann series can approximate ${\mathbf{A}}^{-1}$ both fast and accurate. In case $\mathbf{A}$ is not diagonally dominant, we also provided a way to convert it into a diagonally dominant matrix in [17], such that ${\mathbf{A}}^{-1}$ can still be approximated by Neumann series. When Neumann series given in (11) converges, we can then approximate $\mathbf{A}$ with only the first L terms. The L-term approximation is computed as follows:

$${\tilde{\mathbf{A}}}_L^{-1}=\sum _{n=0}^{L-1}{\left(-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A\right)}^n{\mathbf{D}}_A^{-1},$$

(12)

For instance, when $L=1,2,3$, we have the approximations

$${\tilde{\mathbf{A}}}_L^{-1}=\left\{\begin{array}{cc}{\mathbf{D}}_A^{-1},\hfill & \hfill L=1\\ {\mathbf{D}}_A^{-1}-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1},\hfill & \hfill L=2\\ {\mathbf{D}}_A^{-1}-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}+{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}.\hfill & \hfill L=3\end{array}\right.$$

(13)

3. Uncertainty Quantification in Gaussian Processes

3.1. Uncertainty in Measurements

It is intuitive that noisy measurements would result in less accurate predictions, just as a poor model would do. However, it is not direct from Equations (6) and (7). We will show in detail how the measurement noise would affect the prediction accuracy.

From Equations (6) and (7), we can see that the measurement noise $ϵ$ affects the prediction and the covariance by adding a term ${\sigma }_n^2\mathbf{I}$ to the prior covariance $\mathbf{K}$ in comparison to the noisy free scenario [20]. From the way that they originated, we know that both $\mathbf{K}$ and ${\sigma }_n^2\mathbf{I}$ are symmetrical. Then, a matrix $\mathbf{P}$ exists such that

$$\mathbf{K}={\mathbf{P}}^{-1}{\mathbf{D}}_{\mathbf{K}}\mathbf{P},$$

(14)

where ${\mathbf{D}}_{\mathbf{K}}$ is a diagonal matrix with eigen values of $\mathbf{K}$ along the diagonal. As ${\sigma }_n^2\mathbf{I}$ a diagonal matrix itself, we have

$${\sigma }_n^2\mathbf{I}={\mathbf{P}}^{-1}{\sigma }_n^2\mathbf{I}\mathbf{P}.$$

(15)

Therefore, we have the partial derivative of Equation (6) with respect to ${\sigma }_n^2$ as

$$\frac{\partial {\overline{\mathbf{f}}}_{*}}{\partial {\sigma }_n^2}={\mathbf{K}}_{*}\mathbf{P}{({\mathbf{D}}_{\mathbf{K}}+{\sigma }_n^2\mathbf{I})}^{-2}{\mathbf{P}}^{-1}\mathbf{y},$$

(16)

The element-wise form of Equation (16) can be therefore obtained as

$${\left(\frac{\partial {\overline{\mathbf{f}}}_{*}}{\partial {\sigma }_n^2}\right)}_o=-\sum _{h=1}^n\sum _{i=1}^n\sum _{j=1}^n{p}_{hj}{p}_{ij}{k}_{oh}{\Lambda }_j^{-1}{y}_i,$$

(17)

where ${\Lambda }_j={({\lambda }_j+{\sigma }_n^2)}^2$. $p_{hj}$ and $p_{ij}$ are the entries indexed by the j-th column, h-th and i-th row, respectively. $k_{oh}$ is the o-th row and h-th column entry of ${\mathbf{K}}_{*}$. $y_i$ is the i-th element of $\mathbf{y}$. $o=1,\cdots ,s$ denotes the o-th element of the partial derivation.

We can see that the sign of Equation (17) is determined by $p_{hj}$ and $p_{ij}$. This is because we can actually transform $\mathbf{y}$ to either positive or negative with a linear transformation, which will not be an issue for the GPs model. When we impose no constraints on $p_{hj}$ and $p_{ij}$, Equation (17) could be any real number, indicating that ${\overline{\mathbf{f}}}_{*}$ is multimodal with respect to ${\sigma }_n^2$, which means that one ${\sigma }_n^2$ can lead to different ${\overline{\mathbf{f}}}_{*}$, or equivalently, different ${\sigma }_n^2$ can lead to the same ${\overline{\mathbf{f}}}_{*}$. In such cases, it is difficult to investigate how ${\sigma }_n^2$ affects the prediction accuracy. In this paper, to facilitate the study of the monotonicity of ${\overline{\mathbf{f}}}_{*}$, we constrain $p_{hj}$ and $p_{ij}$ to satisfy

$${\left(\frac{\partial {\overline{\mathbf{f}}}_{*}}{\partial {\sigma }_n^2}\right)}_o\left\{\begin{array}{cc}>0,\hfill & \hfill p_{hj}{p}_{ij}<0,\\ <0,\hfill & \hfill p_{hj}{p}_{ij}>0,\\ =0,\hfill & \hfill p_{hj}{p}_{ij}=0.\end{array}\right.$$

(18)

Then, we can see that ${\overline{\mathbf{f}}}_{*}$ is monotonic. It means that changes of ${\sigma }_n^2$ can cause arbitrarily large/small predictions, whereas a robust method should bound the prediction errors regardless of how ${\sigma }_n^2$ varies.

Similarly, the partial derivative of Equation (7) with respect to ${\sigma }_n^2$ is

$$\frac{\partial \mathrm{cov}\left({\mathbf{f}}_{*}\right)}{\partial {\sigma }_n^2}=\left({\mathbf{K}}_{*}\mathbf{P}\right){({\mathbf{D}}_{\mathbf{K}}+{\sigma }_n^2\mathbf{I})}^{-2}{\left({\mathbf{K}}_{*}\mathbf{P}\right)}^T=\sum _{i=1}^n{\Lambda }_i^{-1}{\overrightarrow{\mathbf{p}}}_i{\overrightarrow{\mathbf{p}}}_i^T,$$

(19)

where we denote the $m\times n$ dimension matrix ${\mathbf{K}}_{*}\mathbf{P}$ as

$${\mathbf{K}}_{*}\mathbf{P}=[{\overrightarrow{\mathbf{p}}}_1,{\overrightarrow{\mathbf{p}}}_2,\cdots ,{\overrightarrow{\mathbf{p}}}_n],$$

(20)

with ${\overrightarrow{\mathbf{p}}}_i$ a $m\times 1$ vector, and $i=1,\cdots ,n$.

As the uncertainty is indicated by the diagonal elements, we only show how these elements change with respect to ${\sigma }_n^2$. The diagonal elements are given as

$$\begin{array}{cc}\hfill \mathrm{diag}\left(\sum _{i=1}^n{\Lambda }_i^{-1}{\overrightarrow{\mathbf{p}}}_i{\overrightarrow{\mathbf{p}}}_i^T\right)=& \mathrm{diag}\left(\sum _{i=1}^n{\Lambda }_i^{-1}{p}_{1i}^2,\sum _{i=1}^n{\Lambda }_i^{-1}{p}_{2i}^2,\cdots ,\sum _{i=1}^n{\Lambda }_i^{-1}{p}_{mi}^2\right)\hfill \\ \hfill =& \mathrm{diag}\left({\Sigma }_{11},{\Sigma }_{22},\cdots ,{\Sigma }_{mm}\right),\hfill \end{array}$$

(21)

with $\mathrm{diag}(·)$ denoting the diagonal elements of a matrix. We see that ${\Sigma }_{jj}⩾0$ stands for $j=1,\cdots ,m$, which implies that $\mathrm{cov}\left({\mathbf{f}}_{*}\right)$ is non-decreasing as ${\sigma }_n^2$ increases. This means that the increase of measurement noise level would cause the non-deceasing of the prediction uncertainty.

3.2. Uncertainty in Hyperparameters

Another factor that affects the prediction of a GPs model is the hyperparameters. In Gaussian processes, the posterior, as shown in Equation (5), is used to do the prediction, while the marginal likelihood is used for hyperparameters selection [18]. The log marginal likelihood as shown in Equation (22) is usually optimised to determine the hyperparameter with a specified kernel function.

$$logp\left(\mathbf{y}\right|\mathbf{X},\mathbf{\theta })=-\frac{1}{2}{\mathbf{y}}^T{(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}\mathbf{y}-\frac{1}{2}log|\mathbf{K}+{\sigma }_n^2\mathbf{I}|-\frac{N}{2}log2\pi .$$

(22)

However, the log marginal likelihood could be non-convex with respect to the hyperparameters, which implies that the optimisation may not converge to the global maxima [22]. A common solution dealing with it is to sample multiple starting points from a prior distribution, then choose the best set of hyperparameters according to the optima of the log marginal likelihood. Let’s assume $\mathbf{\theta }=\{{\theta }_1,{\theta }_2,\cdots ,{\theta }_s\}$ being the hyperparameter set and ${\theta }_s$ denoting the s-th of them, then the derivative of $logp\left(\mathbf{y}\right|\mathbf{X})$ with respect to ${\theta }_s$ is

$$\frac{\partial }{\partial {\theta }_s}logp\left(\mathbf{y}\right|\mathbf{X},\mathbf{\theta })=\frac{1}{2}\mathrm{tr}\left(\left(\mathbf{\alpha }{\mathbf{\alpha }}^T-{(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}\right)\frac{\partial (\mathbf{K}+{\sigma }_n^2\mathbf{I})}{\partial {\theta }_s}\right),$$

(23)

where $\mathbf{\alpha }={(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}\mathbf{y}$, and $\mathrm{tr}(·)$ denotes the trace of a matrix. The derivative in Equation (23) is often multimodal and that is why a fare few initialisations are used when conducting convex optimisation. Chen et al. show that the optimisation process with various initialisations can result in different hyperparameters [22]. Nevertheless, the performance (prediction accuracy) with regard to the standardised root mean square error does not change much. However, the authors do not show how the variation of hyperparameters affects the prediction uncertainty [22].

An intuitive explanation to the fact of different hyperparameters resulting with similar predictions is that the prediction shown in Equation (6) is non-monotonic itself with respect to hyperparameters. To demonstrate this, a direct way is to see how the derivative of (6) with respect to any hyperparameter ${\theta }_s\in \mathbf{\theta }$ changes, and ultimately how it affects the prediction accuracy and uncertainty. The derivatives of ${\overline{\mathbf{f}}}_{*}$ and $\mathrm{cov}\left({\mathbf{f}}_{*}\right)$ of ${\theta }_s$ are as below

$$\frac{\partial {\overline{\mathbf{f}}}_{*}}{\partial {\theta }_s}=\left({\mathbf{K}}_{*}\frac{\partial {(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}}{\partial {\theta }_s}+\frac{\partial {\mathbf{K}}_{*}}{\partial {\theta }_s}{(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}\right)\mathbf{y}.$$

(24)

We can see that Equations (24) and (25) are both involved with calculating ${(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}$, which becomes enormously complex when the dimension increases. In this paper, we focus on investigating how hyperparameters affect the predictive accuracy and uncertainty in general. Therefore, we use the Neumann series to approximate the inverse [21].

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{cov}\left({\mathbf{f}}_{*}\right)}{\partial {\theta }_s}=& \frac{\partial \mathbf{K}({\mathbf{X}}_{*},{\mathbf{X}}_{*})}{\partial {\theta }_s}-\frac{\partial {\mathbf{K}}_{*}}{\partial {\theta }_s}{(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}{\mathbf{K}}_{*}^T-{\mathbf{K}}_{*}\frac{\partial {(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}}{\partial {\theta }_s}{\mathbf{K}}_{*}^T\hfill \\ \hfill & -{\mathbf{K}}_{*}{(\mathbf{K}+{\sigma }_n^2\mathbf{I})}^{-1}\frac{\partial {\mathbf{K}}_{*}^T}{\partial {\theta }_s}.\hfill \end{array}$$

(25)

3.3. Derivatives Approximation with Neumann Series

The approximation accuracy and computationally complexity of Neumann series varies with L. This has been studied in [21,23], as well as in our previous work [17]. This paper aims at providing a way to quantify uncertainties involved in GPs. We therefore choose the 2-term approximation as an example to carry out the derivations. By substituting the 2-term approximation into Equations (24) and (25), we have

$$\begin{array}{c}\hfill \frac{\partial {\overline{\mathbf{f}}}_{*}}{\partial {\theta }_s}\approx \left({\mathbf{K}}_{*}\frac{\partial \left({\mathbf{D}}_A^{-1}-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}\right)}{\partial {\theta }_s}+\frac{\partial {\mathbf{K}}_{*}}{\partial {\theta }_s}\left({\mathbf{D}}_A^{-1}-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}\right)\right)\mathbf{y},\end{array}$$

(26)

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{cov}\left({\mathbf{f}}_{*}\right)}{\partial {\theta }_s}\approx & \frac{\partial \mathbf{K}({\mathbf{X}}_{*},{\mathbf{X}}_{*})}{\partial {\theta }_s}-\frac{\partial {\mathbf{K}}_{*}}{\partial {\theta }_s}\left({\mathbf{D}}_A^{-1}-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}\right){\mathbf{K}}_{*}^T\hfill \\ \hfill & -{\mathbf{K}}_{*}\frac{\partial \left({\mathbf{D}}_A^{-1}-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}\right)}{\partial {\theta }_s}{\mathbf{K}}_{*}^T-{\mathbf{K}}_{*}\left({\mathbf{D}}_A^{-1}-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}\right)\frac{\partial {\mathbf{K}}_{*}^T}{\partial {\theta }_s}.\hfill \end{array}$$

(27)

Due to the simple structure of matrices ${\mathbf{D}}_A$ and ${\mathbf{E}}_A$, we can get the element-wise form of Equation (26) as

$${\left(\frac{\partial {\overline{\mathbf{f}}}_{*}}{\partial {\theta }_s}\right)}_o=\sum _{i=1}^n\sum _{j=1}^n\left(k_{oj}\frac{\partial d_{ji}}{\partial {\theta }_s}+\frac{\partial k_{oj}}{\partial {\theta }_s}{d}_{ji}\right)y_i.$$

(28)

Similarly, the element-wise form of Equation (27) is

$${\left(\frac{\partial \mathrm{cov}\left({\mathbf{f}}_{*}\right)}{\partial {\theta }_s}\right)}_{oo}=\frac{\partial \mathbf{K}{({\mathbf{X}}_{*},{\mathbf{X}}_{*})}_{oo}}{\partial {\theta }_s}-\sum _{i=1}^n\sum _{j=1}^n\left(\frac{\partial k_{oj}}{\partial {\theta }_s}{d}_{ji}{k}_{oi}+k_{oj}\frac{\partial d_{ji}}{\partial {\theta }_s}{k}_{oi}-k_{oj}{d}_{ji}\frac{\partial k_{oi}}{\partial {\theta }_s}\right),$$

(29)

where $o=1,\cdots ,m$ denotes the o-th output, $d_{ji}$ is the j-th row and i-th column entry of ${\mathbf{D}}_A^{-1}-{\mathbf{D}}_A^{-1}{\mathbf{E}}_A{\mathbf{D}}_A^{-1}$, $k_{oj}$ and $k_{oi}$ are the o-th row, j-th and i-th entries of matrix ${\mathbf{K}}_{*}$, respectively. When the kernel function is determined, Equations (26)–(29) can be used for GPs uncertainty quantification.

3.4. Impacts of Noise Level and Hyperparameters on ELBO and UBML

The minimisation of $\mathrm{KL}\left(q(\mathbf{f},\mathbf{u})\parallel p(\mathbf{f},\mathbf{u}|\mathbf{y})\right)$ is equivalent to maximise the ELBO [18,24] as shown in

$${\mathcal{L}}_{\mathrm{lower}}=-\frac{1}{2}{\mathbf{y}}^T{\mathbf{G}}_n^{-1}\mathbf{y}-\frac{1}{2}log\left|{\mathbf{G}}_n\right|-\frac{N}{2}log\left(2\pi \right)-\frac{t}{2{\sigma }_n^2},$$

(30)

where ${\mathbf{G}}_n={\mathbf{G}}_{\mathbf{xx}}+{\sigma }_n^2\mathbf{I}$, and $t=\mathrm{Tr}({\mathbf{K}}_{\mathbf{xx}}-{\mathbf{G}}_{\mathbf{xx}})$. Combining it with UBML, as shown in Equation (31), an interval can be given to quantify the uncertainty in marginal likelihood.

$${\mathcal{L}}_{\mathrm{upper}}=\frac{1}{2}{\mathbf{y}}^T{\left({\mathbf{G}}_n+t\mathbf{I}\right)}^{-1}\mathbf{y}-\frac{1}{2}log\left|{\mathbf{G}}_n\right|-\frac{N}{2}log\left(2\pi \right).$$

(31)

This paper, however, focuses on investigating how ELBO and UBML change according to ${\sigma }_n^2$ only. Because the investigation of how ELBO and UBML change with respect to kernel hyperparameters involves multiple Neumann series approximations, which makes the analysis less convincing. We shall leave it as an open problem for future study. The derivatives of Equations (30) and (31) with respect to ${\sigma }_n^2$ are as follows,

$$\begin{array}{c}\hfill \frac{\partial {\mathcal{L}}_{\mathrm{lower}}}{\partial {\sigma }_n^2}=\frac{1}{2}\left[\sum _{i=1}^n{({\lambda }_i+{\sigma }_n^2)}^{-2}{\left(\sum _{j=1}^n{y}_j{v}_{ji}\right)}^2-\sum _{i=1}^n\frac{1}{{\lambda }_i+{\sigma }_n^2}+\frac{t}{{\sigma }_n^4}\right],\end{array}$$

(32)

$$\begin{array}{c}\hfill \frac{\partial {\mathcal{L}}_{\mathrm{upper}}}{\partial {\sigma }_n^2}=-\frac{1}{2}\left[\sum _{i=1}^n{({\lambda }_i+{\sigma }_n^2+t)}^{-2}{\left(\sum _{j=1}^n{y}_j{v}_{ji}\right)}^2+\sum _{i=1}^n\frac{1}{{\lambda }_i+{\sigma }_n^2}\right].\end{array}$$

(33)

Figure 1 shows how ${\sigma }_n^2$ affects ELBO and UBML. We set ${\sigma }_n^2$ to increase from 0.1 to 200.0 with a step of 0.01. Both ELBO and UBML are recorded step by step. From the figure, we can see that when ${\sigma }_n^2$ is small (${\sigma }_n^2\in [0.1,1.5]$), ELBO increases with different speeds, however, UBML fluctuates as the derivative of UBML jumps between positive and negative. When ${\sigma }_n^2$ is in [1.5, 3.0], ELBO still increases, but the speeds slow down significantly. In comparison, UBML keeps decreasing with reducing speeds. The decrements of UBML mean that when ${\sigma }_n^2$ increases, though ELBO could be increased still, but the maximum (which is the UBML) can decrease. When ${\sigma }_n^2\in [3.0,20.0]$, ELBO starts to decrease when ${\sigma }_n^2\approx 3.2$, while UBML keeps decreasing. This means that as ${\sigma }_n^2$ increases, both ELBO and UBML decrease, which indicates that the model becomes less and less effective to explain the data. When ${\sigma }_n^2$ keeps increasing (${\sigma }_n^2\in [20.0,200.0]$), the decreasing speeds of ELBO and UBML becomes similar and approaches zero. This means that UBML and ELBO both converge and together define an interval for the marginal likelihood, which however, can result in non-optimal hyperparameters. Our conclusion is that when ${\sigma }_n^2$ increases, UBML tends to decrease, which decreases the maximum that ELBO can reach. ELBO, on the other hand, is robust to the change of ${\sigma }_n^2$ (as it keeps increasing when ${\sigma }_n^2$ is below ∼3.2). However, when ${\sigma }_n^2$ exceeds a certain threshold, ELBO turns to decrease, indicating that the GPs model becomes less and less reliable. However, both ELBO and UBML converge, even when ${\sigma }_n^2$ becomes very significant, though we can no longer trust the model.

4. Experiments and Analysis

To verify that the proposed solution can help to identify the impacts of ${\sigma }_n^2$ and $\mathbf{\theta }$ on the predition accuracy and uncertainty of GPs model and its sparse variants such as the fully independent training conditional (FITC) [25] and variational free energy (VFE) [24] models, we conduct various experiments to process air quality data collected from Sheffield, UK, and Pershawar, Pakistan (see Appendix A), during the time period of 24 June 2019–14 July 2019 for three weeks, which will be denoted as W1, W2, and W3 hereafter. The data were collected with digital sensors called AQMesh pod with a 15 min time interval. Though the sensor itself is able to measure the concentrations of quite a few atmospheric pollutants, here we only analyse the concentrations of NO, NO${}_2$, SO${}_2$, and PM${}_{2.5}$. Figure 2 shows the raw data. We can see directly that the air quality of Sheffield is much better than Pershawar on average. Especially during daytime, concentrations of NO${}_2$ and PM${}_{2.5}$ in Pershawar exceed the WHO criteria (see Appendix B). Meanwhile, those in Sheffield are much lower than the criteria. Being a postindustrial city itself, Sheffield has improved air quality significantly. The experience can be spread to help cities like Pershwar to improve air quality.

4.1. Air Quality Prediction

Figure 3 and Figure 4 show Sheffield and Pershawar forecasting results of GPs, FITC, and VFE, with $3\sigma$ confidence intervals (denoted as Conf in the figures) indicated by the shaded area. We can see that the GPs model reports the best results in general, in terms of absolute error between predicts and measurements (denoted as Meas in the figures). However, the performance of all the models varies from pollutant types to cities. This is actually one of the reasons why the investigation of how measurement noise level and hyperparameters affect prediction accuracy and uncertainty is necessary. To make the results more convincing, we normalise the data from both cities for uncertainty quantification studies.

4.2. Impacts of Measurement Noise Level and Hyperparameters

To demonstrate how noise level ${\sigma }_n^2$ and hyperparameters affect prediction accuracy and uncertainty, three sets of experiments are conducted. This paper adopts the squared exponential (SE) kernel, with hyperparameters $s_f$ and l. The analytical derivation can be found in Appendix C. The prediction accuracy is identified by the root mean square error (RMSE), as shown in Equation (34), while the uncertainty is identified by $\frac{1}{2}\sigma$ confidence bound. Configurations of the experiments are as follows.

Experiment 1: Impacts of ${\sigma }_n^2$ on prediction accuracy and uncertainty. Both $s_f$ and l are fixed to be the optimised values. ${\sigma }_n^2$ varies from 0.1 through to 20.0. NO, ${\mathrm{NO}}_2$, ${\mathrm{SO}}_2$, and ${\mathrm{PM}}_{2.5}$ data from both cities are processed. Six inducing points are applied to both FITC and VFE.

Experiment 2: Impacts of $s_f$ on prediction accuracy and uncertainty. l is set to the optimised value. $s_f$ varies from 0.1 through to 30.0. ${\sigma }_n^2$ is set to 0.5 and 1.5, respectively. NO data from both cities are processed. Six inducing points are applied to both FITC and VFE.

Experiment 3: Impacts of l on prediction accuracy and uncertainty. $s_f$ is set to the optimised value. l varies from 0.1 through to 30.0. ${\sigma }_n^2$ is set to 0.5 and 1.5, respectively. NO data from both cities are processed. Six inducing points are applied to both FITC and VFE.

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{i=1}^{Num}{(y_i-\widehat{y_i})}^2}{Num}},$$

(34)

where $y_i$ is the ground truth value and $\widehat{y_i}$ represents predicted meant. $Num$ is the sample number in testing set.

Figure 5 and Figure 6 show the results from Experiment 1. To make the results more distinguishable, the horizontal axes of the figures are set to $log\left({\sigma }_n^2\right)$. We can see from Figure 5 that when ${\sigma }_n^2$ is small, GPs perform the best in general, while the performance of FITC and VFE varies. We can also observe that as ${\sigma }_n^2$ keeps increasing, the RMSE becomes very significant for all methods/pollutants. Similar results can be observed from Figure 6 as well. Both comply with our theoretical conclusions, despite the fact that the Neumann series is used to approximate the matrix inverse. We also notice that ${\sigma }_n^2$ has a more significant impact on Sheffield data as RMSE increases ealier after $log\left({\sigma }_n^2\right)$ reaches zero. From Figure 6b,c, we also see that the uncertainty bounds of Sheffield data are greater after $log\left({\sigma }_n^2\right)$ reaches zero. We think the reason is that Sheffield data are generally less periodical than Pershawar data (see Figure 2), which influences the performance of the models.

4.3. Impacts of Noise Level on ELBO and UBML

Figure 7 shows the results from Experiment 2. According to our theoretical results, the impact of $s_f$ on the uncertainty should become greater as $s_f$ increases. This is verified by the results shown in Figure 7b,d. Our theoretical results also suggest that the variation of $s_f$ would not affect the prediction accuracy. We can see from Figure 7a,c that when $s_f$ is smaller, it does affect the prediction accuracy, but when it exceeds a certain value, the impacts become negligible. Considering the Neumann series approximation, we would say that the experimental results comply with the theoretical conclusion.

The results of Experiment 3 are shown in Figure 8. We can see that when l is smaller, both RMSE and the uncertainty bounds change rapidly. While after it exceeds certain values, both converge. This again complies with our theoretical conclusions and simulation results. We should also notice from Figure 7 and Figure 8 that the increment of $s_f$ tends to increase the uncertainty, whereas the increment of l tends to decrease the uncertainty. Taking both into consideration, an optimised uncertainty bound can be obtained.

We also conduct an experiment to demonstrate how the noise level ${\sigma }_n^2$ affects the ELBO and UBML. In our experiment, we set ${\sigma }_n^2$ to vary from 0.5 to 4.5. The results are shown in Figure 9. To make the results distinguishable, we set the vertical axes to $log(-ELBO/UBML)$. To make the logrithm work, we reverse the signs of both ELBO and UBML. This is the reason why ELBO is ‘greater’ than UBML in Figure 9. The full GPs model is trained by setting ${\sigma }_n^2$ to $\{1,7,13,19,25,31,37,43,49\}$ to obtain 9 sets of hyperparameters. For each set of them, we then set ${\sigma }_n^2$ to vary from 0.5 to 4.5. The darker the colour in Figure 9, the smaller ${\sigma }_n^2$ is for model training. We can see that generally, greater ${\sigma }_n^2$ can slow down the convergence speed of both ELBO and UBML, while training a model. When the model is trained, the increment of ${\sigma }_n^2$ can lower down UBML, which is the maximum that ELBO can reach. This implies that the increment of ${\sigma }_n^2$ can cause the failure of a sparse GPs model, as ELBO is deeply related to determine a sparse GPs model. Nevertheless, the experimental results again comply with our theoretical conclusions.

5. Conclusions

This paper proposes a general method to investigate how the performance variation of a Gaussian process model can be attributed to hyperparameters and measurement noises, etc. The method is demonstrated by applying it to process particulate matter (e.g., ${\mathrm{PM}}_{2.5}$) and gaseous pollutants (e.g., NO, ${\mathrm{NO}}_2$, and ${\mathrm{SO}}_2$) from both Sheffield, UK, and Peshawar, Pakistan. Experimental results show that the proposed method provides insights on how measurement noises and hyperparameters, etc. affect the prediction performance of a Gaussian process. The results align with the analytical derivations, which is enabled by adopting Neuman series to approximate matrix inversions in Gaussian process models. The theoretical findings and experimental results combined demonstrate that the proposed method can generate air quality forecasting results. In the meantime, it provides a way to link uncertainties in measurements and hyperparameters, etc. with the forecasting results. This will help with forecasting performance analysis when measurement noise level or model hyperparameters vary, making the method more general.