Gaussian Processes Over Graphs

TL;DR

This paper introduces Gaussian processes tailored for signals on graphs, leveraging graph Laplacian regularization to improve prediction accuracy, especially with limited or noisy data, by enforcing specific spectral properties.

Contribution

It develops a novel Gaussian process framework over graphs that incorporates spectral regularization, reducing predictive variance compared to traditional GPs.

Findings

GPG has strictly smaller predictive variance than standard GP.

GPG outperforms GP with small training datasets.

GPG demonstrates robustness under noisy training conditions.

Abstract

We propose Gaussian processes for signals over graphs (GPG) using the apriori knowledge that the target vectors lie over a graph. We incorporate this information using a graph- Laplacian based regularization which enforces the target vectors to have a specific profile in terms of graph Fourier transform coeffcients, for example lowpass or bandpass graph signals. We discuss how the regularization affects the mean and the variance in the prediction output. In particular, we prove that the predictive variance of the GPG is strictly smaller than the conventional Gaussian process (GP) for any non-trivial graph. We validate our concepts by application to various real-world graph signals. Our experiments show that the performance of the GPG is superior to GP for small training data sizes and under noisy training.

Figures

Tables

Table 1. TABLE I: Names of different proteins that represent the nodes of the graph considered in Section IV-D.

Number	1	2	3	4	5	6	7	8	9	10	11
Name	praf	pmek	plcg	PIP2	PIP3	p44/42	pjnk	pakts473	PKA	PKC	P38

Equations

$$\mathbf{L}=\mathbf{V}{\mathbf{J}}_L{\mathbf{V}}^{\top },$$

$$\mathbf{V}$$

$$\hat{\mathbf{x}}={\mathbf{V}}^{\top }\mathbf{x}.$$

$${\mathbf{y}}_g=\arg \underset{\mathbf{z}}{\min }\Vert \mathbf{y}-\mathbf{z}{\Vert }_2^2+\alpha \Vert {\mathbf{J}}_p\hat{\mathbf{z}}{\Vert }_2^2,\alpha \ge 0$$

$$\Vert {\mathbf{J}}_p\hat{\mathbf{z}}{\Vert }_2^2={\hat{\mathbf{z}}}^{\top }{\mathbf{J}}_p^2\hat{\mathbf{z}}={\mathbf{z}}^{\top }\mathbf{V}{\mathbf{J}}_p^2{\mathbf{V}}^{\top }\mathbf{z}={\mathbf{z}}^{\top }\mathbf{G}\mathbf{z},$$

$${\mathbf{y}}_g=\arg \underset{\mathbf{z}}{\min }\Vert \mathbf{y}-\mathbf{z}{\Vert }_2^2+\alpha {\mathbf{z}}^{\top }\mathbf{G}\mathbf{z},\alpha \ge 0.$$

$${\mathbf{y}}_g={\left({\mathbf{I}}_M+\alpha \mathbf{G}\right)}^{-1}\mathbf{y}.$$

$${\mathbf{t}}_n=\mathbf{y}({\mathbf{x}}_n,\mathbf{W})+{\mathbf{e}}_n,$$

$$p({\mathbf{e}}_n)=\mathcal{N}({\mathbf{e}}_n|0,{\beta }^{-1}{\mathbf{I}}_M).$$

$$\mathbf{y}({\mathbf{x}}_n,\mathbf{W})={\mathbf{W}}^{\top }\mathbf{\varphi }({\mathbf{x}}_n)\triangleq {\mathbf{y}}_n,$$

$$p(\tilde{\mathbf{w}})=\mathcal{N}(\tilde{\mathbf{w}}|0,{\gamma }^{-1}{\mathbf{I}}_{KM}).$$

$${\mathbf{y}}_{g,n}$$

$${\mathbf{W}}_g=\mathbf{WB}.$$

$${\tilde{\mathbf{w}}}_g\triangleq \text{mbox}vec({\mathbf{W}}_g)=\text{mbox}vec(\mathbf{WB})=(\mathbf{B}\otimes {\mathbf{I}}_K)\tilde{\mathbf{w}},$$

$$p({\tilde{\mathbf{w}}}_g)=\mathcal{N}(\widetilde{{\mathbf{w}}_g}|0,{\gamma }^{-1}({\mathbf{B}}^2\otimes {\mathbf{I}}_K)).$$

$${\mathbf{t}}_n={\mathbf{y}}_g({\mathbf{x}}_n,{\mathbf{W}}_g)+{\mathbf{e}}_n^{\prime },$$

$$p({\mathbf{t}}_n|\mathbf{X},{\mathbf{W}}_g,\beta )=\mathcal{N}({\mathbf{t}}_n|{\mathbf{W}}_g^{\top }\mathbf{\varphi }({\mathbf{x}}_n),{\beta }^{-1}{\mathbf{I}}_M).$$

$$\mathbf{\Phi }$$

$${\mathbf{Y}}_g$$

$$\mathbf{T}$$

$$\tilde{\mathbf{t}}={\tilde{\mathbf{y}}}_g+\tilde{\mathbf{e}}$$

$$p(\tilde{\mathbf{t}})=\mathcal{N}(\tilde{\mathbf{t}}|\tilde{\mathbf{y}},{\beta }^{-1}{\mathbf{I}}_{MN}).$$

$$\tilde{\mathbf{y}}$$

$$\mathbb{E}\{\tilde{\mathbf{y}}\}$$

$$\mathbb{E}\{\tilde{\mathbf{y}}{\tilde{\mathbf{y}}}^{\top }\}$$

$${\mathbf{C}}_{g,N}=({\mathbf{B}}^2\otimes \mathbf{K})+{\beta }^{-1}{\mathbf{I}}_{MN},$$

$$\displaystyle\mathbf{C}_{g,N+1}=\left[\begin{array}[]{cc}\mathbf{C}_{g,N}&\mathbf{D}\\ \mathbf{D}^{\top}&\mathbf{F}\end{array}\right],$$

$$p({\mathbf{t}}_{N+1}|{\mathbf{t}}_1,\cdots ,{\mathbf{t}}_N)=\mathcal{N}({\mathbf{t}}_{N+1}|{\mathbf{\mu }}_{N+1},{\mathbf{\Sigma }}_{N+1}),$$

$${\mathbf{C}}_{c,N}=({\mathbf{I}}_M\otimes \mathbf{K})+{\beta }^{-1}{\mathbf{I}}_{MN},$$

$$p({\mathbf{t}}_{N+1}|{\mathbf{t}}_1,\cdots ,{\mathbf{t}}_N)=\mathcal{N}({\mathbf{t}}_{N+1}|{\mathbf{\mu }}_{c,N+1},{\mathbf{\Sigma }}_{c,N+1}),$$

$$\mathrm{tr}\left({\mathbf{C}}_{c,N}\right)$$

$$\mathrm{tr}(\Delta {\mathbf{C}}_N)$$

Full text

I Introduction

Gaussian processes are a natural extension of the ubiquitous kernel regression to the Bayesian setting where the regression parameters are modelled as random variables with a Gaussian prior distribution [1]. Given the training observations, Gaussian processes generate posterior probabilities of the target or output for new inputs or observations, as a function of the training data and the input kernel function [2]. Gaussian process models and its variants have been applied in a number of diverse fields such as model predictive control and system analysis [3, 4, 5, 6, 7], latent variable models [8, 9, 10, 11], multi-task learning [10, 12, 13], image analysis and synthesis [14, 15, 16, 17], speech processing [18, 19, 20], and magnetic resonance imaging (MRI) [21, 22]. Gaussian processes have also been extended to a non-stationary regression setting [23, 24, 25] and for regression over complex-valued data[26]. Recently, Gaussian processes were shown to be useful in training and analysis of deep neural networks, and that a Gaussian process can be viewed as a neural network with a single infinite-dimensional layer of hidden units [27, 28]. The prediction performance of the Gaussian process depends on the availability of training data, progressively improving as the training data size is increased. In many applications, however, we are required to make predictions using a limited number of observations which may be further corrupted with observation noise. For such cases, providing additional structure helps improve the prediction performance of the GP significantly. In this article, we advocate the use of graph signal processing for improving prediction performance in the lack of sufficient and reliable training data.

We propose Gaussian processes by incorporating the apriori knowledge that the vector-valued target or output vectors lie over an underlying graph. This forms a natural Bayesian extension of the kernel regression for graph signals proposed recently in [29]. In particular, the target vectors are enforced to follow a pre-specified profile in terms of the graph Fourier coefficients, such being lowpass, bandpass, or high-pass. We show that this in turn translates to a specific structure on the prior distribution of the target vectors. We derive the predictive distribution for a general input given the training observations, and prove that graph signal structure leads to decrease in variance of the predictive distribution. Our hypothesis is that incorporating the graph structure would boost the prediction performance. We validate our hypothesis on various real-world datasets such as temperature measurements, flow-cytometry data, functional MRI data, and tracer diffusion experiment for air pollution studies. Though we consider GPG mainly for undirected graphs characterized by the graph-Laplacian in our analysis, we also discuss how our approach may be extended to handle directed graphs using an appropriate regularization.

II Preliminaries on Graph Signal Processing

Graph signal processing or signal processing over graphs deals with extension of several traditional signal processing methods while incorporating the graph structural information [30, 31]. This includes signal analysis concepts such as Fourier transforms[32], filtering [30, 33], wavelets [34, 35], filterbanks [36, 37], multiresolution analysis [38, 39, 40], denoising [41, 42], and dictionary learning [43, 44], and stationary signal analysis [45, 46]. Spectral clustering and principal component analysis approaches based on graph signal filtering have also been proposed recently [47, 48]. Several approaches have been proposed for learning the graph directly from data [49, 50]. Recently, kernel regression based approaches have also been developed for graph signals [29, 51, 52]. We now briefly review the relevant concepts from graph signal processing which will be used in our analysis and development. This review has been included to keep our article self-contained.

II-A Graph Laplacian and graph Fourier spectrum

Consider an undirected graph with $M$ nodes and adjacency matrix $\mathbf{A}$. The $(i,j)$th entry of $\mathbf{A}$ denotes the strength of the edge between the $i$th and $j$th nodes, $\mathbf{A}(i,j)=0$ denoting absence of edge. Since the graph is undirected, we have symmetric edge-weights or $\mathbf{A}=\mathbf{A}^{\top}$. The graph-Laplacian matrix $\mathbf{L}$ is defined as $\mathbf{L=D-A}$, where $\mathbf{D}$ is the diagonal degree matrix with $i$th diagonal element given by the sum of the elements in the $i$th row of $\mathbf{A}$[53]. $\mathbf{L}$ is symmetric for an undirected graph and by construction has nonnegative eigenvalues with the smallest eigenvalue being equal to zero. A graph signal $\mathbf{y}=[y(1)\,y(2)\,\cdots y(M)]^{\top}\in\mathbb{R}^{M}$ is an $M$-dimensional vector such that $y(i)$ denotes the value of the signal at the $i$th node, where $\top$ denotes transpose operation. The smoothness of $\mathbf{y}$ is measured using $l(\mathbf{y})=\displaystyle\frac{\mathbf{y}^{\top}\mathbf{L}\mathbf{y}}{\mathbf{y}^{\top}\mathbf{y}}=\frac{1}{\mathbf{y}^{\top}\mathbf{y}}{\sum_{\mathbf{A}(i,j)\neq 0}\mathbf{A}(i,j)(y(i)-y(j))^{2}}$. The quantity $l(\mathbf{y})$ is small when $\mathbf{y}$ takes the similar value across all connected nodes, in agreement with what one intuitively expects of a smooth signal. Similarly, $\mathbf{y}$ is a high-frequency or non-smooth signal if it has dissimilar values across connected nodes, or equivalently a large value of $l(\mathbf{y})$.

The eigenvectors of the graph-Laplacian are used to define the notion of frequency of graph signals. Let $\{\lambda_{i}\}_{i=1}^{M}$ and $\{\mathbf{v}_{i}\}_{i=1}^{M}\in\mathbb{R}^{M}$ denote the eigenvalues and eigenvectors of $\mathbf{L}$, respectively. Then, the eigenvalue decomposition of $\mathbf{L}$ is

[TABLE]

where $\mathbf{J}_{L}$ is the diagonal eigenvalue matrix of $\mathbf{L}$, such that

[TABLE]

It is standard practice to order the eigenvectors $\mathbf{v}_{i}$ of the graph-Laplacian according to their smoothness in terms of $l(\mathbf{v}_{i})$. The eigenvectors $\mathbf{v}_{i}$ are referred to as the graph Fourier transform (GFT) basis vectors since they generalize the notion of the discrete Fourier transform [31]. Let $\lambda_{i}$ denote the $i$th eigenvalue of $\mathbf{L}$ ordered as $0=\lambda_{1}\leq\lambda_{2}\cdots\leq\lambda_{M}$. Then, we observe that $l(\mathbf{v}_{i})=\lambda_{i}$. In other words, the eigenvectors corresponding to smaller $\lambda_{i}$ vary smoothly over the graph and those with large $\lambda_{i}$ exhibit more variation across the graph. This in turn gives an intuitive frequency ordering of the GFT basis vectors. The GFT coefficients of a graph signal $\mathbf{x}$ are defined as the inner product of the signal with the GFT basis vectors $\mathbf{V}$, that is, the GFT coefficient vector $\hat{\mathbf{x}}$ is defined as:

[TABLE]

The inverse GFT is then given by $\mathbf{x}=\mathbf{V}\hat{\mathbf{x}}$. A graph signal may then be low-pass or smooth, high-pass, band-pass, or band-stop according to the distribution of its GFT coefficients in $\hat{\mathbf{x}}$.

II-B Generative model for graph signals

We now derive the expressions for the graph signal of a particular GFT profile closest to any given signal, which shall be used in our analysis later. For any signal $\mathbf{y}$, a graph signal with a specified graph frequency profile $\mathbf{y}_{g}$ closest to $\mathbf{y}$ can be obtained by solving the following generative model

[TABLE]

where $\hat{\mathbf{z}}=\mathbf{V}^{\top}\mathbf{z}$ is the GFT of $\mathbf{z}$ and $\mathbf{J}_{p}=\mathrm{diag}(J(1),J(2),\cdots,J(M))$ is the diagonal matrix whose values penalize or promote specific GFT coefficients. For example, if $\mathbf{z}$ is to have energy predominantly in GFT vectors with indices $\{i,j,k,l\}$, we set the corresponding diagonal entries of $\mathbf{J}_{p}$ to zero and assign large values to the remaining diagonal entries. From the properties of the GFT, we have that

[TABLE]

where $\mathbf{G}=\mathbf{V}\mathbf{J}^{2}_{p}\mathbf{V}^{\top}$. Then, the generative model is equivalently expressible as:

[TABLE]

Since all the eigenvalues of $\mathbf{G}$ are nonnegative, $\mathbf{G}$ is positive semidefinite giving the unique closed-form solution for $\mathbf{y}_{g}$ as

[TABLE]

We note that most graph signals encountered in practice are usually smooth over the associated graph. One of the simplest cases of generating a smooth graph signal is to penalize the frequencies linearly by setting ${J}(i)=\sqrt{\lambda_{i}}$. This in turn corresponds to $\mathbf{G}=\mathbf{V}\mathbf{J}_{L}\mathbf{T}^{\top}=\mathbf{L}$ and the smooth graph signal then becomes $\mathbf{y}_{g}=\left(\mathbf{I}_{M}+\alpha\mathbf{L}\right)^{-1}\mathbf{y}$.

III Gaussian process over graph

We first develop Bayesian linear regression over graph and then proceed to develop a Gaussian process over graph.

III-A Bayesian linear regression over graph

We now propose linear regression over graphs. Let $\{\mathbf{x}_{n}\}_{n=1}^{N}$ denote the set of $N$ input observations, $n=1,\cdots,N$ and $\mathbf{t}_{n}\in\mathbb{R}^{M}$ denote the corresponding vector target values. We model the target as the output of a linear regression such that

[TABLE]

where $\mathbf{e}_{n}$ denotes the additive noise vector following a zero mean isotropic Gaussian distribution, that is,

[TABLE]

Here $\mathbf{I}_{M}$ denotes the identity matrix of dimension $M$ and $\beta$ is the precision parameter for the noise. The signal $\mathbf{y}$ is

[TABLE]

where $\boldsymbol{\phi}(\mathbf{x})\in\mathbb{R}^{K}$ is a function of $\mathbf{x}$, $\mathbf{W}\in\mathbb{R}^{K\times M}$ denotes the regression coefficient matrix. We further consider an isotropic Gaussian prior with precision $\gamma$ for the entries of $\mathbf{W}$. Let us use the notation $\tilde{\mathbf{w}}\triangleq\mbox{vec}(\mathbf{W})$ where $\mbox{vec}(\mathbf{W})$ denotes the vectorization of $\mathbf{W}$ obtained by concatenating the columns of $\mathbf{W}$ into a single vector. Then, we have that

[TABLE]

Our principal assumption is that the predicted target vector or regression output $\mathbf{y}_{n}$ has a graph Fourier spectrum as characterized by the diagonal spectrum matrix $\mathbf{J}_{p}$. However, $\mathbf{y}_{n}$ is not necessarily satisfy this requirement for an arbitrary choice of $\mathbf{W}$, since $\boldsymbol{\phi}(\cdot)$ is also fixed apriori and has not been assumed to be graph-specific. It becomes clear then that the prior on $\tilde{\mathbf{w}}$ should be chosen in a way that it promotes $\mathbf{y}_{n}$ to have the required graph Fourier spectrum. We next discuss our strategy for formulating such a prior distribution. Given the regression output $\mathbf{y}_{n}$ generated using a fixed $\mathbf{W}$ drawn from (6), the graph signal $\mathbf{y}_{g,n}$ closest to $\mathbf{y}_{n}$ is obtained by using the generative model (2) with the following result

[TABLE]

where $\mathbf{W}_{g}\triangleq\mathbf{W}\left(\mathbf{I}_{M}+\alpha\mathbf{G}\right)^{-1}$. Let $\mathbf{B}=\left(\mathbf{I}_{M}+\alpha\mathbf{G}\right)^{-1}$, then we have that

[TABLE]

Thus, we note that regression of $\boldsymbol{\phi}(\mathbf{x})$ with regression coefficient matrix $\mathbf{W}_{g}$ produces graph signals possessing the desired graph Fourier profile. In the case of $\mathbf{G}=\mathbf{L}$, this is the same as regression output being smooth over the graph. On vectorizing both sides of (7) and using properties of the Kronecker product [54], we get that

[TABLE]

where $\otimes$ denotes the Kronecker product operation. Note that $\mathbf{B}$ is symmetric as $\mathbf{L}$ and hence, $\mathbf{G}$ is symmetric. Since $p(\tilde{\mathbf{w}})=\mathcal{N}(\mbox{vec}(\tilde{\mathbf{w}})|\mathbf{0},\gamma^{-1}\mathbf{I}_{KM})$, we get that $\tilde{\mathbf{w}}_{g}$ is distributed according to

[TABLE]

This in turn implies that choosing a prior for the regression coefficients of the form (8) yields graph signals with specified graph Fourier spectrum over the graph with the Laplacian matrix $\mathbf{L}$. We then pose the regression problem (III-A) in the following form

[TABLE]

We assume that $\mathbf{e}^{\prime}_{n}$ follows the same distribution of $\mathbf{e}_{n}$ in (III-A). Then, the conditional distribution of $\mathbf{t}_{n}$ is given by

[TABLE]

We define the matrices $\mathbf{\Phi}$, $\mathbf{Y}_{g}$, and $\mathbf{T}$ as follows:

[TABLE]

On accumulating all the observations we have

[TABLE]

where $\tilde{\mathbf{t}}\triangleq\mathrm{vec}(\mathbf{T})$, $\tilde{\mathbf{y}}_{g}\triangleq\mathrm{vec}(\mathbf{Y}_{g})$ and $\tilde{\mathbf{e}}$ is the noise term. Dropping the dependency of the distribution on $\mathbf{X},\mathbf{W}_{g},\beta$ for brevity, we have that

[TABLE]

Since $\mathbf{Y}_{g}=\mathbf{\Phi}\mathbf{W}_{g}$, we have that

[TABLE]

The prior distribution of $\tilde{\mathbf{y}}$ is then a multivariate Gaussian distribution with mean and covariance given by

[TABLE]

We refer to this model as Bayesian linear regression over graph.

III-B Gaussian process over graph

We next develop Gaussian process over graphs (GPG). We first note the existence of $\mathbf{\Phi}\mathbf{\Phi}^{\top}$ in (12) and that the input $\mathbf{x}_{n}$ enters the equation only in the form of inner products of $\boldsymbol{\phi}(\cdot)$. The inner product $\boldsymbol{\phi}(\mathbf{x}_{n})^{\top}\boldsymbol{\phi}(\mathbf{x}_{m})$ is a measure of similarity between the $m$th and $n$th inputs. Keeping this in mind, one may generalize the inner-product $\boldsymbol{\phi}(\mathbf{x}_{m})^{\top}\boldsymbol{\phi}(\mathbf{x}_{n})$ to any valid kernel function[2] $k(\mathbf{x}_{m},\mathbf{x}_{n})$ of the inputs $\mathbf{x}_{m}$ and $\mathbf{x}_{n}$. The kernel matrix is denoted by $\mathbf{K}=\gamma^{-1}\mathbf{\Phi}\mathbf{\Phi}^{\top}$, and its $(m,n)$th entry is $k(\mathbf{x}_{m},\mathbf{x}_{n})$. Following (10), (11), and (III-A) we find that $\tilde{\mathbf{t}}$ is Gaussian distributed with zero mean and following covariance

[TABLE]

where we use the first subscript $g$ to denote the graph, and the second subscript $N$ to show explicit dependency on $N$ observations. With $(N+1)$ samples, we have that

[TABLE]

where

• $\mathbf{D}=\mathbf{B}^{2}\otimes\mathbf{k}$,
• $\mathbf{k}=[k(\mathbf{x}_{N+1},\mathbf{x}_{1}),k(\mathbf{x}_{N+1},\mathbf{x}_{2}),\cdots,k(\mathbf{x}_{N+1},\mathbf{x}_{N})]^{\top}$, and
• $\mathbf{F}=k(\mathbf{x}_{N+1},\mathbf{x}_{N+1})\mathbf{B}^{2}+\beta^{-1}\mathbf{I}_{M}$.

Then, using properties of conditional probability for jointly Gaussian vectors [2], we have the predictive distribution of GPG as follows

[TABLE]

where

• $\mathbf{\boldsymbol{\mu}}_{N+1}=\mathbf{D}^{\top}\mathbf{C}_{g,N}^{-1}\tilde{\mathbf{t}}$ and $\mathbf{\Sigma}_{N+1}=\mathbf{F}-\mathbf{D}^{\top}\mathbf{C}_{g,N}^{-1}\mathbf{D}$.

We note that in the case of completely disconnected graph, which corresponds to the conventional Gaussian process (GP), the graph adjacency matrix is equal to the identity matrix and correspondingly, $\mathbf{L=0}$. Then, the covariance of the joint distribution of $\tilde{\mathbf{t}}$ is given by

[TABLE]

where the subscript $c$ denotes it being the conventional Gaussian process. Correspondingly, the predictive distribution of $\mathbf{t}_{N+1}$ is given by

[TABLE]

where

• $\mathbf{\boldsymbol{\mu}}_{c,N+1}=\mathbf{D}_{c}^{\top}\mathbf{C}_{c,N}^{-1}\tilde{\mathbf{t}}$,
• $\mathbf{\Sigma}_{c,N+1}=\mathbf{F}_{c}-\mathbf{D}_{c}^{\top}\mathbf{C}_{c,N}^{-1}\mathbf{D}_{c}$,
• $\mathbf{D}_{c}=\mathbf{I}_{M}\otimes\mathbf{k}$,
• $\mathbf{F}_{c}=(k(\mathbf{x}_{N+1},\mathbf{x}_{N+1})+\beta^{-1})\mathbf{I}_{M}$.

We note that since Bayesian linear regression on graphs developed in Section A is a special case of the Gaussian process on graphs with kernels, we shall refer to the former as Gaussian process-Linear (GPG-L) and the latter as Gaussian process-Kernel (GPG-K). Correspondingly, we refer to the conventional versions as GP-L, and GP-K, respectively.

We next show that use of graph information reduces the variance of the predictive distribution. This implies that the GPG models the observed training samples better than conventional GP.

Theorem 1** (Reduction in variance).**

Graph structure reduces the variance of the marginal distribution $p(\tilde{\mathbf{t}})$, that is, GPG results in smaller variance of $p(\tilde{\mathbf{t}})$ than GP:

[TABLE]

Proof.

In order to prove the result, we need to show that the trace of $\Delta\mathbf{C}_{N}:=\mathbf{C}_{c,N}-\mathbf{C}_{g,N}$ is nonnegative. Using the properties of trace operation, we have that

[TABLE]

Since $\mathbf{K}$ is positive semidefinite by construction as a kernel matrix, we have that $\mathrm{tr}(\Delta\mathbf{C}_{N})\geq 0$ if and only if $\mathrm{tr}(\mathbf{I}_{M}-\mathbf{B}^{2})\geq 0$. Let $\{\theta_{i}\}_{i=1}^{N}$ and $\{\mathbf{u}_{i}\}_{i=1}^{N}\in\mathbb{R}^{N}$ denote the eigenvalues and eigenvectors of $\mathbf{K}$, respectively. Then, the eigenvalue decomposition of $\mathbf{G}$ and $\mathbf{K}$ are given by

[TABLE]

where $\mathbf{J}_{K}$ and $\mathbf{J}_{G}$ denote the diagonal eigenvalue matrix of $\mathbf{K}$ and $\mathbf{G}$, respectively, such that

[TABLE]

Since $\mathbf{J}_{G}=\mathbf{J}^{2}_{p}=\mathrm{diag}(J^{2}(1),J^{2}(2),\cdots J^{2}(M))$, we have

[TABLE]

Let $s$ denote the graph eigenvalue index corresonding to the denote the smallest non-zero diagonal value in $\mathbf{J}_{p}$, the value given by $J(s)$. Then, we have that $\displaystyle\frac{1}{(1+\alpha J^{2}(i))}\leq\frac{1}{(1+\alpha J^{2}(s))}$ $\forall J^{2}(i)\neq 0$. Then, a sufficient condition to ensure $\mathrm{tr}(\mathbf{I}_{M}-\mathbf{B}^{2})\geq 0$ is to impose that

[TABLE]

A graph with atleast one connected subgraph has atleast one eigenvalue of $\mathbf{L}$ strictly greater than 0 [53]. Thus, barring the completely disconnected graph and the pathological case of $\mathbf{J}_{p}=\mathbf{0}$ (or $\mathbf{G}=\mathbf{0}$), any graph with connections across nodes results in $\alpha J^{2}(s)>0$ or in other words that the variance strictly reduces in comparison with the conventional regression, that is,

[TABLE]

In other words, introducing non-zero connections among the nodes ensures that the marginal variance is strictly lesser than that obtained for the conventional regression case. We also note that in order for the reduction in variance to be significant, $\lambda_{s}$ must be large, which in turn implies that the connected communities within the graph must have strong algebraic connectivity [53]. In the case of regression output being smooth graph signals in the sense of $\mathbf{G}=\mathbf{L}$, we have that $J^{2}(i)=\lambda_{i}$ and $s=K$ for a graph with $K$-connected components or disjoint subgraphs. This is because the number of zero eigenvalues of $\mathbf{L}$ is equal to the nunber of connected components in the graph [53, 55]. ∎

An immediate consequence of the Theorem is the following important Corrollary which informs us that the variance of the predicted target reduces when the graph signal structure is employed:

Corollary 1.1** (Reduction in predictive variance).**

GPG-K with a non-trivial graph has strictly smaller predictive variance than GPG-K, that is,

[TABLE]

Proof.

We are required to prove that $\mathrm{tr}\left(\mathbf{\Sigma}_{c,N+1}-\mathbf{\Sigma}_{g,N+1}\right)>0$, or equivalently that $\Delta\boldsymbol{\Sigma}:=\mathbf{\Sigma}_{c,N+1}-\mathbf{\Sigma}_{g,N+1}$ is a positive semidefinite matrix. We notice that $\Delta\mathbf{C}_{N+1}$ is given by

[TABLE]

We observe that $\Delta\boldsymbol{\Sigma}_{N+1}$ is then the Schur complement of $\mathbf{F}_{c}-\mathbf{F}$ in $\Delta\mathbf{C}_{N+1}$. Since the Schur complement of a positive-definite matrix is also positive-definite [56], and we have already proved that $\Delta\mathbf{C}_{N+1}$ is positive-definite from Theorem 1, it follows that

[TABLE]

∎

By the preceding analysis, we also observe that the variance of the joint distribution $\mathbf{C}_{g,N+1}$, and hence, that of $\boldsymbol{\Sigma}_{g,N+1}$ is inversely related to the graph regularization parameter $\alpha$. This is because a large $\alpha$ implies a small value of $\displaystyle\frac{1}{1+\alpha J^{2}(i)}$ for all $i$, and therefore, a small value of $\mbox{tr}(\mathbf{C}_{g,N+1})$.

III-C On the mean vector of predictive distribution

We now show that the mean of predictive distribution is a graph signal with a graph Fourier spectrum which adheres to the condition imposed by $\mathbf{J}_{p}$. We demonstrate this by computing the graph Fourier spectrum of the predictive mean. Using the eigendecompostions of $\mathbf{G}$ and $\mathbf{K}$, we have that

[TABLE]

where

1. $\mathbf{J}=(\left(\mathbf{I}+\alpha\mathbf{J}_{G}\right)^{-2}\otimes\mathbf{J}_{K}+\beta^{-1}\mathbf{I}_{MN})^{-1}$,
2. $\eta_{i}$ is the $i$th diagonal element of $\mathbf{J}$,
3. $\mathbf{Z}=\mathbf{V}\otimes\mathbf{U}$,
4. $\mathbf{z}_{i}$ denotes the $i$th column vector of $\mathbf{Z}$, and
5. $\rho_{i}=\mathbf{z}_{i}^{\top}\tilde{\mathbf{t}}$.

Since the $\eta_{i}$ is a function of some $i1$th eigenvalue of $\mathbf{L}$ and $i2$th eigenvalue of $\mathbf{K}$, we shall alternatively use the notation $\eta_{i1,i2}$ to be more specific in the following analysis. Similarly, $\rho_{i}$ expressed as $\rho_{i1,i2}$ and the eigenvectors $\mathbf{z}_{i}=\mathbf{v}_{i1}\otimes\mathbf{u}_{i2}$. The component of the prediction mean along the graph eigenvector $\mathbf{v}_{k}$ (or the $k$th graph frequency) is then given by

[TABLE]

We observe that a nonzero $\alpha$ reduces or shrinks the contribution from the graph-frequencies corresponding to larger $J^{2}(k)$, in comparison with the conventional mean obtained by setting $\alpha=0$. In the case of smooth graph signals, we have $J^{2}(k)=\lambda_{k}$ implying that the prediction mean has lower contributions from higher graph-frequencies, which in turn shows that GPG performs a noise-smoothening by making use of the graph topology. However, this does not imply that the value of $\alpha$ may be set to be arbitrarily large. A large $\alpha$ will in turn force the resulting target predictions to lie close to the smoother eigenvectors of $\mathbf{L}$, which is not desirable since it reduces the learning ability of the GP. We note that since GPG-L is a special case of the GPG-K, Theorem 1 and the analysis following it apply directly also to GPG-L.

III-D Extension for directed graphs

In our analysis, we have assumed the underlying graph for the target vectors to be undirected and hence, characterized by the symmetric positive semidefinite graph-Laplacian matrix $\mathbf{L}$. We now discuss how GPG may be derived for the case of directed graphs, that is, when the adjacency matrix of the graph is assymetric. In the case of directed graphs, the following metric is popular for quantifying the smoothness of the graph signals:

[TABLE]

where $\mathbf{y}$ and $\mathbf{A}$ denote the graph signal and adjacency matrix, respectively, the adjacency matrix assumed to be normalized to have maximum eigenvalue modulus equal to unity [30]. A signal $\mathbf{y}$ is smooth over a directed graph if it has a small $\mbox{MS}_{g}$ value. This is because $\mbox{MS}_{g}$ measures the difference between the signal and its one-step diffused or ’graph-shifted’ version, thereby measuring the difference between the signal value at each node and its neighbours, weighted by the strength of the edge between the nodes. In the case of directed graphs, we may adopt the same approach as that of the undirected graph case with one important distinction: instead of solving (2), we now solve the following problem for each observation:

[TABLE]

Then, we have that

[TABLE]

Using a similar analysis as with the undirected graph case, we arrive at the following Gaussian process model:

[TABLE]

where $\tilde{\mathbf{y}}$ is distributed according to the prior:

[TABLE]

where $\mathbf{B}_{d}=\left(\mathbf{I}_{M}+\alpha(\mathbf{I-A})^{\top}(\mathbf{I-A})\right)^{-1}$. By following similar analysis as with the undirected graphs, it is possible to show that the variance of the predictive distribution reduces using GPG and that the predictive mean is smooth over the graph. In the interest of space and to avoid repetition, we do not include the analysis for directed graphs here.

IV Experiments

We consider application of GPG to various real-world signal examples. For the examples, we consider undirected graphs. Our interest is to compute the predictive distribution (17) given the noisy targets $\mathbf{T}$ and the corresponding inputs $\mathbf{X}$. Our assumption is that a target vector is smooth over an underlying graph. We use the graph-regularization with $\mathbf{G}=\mathbf{L}$. To evaluate the prediction performance, we use the normalized-mean-square-error (NMSE) defined as follows:

[TABLE]

where $\mathbf{Y}$ denotes the mean of the predictive distribution and $\mathbf{T}_{0}$ the true value of target matrix, that means $\mathbf{T}_{0}$ does not contain any noise. The noisy target matrix $\mathbf{T}$ is generated obtained by adding white Gaussian noise with precision parameter $\beta$ to $\mathbf{T}_{0}$. In the case of real-world examples, we compare the performance of the following cases:

GP with linear regression (GP-L): ${k}_{m,n}=\gamma^{-1}\boldsymbol{\phi}(\mathbf{x}_{m})^{\top}\boldsymbol{\phi}(\mathbf{x}_{n})$, where $\boldsymbol{\phi}(\mathbf{x})=\mathbf{x}$ and $\alpha=0$, 2. 2.

GPG with linear regression (GPG-L): ${k}_{m,n}~{}=~{}\gamma^{-1}\boldsymbol{\phi}(\mathbf{x}_{m})^{\top}\boldsymbol{\phi}(\mathbf{x}_{n})$ and $\alpha>0$, where $\boldsymbol{\phi}(\mathbf{x})=\mathbf{x}$, 3. 3.

GP with kernel regression (GP-K): Using radial basis function (RBF) kernel ${k}_{m,n}~{}=~{}\gamma^{-1}\displaystyle\exp\left(-\frac{\|\mathbf{x}_{m}-\mathbf{x}_{n}\|_{2}^{2}}{\sigma^{2}}\right)$ and $\alpha=0$, and 4. 4.

GPG with kernel regression over graphs (GPG-K): Using RBF kernel ${k}_{m,n}~{}=~{}\gamma^{-1}\displaystyle\exp\left(-\frac{\|\mathbf{x}_{m}-\mathbf{x}_{n}\|_{2}^{2}}{\sigma^{2}}\right)$ and $\alpha>0$.

We use five-fold cross-validation to obtain the values of regularization parameters $\alpha$ and $\gamma$. We set the kernel parameter $\sigma^{2}=\displaystyle\sum_{m,n}\|\mathbf{x}_{m}-\mathbf{x}_{n}\|_{2}^{2}$. We perform experiments under various signal-to-noise ratio (SNR) levels, and choose the precision parameter $\beta$ accordingly.

IV-A Prediction for fMRI in cerebellum graph

We first consider the functional magnetic resonance imaging (fMRI) data obtained for the cerebellum region of brain [57]111The data is available publicly at https://openfmri.org/dataset/ds000102.. The original graph consists of 4465 nodes corresponding to different voxels of the cerebellum region. The voxels are mapped anatomically following the atlas template [58, 59]. We refer to [57] for details of graph construction and associated signal extraction. We consider a subset of the first 100 vertices in our analysis. Our goal is to use the first ten vertices as input $\mathbf{x}\in\mathbb{R}^{10}$ to make predictions for remaining 90 vertices, which forms the output $\mathbf{t}\in\mathbb{R}^{90}$. Thus, the target signals lie over a graph of dimension $M=90$. The corresponding adjacency matrix is shown in Figure 1(a). We have a total of 295 graph signals corresponding to different measurements from a single subject. We use a portion of the signals for training and the remaining for testing. We construct noisy training targets by adding white Gaussian noise at SNR-levels of 10 dB and 0 dB. The NMSE of the prediction mean for testing data, averaged over 100 different random choices of training and testing sets is shown in 1(b) and (c); this is Monte Carlo simulation to check robustness. We observe that for both linear and kernel regression cases, GPG outperforms GP by a significant margin, particularly at small training data sizes as expected. The trend is also similar when larger subsets of nodes from the full set are considered. The results are not reported here for brevity and to avoid repetition.

IV-B Prediction for temperature data

We next apply GPG on temperature measurements from the 45 most populated cities in Sweden for a period of three months from October to December 2017. The data is available publicly from the Swedish Meteorological and Hydrological Institute [60]. Our goal is to perform one-day temperature prediction: given the temperature data for a particular day $\mathbf{x}_{n}$, we predict the temperature for the next day $\mathbf{t}_{n}$. We have 90 input-target data pairs in total, of which one half is used for training and the rest for testing. We consider the geodesic graph in our analysis with the $(i,j)$th entry of the graph adjacency matrix $\mathbf{A}$ is given by

[TABLE]

where $d_{ij}$ denotes the geodesic distance between the $i$th and $j$th cities. In order to remove self loops, the diagonal of $\mathbf{A}$ is set to zero. We generate noisy training data by adding zero-mean white Gaussian noise at SNR of $5$ dB and [math] dB to the true temperature measurements. In Figure 2, we show the NMSE obtained for testing data by averaging over 100 different random partitioning of the total dataset into training and testing datasets. We observe that GPG outperforms GP for both linear and kernel regression cases.

IV-C Prediction for flow-cytometry data

We now consider the application of GPG to flow-cytometry data considered by Sachs et al. [61] which consists of response or signalling level of 11 proteins in different experiment cells. Since the protein signal values have a large dynamic range and are all positive values, we perform experiments on signals obtained by taking a logarithm with the base of 10 for reducing the dynamic range. We use the first 1000 measurements in our analysis. We use the symmetricized version of the directed unweighted acyclic graph proposed by Sachs et al. [61]. Among the 11 proteins, we choose proteins 10 and 11 arbitrarily as input to make predictions for the remaining 9 proteins which forms the target vector $\mathbf{y}_{n}\in\mathbb{R}^{7}$ lying on a graph of $M=7$ nodes. We perform the experiment 100 times in Monte Carlo simulation, where we random divide the total dataset into training and testing datasets. The average NMSE for testing datasets is shown in Figure 3 at SNR levels of $5$ dB and [math] dB. We once again observe the same trend that GPG outperforms GP for both linear and kernel cases at low sample sizes corrupted with noise.

IV-D Prediction for atmospheric tracer diffusion data

Our next experiment is on the atmospheric tracer diffusion measurements obtained from the European Tracer Experiment (ETEX) which tracked the concentration of perfluorocarbon tracers released into the atmosphere starting from a fixed location (Rennes, France)[62]. The observations were collected from two experiments over a span of 72 hours at 168 ground stations over Europe, giving two sets of 30 measurements, in total 60 measurements. Our goal is to predict the tracer concentrations on one half (84) of the total ground stations, using the concentrations at the remaining locations. The target signal $\mathbf{t}_{n}$ is then a graph signal over a graph of $M=84$ nodes where $n$ denotes the measurement index. The corresponding input vector $\mathbf{x}_{n}$ is also of length 84. We illustrate the ground station locations in a schematic in Figure 4 (a). The output nodes which correspond to the the target are shown in red markers with corresponding edges, whereas the rest of the markers denote the input. We simulate noisy training by adding white Gaussian noise at different SNR levels to the training data. We consider a geodesic distance based graph. The graph is constructed in the same manner like the graph used in the temperature data experiment previously. We randomly divide the total dataset of 60 samples equally into training and test datasets. We compute the NMSE for the different approaches by averaging over 100 different randomly drawn training subsets of size $N$ from the full training set of size $N_{ts}=30$. We plot the NMSE as a function of $N$ in Figures 4(b)-(c) at SNR levels of $5$dB and [math]dB. We observe that graph structure enhances the prediction performance signficantly under noisy and low sample size conditions.

V Reproducible research

In the spirit of reproducible research, all the codes relevant to the experiments in this article are made available at https://www.researchgate.net/profile/Arun_Venkitaraman and https://www.kth.se/ise/research/reproducibleresearch-1.433797.

VI Conclusions

We developed Gaussian processes for signals over graphs by employing a graph-Laplacian based regularization. The Gaussian process over graphs was shown to be a consistent generalization of the conventional Gaussian process, and that it provably results in a reduction of uncertainty in the output prediction. This in turn implies that the Gaussian process over graph is a better model for target vectors lying over a graph in comparison to the conventional Gaussian process. This observation is important in cases when the available training data is limited in both quantity and quality. Our expectation and motivation was that incorporating the graph structural information would help the Gaussian process make better predictions, particularly in absence sufficient and reliable training data. The experimental results with the real-world graph signals illustrated that this is indeed the case.