Sampling And Reconstruction Of Diffusive Fields On Graphs

TL;DR

This paper presents methods for reconstructing diffusive fields and localizing sources on graphs from limited, time-sampled observations, achieving exact recovery in noiseless cases and robustness with noise, with applications to physical systems.

Contribution

It introduces novel sampling and reconstruction techniques for diffusive fields on graphs, enabling source localization without sparsity constraints and handling noisy data.

Findings

Exact recovery possible with noiseless data

Robust performance of least squares estimators with noise

Application demonstrated on heat diffusion in a metal plate

Abstract

In this paper, the focus is on the reconstruction of a diffusive field and the localization of the underlying driving sources on arbitrary graphs by observing a significantly smaller subset of vertices of the graph uniformly in time. Specifically, we focus on the heat diffusion equation driven by an initial field and an external time-invariant input. When the underlying driving sources are modeled as an initial field or external input, the sources (hence the diffusive field) can be recovered from the subsampled observations without imposing any band-limiting or sparsity constraints. When the diffusion is induced by both the initial field and external input, then the field and sources can be recovered from the subsampled observations, however, by imposing band-limiting constraints on either the initial field or external input. For heat diffusion on graphs, we can compensate for the unobserved vertices with the temporal samples at the observed vertices. If the observations are noiseless, then the recovery is exact. Nonetheless, the developed least squares estimators perform reasonably well with noisy observations. We apply the developed theory for localizing and recovering hot spots on a rectangular metal plate with a cavity.

Equations

$${\mathbf{x}}_f={\mathbf{U}}^T\mathbf{x}\iff \mathbf{x}=\mathbf{U}{\mathbf{x}}_f.$$

$$\mathbf{H}=\mathbf{U}\mathrm{diag}[{\mathbf{h}}_f]{\mathbf{U}}^T\in {\mathbb{R}}^{N\times N},$$

$$\frac{\partial x(\mathbb{D},t)}{\partial t}=\alpha {\nabla }^{2}x(\mathbb{D},t)+q(\mathbb{D}),$$

$$\frac{\partial \mathbf{x}(t)}{\partial t}=-\mathbf{L}\mathbf{x}(t)+\mathbf{q},$$

$$\mathbf{x}(t)$$

$$=\mathbf{U}{e}^{-t\mathbf{\Lambda }}{\mathbf{U}}^T\mathbf{x}(0)+\mathbf{U}\left({\int }_0^t{e}^{-s\mathbf{\Lambda }}ds\right){\mathbf{U}}^T\mathbf{q}$$

$$=\mathbf{U}{e}^{-t\mathbf{\Lambda }}{\mathbf{x}}_f(0)+\mathbf{U}\left({\int }_0^t{e}^{-s\mathbf{\Lambda }}ds\right){\mathbf{q}}_f$$

$$f_t(\lambda )={\int }_0^t{e}^{-\lambda s}ds=\frac{1-e^{-t\lambda }}{\lambda }$$

$$\mathbf{x}(t)=\mathbf{U}\mathrm{diag}[{\mathbf{x}}_f(0)]\mathbf{a}(t)+\mathbf{U}\mathrm{diag}[{\mathbf{q}}_f]\mathbf{b}(t).$$

$$\mathbf{X}=\mathbf{U}\mathrm{diag}[{\mathbf{x}}_f(0)]{\mathbf{A}}^T+\mathbf{U}\mathrm{diag}[{\mathbf{q}}_f]{\mathbf{B}}^T,$$

$$\mathbf{Y}=[\mathbf{y}(t_1),\mathbf{y}(t_2),\cdots ,\mathbf{y}(t_T)]=\mathbf{\Phi }\mathbf{X}.$$

$$\mathbf{Y}=\mathbf{\Phi }\mathbf{X}=\mathbf{\Phi }\mathbf{U}\mathrm{diag}[{\mathbf{x}}_f(0)]{\mathbf{A}}^T.$$

$$\mathbf{y}=\mathrm{vec}(\mathbf{Y})=\left(\mathbf{A}\circ \mathbf{\Phi }\mathbf{U}\right){\mathbf{x}}_f(0),$$

$${\hat{\mathbf{x}}}_f(0)={\left(\mathbf{A}\circ \mathbf{\Phi }\mathbf{U}\right)}^{†}\mathbf{y},$$

$$\hat{\mathbf{x}}(0)=\mathbf{U}{\hat{\mathbf{x}}}_f(0).$$

$$\mathrm{rank}(\mathbf{A}\circ \mathbf{B})\ge \max \{\mathrm{rank}(\mathbf{A}),\mathrm{rank}(\mathbf{B})\}.$$

$$\mathbf{q}=\sum _{i=1}^P{\mathbf{u}}_i{q}_{f,i}={\mathbf{U}}_P{\mathbf{q}}_{f,P},$$

$${\boldsymbol{y}}={\rm vec}({\boldsymbol{Y}})=\left[\begin{array}[]{cc}{\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}&{\boldsymbol{B}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}\end{array}\right]\left[\begin{array}[]{c}{\boldsymbol{x}}_{f}(0)\\ {\boldsymbol{q}}_{f}\end{array}\right].$$

$$\mathbf{y}$$

$$\displaystyle=\left[\begin{array}[]{cc}{\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}&({\boldsymbol{B}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}){\boldsymbol{U}}^{T}{\boldsymbol{U}}_{P}\end{array}\right]\left[\begin{array}[]{c}{\boldsymbol{x}}_{f}(0)\\ {\boldsymbol{q}}_{f,P}\end{array}\right].$$

$$\left[\begin{array}[]{c}\widehat{{\boldsymbol{x}}}_{f}(0)\\ \widehat{{\boldsymbol{q}}}_{f,P}\end{array}\right]={\boldsymbol{\Psi}}^{\dagger}{\boldsymbol{y}},$$

$$\hat{\mathbf{x}}(0)=\mathbf{U}{\hat{\mathbf{x}}}_f(0);\hat{\mathbf{q}}={\mathbf{U}}_P{\hat{\mathbf{q}}}_{f,P}.$$

Full text

Index Terms—  Graph signal processing, graph sampling, heat diffusion, non-bandlimited signals, source localization on graphs.

1 Introduction

Graph signal processing extends tools from classical signal processing to deal with data defined on networks and other irregular domains [1, 2, 3]. We often come across such datasets in many diverse applications such as environmental sensing, traffic monitoring, mapping the human brain [4], cybersecurity [5], and social networks, to list a few.

Similar to how we understand many physical phenomena by a partial differential equation that explains the evolution of a spatiotemporal field and relates it to the inducing sources, we can also understand the temporal evolution of data over a network or an irregular domain using a partial differential equation. For example, the heat equation is often used to model the traffic movement, infection or virus spread, or rumor propagation [6, 5].

In this work, we focus on heat diffusion over networks. Specifically, we are interested in recovering diffusive signals on a graph by sampling a significantly smaller subset of vertices of the graph. This essentially amounts to localizing the underlying sources that drive the diffusion process from the observations that are collected at a few nodes. Oftentimes, the sources (e.g., traffic bottleneck or rumor sources) that induce the diffusion process are highly localized in the network or sparse in the vertex domain and hence are not usually bandlimited. Therefore, we require new sampling and recovery methods for graph signals that do not impose any structural or band-limiting constraints, unlike some of the existing graph sampling methods [7, 8, 9]. Although band-limiting constraints are not needed for recovering the second-order statistics of a signal defined on a graph from the subsampled observations, the framework developed in [10] cannot be used for localizing diffusive sources. Spatio-temporal sampling and reconstruction of diffusive fields on a regular domain under the assumption that the inducing sources are sparse are studied in [11, 12]. Assuming that the underlying sources are known, [13] focuses on estimating the time instance when the sources appear. In contrast, we will assume that the start time of the sources are is known.

In this work, we develop a graph sampling method to recover diffusive fields induced by an initial field and/or an external input that does not vary with time. The main results of this paper are as follows. When the underlying driving sources are modeled as an initial field or external input, we can localize and recover the sources by sampling a significantly smaller subset of vertices of the graph uniformly in time and by using a simple least squares estimator. To do so, we do not impose any constraints on the sources such as sparsity or bandlimitedness. Since we can compensate for the unobserved vertices with the temporal samples at the observed vertices, we can recover the sources without imposing any constraints. However, when the diffusion field is due to both the initial field and external input, to reconstruct the diffusive fields from the subsampled observations, we require either the initial field or external input to be bandlimited. If the observations are noiseless, then the recovery is exact. Nonetheless, the developed estimators perform reasonably well with noisy observations.

Throughout this paper, we will use upper (lower) case boldface letters to denote matrices (column vectors), and we will denote sets using calligraphic letters.

2 Graph signals

Consider an undirected graph ${\cal G}=\{{\cal V},{\cal E}\}$ with $N$ vertices (or nodes), where ${\cal V}=\{v_{1},v_{2},\ldots,v_{N}\}$ and ${\cal E}$ represent the vertex set and edge set, respectively. Let us denote the graph Laplacian matrix associated with ${\cal G}$ as ${\boldsymbol{L}}\in\mathbb{R}^{N\times N}$. A graph signal is a function $x:\mathcal{V}\rightarrow\mathbb{C}$ with $x(v)$ being the value of the function at vertex $v\in{\cal V}$. Let us collect the function values $\{x(v_{n})\}_{n=1}^{N}$ in a length-$N$ vector ${\boldsymbol{x}}=[x_{1},x_{2},\ldots,x_{N}]^{T}$.

For undirected graphs ${\boldsymbol{L}}$ is real symmetric, and hence admits an eigendecomposition ${\boldsymbol{L}}={\boldsymbol{U}}{\boldsymbol{\Lambda}}{\boldsymbol{U}}^{T}$ with ${\boldsymbol{U}}=[{\boldsymbol{u}}_{1},\cdots,{\boldsymbol{u}}_{N}]$ being the eigenvector matrix collecting the eigenvectors $\{{\boldsymbol{u}}_{n}\}_{n=1}^{N}$ and ${\boldsymbol{\Lambda}}={\rm diag}[\lambda_{1},\cdots,\lambda_{N}]$ being the diagonal matrix containing the corresponding eigenvalues $\{\lambda_{n}\}_{n=1}^{N}$. Here, $\operatorname*{\mathrm{diag}}[\cdot]$ refers to a diagonal matrix with its argument on the main diagonal. The eigenvectors and eigenvalues of ${\boldsymbol{L}}$ provide the notion of frequency in the graph setting [2, 3]. Specifically, $\{{\boldsymbol{u}}_{n}\}_{n=1}^{N}$ forms an orthonormal Fourier-like basis for graph signals with the graph frequencies denoted by $\{\lambda_{n}\}_{n=1}^{N}$. The graph Fourier transform of ${\boldsymbol{x}}$, denoted by ${\boldsymbol{x}}_{f}$, is given by

[TABLE]

We say that a graph signal ${\boldsymbol{x}}$ is bandlimited, if its graph Fourier transform ${\boldsymbol{x}}_{f}$ is sparse (i.e., contains a very few nonzero entries). Due to the uncertainty principle [14], a sparse graph signal ${\boldsymbol{x}}$ is not bandlimited in general.

The frequency content of graph signals may be modified using linear shift-invariant graph filters [15] of the form

[TABLE]

where ${\boldsymbol{h}}_{f}$ is the frequency response of the graph filter.

3 Data model

Let us consider a signal $x(\mathbb{D},t)$ in a physical domain $\mathbb{D}$ and temporal domain $t$. We will assume that $x(\mathbb{D},t)$ obeys the heat equation

[TABLE]

where $\nabla^{2}$ is the Laplace operator, $\alpha$ is the diffusion constant, and $q(\mathbb{D})$ is the external time-invariant input. When $t=0$, $x(0)=x(\mathbb{D},0)$ represents the initial field distribution. Without loss of generality, from now on we will assume $\alpha=-1$.

To solve such a differential equation on a surface or manifold, the manifold is discretized (e.g., using a Delaunay mesh), and the Laplace operator is replaced with a discrete Laplacian matrix (more specifically, a cotan-Laplacian matrix) denoted by ${\boldsymbol{L}}$. Thus, approximating (3) to

[TABLE]

where ${\boldsymbol{x}}(t)=[x_{1}(t),\ldots,x_{N}(t)]^{T}\in\mathbb{R}^{N}$ and ${\boldsymbol{q}}=[q_{1},\ldots,q_{N}]^{T}\in\mathbb{R}^{N}$ are signals defined on the graph represented by the Laplacian matrix ${\boldsymbol{L}}$. The differential equation (4) models heat diffusion on graphs, where the diffusion field is induced by ${\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}$.

The solution to the non-homogenous differential equation (4) is given by [16]

[TABLE]

where $e^{{\boldsymbol{L}}}={\boldsymbol{U}}e^{{\boldsymbol{\Lambda}}}{\boldsymbol{U}}^{T}\in\mathbb{R}^{N\times N}$ denotes the matrix exponential of ${\boldsymbol{L}}\in\mathbb{R}^{N\times N}$, ${\boldsymbol{x}}(0)$ is the initial field distribution at $t=0$. Here, ${\boldsymbol{x}}_{f}(0)={\boldsymbol{U}}^{T}{\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}_{f}={\boldsymbol{U}}^{T}{\boldsymbol{q}}$ are, respectively the graph Fourier transforms of ${\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}$. From (2), we can see that the diffusive field ${\boldsymbol{x}}(t)$ is obtained by filtering ${\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}$ with graph filters having frequency responses $e^{-t{\boldsymbol{\Lambda}}}$ and $\int_{0}^{t}e^{-s{\boldsymbol{\Lambda}}}\,\,ds$, respectively.

Let us introduce the vectors ${\boldsymbol{a}}(t)=[e^{-\lambda_{1}t},\ldots,e^{-\lambda_{N}t}]^{T}$ and ${\boldsymbol{b}}(t)=[f_{t}(\lambda_{1}),\ldots,f_{t}(\lambda_{N})]^{T}$, where

[TABLE]

with $f_{t}(0)=t$, and $f_{0}(\lambda)=0$. We can now express (5) compactly as

[TABLE]

Next, let us sample ${\boldsymbol{x}}(t)$ uniformly in time at instances $\{t_{k}=\Delta k,k=1,2,\cdots,T\}$ with step size $\Delta$ to obtain the data matrix ${\boldsymbol{X}}=[{\boldsymbol{x}}(t_{1}),{\boldsymbol{x}}(t_{2}),\cdots,{\boldsymbol{x}}(t_{T})]\in\mathbb{R}^{N\times T}$, which is given by

[TABLE]

where ${\boldsymbol{A}}=[{\boldsymbol{a}}(t_{1}),{\boldsymbol{a}}(t_{2}),\ldots,{\boldsymbol{a}}(t_{T})]^{T}\in\mathbb{R}^{T\times N}$ and ${\boldsymbol{B}}=[{\boldsymbol{b}}(t_{1}),{\boldsymbol{b}}(t_{2}),\ldots,{\boldsymbol{b}}(t_{T})]^{T}\in\mathbb{R}^{T\times N}$. Also, let us observe a subset of $K$ out of $N$ mesh points and denote this subset with ${\cal K}\subseteq{\cal V}$, where $|{\cal K}|=K$. By introducing a selection matrix ${\boldsymbol{\Phi}}\in\{0,1\}^{K\times N}$ that selects the field values at vertices indicated by ${\cal K}$, we can mathematically relate the subsampled observations to ${\boldsymbol{X}}$ as

[TABLE]

In what follows, we will develop estimators to recover ${\boldsymbol{x}}(0)$ and/or ${\boldsymbol{q}}(0)$ from ${\boldsymbol{Y}}$.

4 Diffusion field induced by ${\boldsymbol{x}}(0)$ or ${\boldsymbol{q}}$

In this section, we will develop a simple least squares estimator for reconstructing the diffusion field induced by ${\boldsymbol{x}}(0)$ or ${\boldsymbol{q}}$ from the subsampled data matrix ${\boldsymbol{Y}}$. More importantly, we do not impose any band-limiting constraints on the sources. This means that the sources may be sparse in the vertex domain and can model localized events such as rumor or infection sources in a complex network, traffic accidents in a road network, or diffusion of hot spots on a surface, to list a few .

Consider the case in which the diffusion field (5) is induced by only the initial field ${\boldsymbol{x}}(0)$ and the external input ${\boldsymbol{q}}={\bf 0}$. From (5) and (7), we have

[TABLE]

Vectorizing ${\boldsymbol{Y}}$, we get a system of $KT$ equations in $N$ unknowns given by

[TABLE]

where $\circ$ denotes the Khatri-Rao (i.e., columnwise Kronecker) product, ${\rm vec}(\cdot)$ refers to the matrix vectorization operator. Here, we have used the property ${\rm vec}({\boldsymbol{A}}\operatorname*{\mathrm{diag}}[{\boldsymbol{b}}]{\boldsymbol{C}})=({\boldsymbol{C}}^{T}\circ{\boldsymbol{A}}){\boldsymbol{b}}$.

Suppose we choose $K$ and $T$ such that $KT\geq N$, and if the matrix ${\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}$ has full-column rank, then we can estimate ${\boldsymbol{x}}_{f}(0)$ using least squares as

[TABLE]

and localize the sources as

[TABLE]

Using this in (5) allows us to compute the diffusive field at any time $t$ and at all the vertices. When the diffusion field is induced by ${\boldsymbol{q}}$ with ${\boldsymbol{x}}(0)={\bf 0}$, the least squares estimator for ${\boldsymbol{q}}$ may be developed along the similar lines.

The rank of the Khatri-Rao product of two matrices ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$ (of appropriate dimensions) with no all-zero column satisfies [17]

[TABLE]

When the sampling time instances $\{t_{1},\cdots,t_{T}\}$ and the eigenvalues of ${\boldsymbol{L}}$ are distinct, then the $T\times N$ Vandermonde matrix ${\boldsymbol{A}}$ will have full column rank of $N$ for $T\geq N$ and by construction does not have an all-zero column. Therefore, selecting rows of ${\boldsymbol{U}}$ such that there are no all-zero columns ensures that the rank of the matrix ${\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}$ will be $N$. In fact, observing only one node uniformly in time might result in the matrix ${\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}$ that has full column rank. However, in practice, depending on the observation time window, diffusion constant and the spectrum of ${\boldsymbol{L}}$, ${\boldsymbol{A}}$ might be ill-conditioned. In such cases, ${\boldsymbol{\Phi}}$ may be designed using sparse sensing (or sensor selection) techniques (e.g., see [18, 19]) to obtain a full column rank matrix ${\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}$.

5 Diffusion field induced by ${\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}$

In this section, we consider the case in which the diffusion field is induced by ${\boldsymbol{x}}(0)$ and a bandlimited time-invariant input ${\boldsymbol{q}}$, and provide a simple least squares estimator to recover the underlying sources from the subsampled data matrix ${\boldsymbol{Y}}$. Although we restrict ${\boldsymbol{q}}$ to be bandlimited, we do not impose any band-limiting or other structural constraints on ${\boldsymbol{x}}(0)$. The heat diffusion equation may be used to understand the movement of traffic in cities. Although the usual traffic movement may be assumed to be a smooth signal on a road network, there could exist a localized traffic bottleneck (e.g., due to an accident), which is a sparse non-bandlimited graph signal. Such diffusive fields may be modeled using (4) with a sparse ${\boldsymbol{x}}(0)$ representing the localized events and a bandlimited ${\boldsymbol{q}}$ representing the usual activity.

Recall that if ${\boldsymbol{q}}$ is bandlimited, then ${\boldsymbol{q}}_{f}$ will be sparse.Without loss of generality, let us assume that the first $P$ entries of ${\boldsymbol{q}}_{f}=[q_{f,1},q_{f,2},\ldots,q_{f,N}]^{T}$ are nonzero. Then, the bandlimited (or smooth) signal ${\boldsymbol{q}}$ may be expressed as a linear combination of the first few eigenvectors as

[TABLE]

where ${\boldsymbol{U}}_{P}\in\mathbb{R}^{N\times P}$.

Vectorizing ${\boldsymbol{Y}}$ in (7), we have

[TABLE]

Substituting (9), we get a linear system of $KT$ equations in $N+P$ unknowns

[TABLE]

If the matrix ${\boldsymbol{\Psi}}=\left[{\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}\quad({\boldsymbol{B}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}){\boldsymbol{U}}^{T}{\boldsymbol{U}}_{P}\right]$ has full column rank, which requires $KT\geq N+P$, we can use least squares to obtain

[TABLE]

and subsequently localize the underlying sources as

[TABLE]

In the previous subsection, we have seen that by appropriately selecting $K<N$ rows of ${\bf U}$ we may obtain a full column rank matrix ${\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}$ as ${\boldsymbol{A}}$ has full-column rank. As a consequence, by appropriately selecting $P>K$ rows of ${\bf U}$ will only increase the rank of $\left[{\boldsymbol{A}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}\quad{\boldsymbol{B}}\circ{\boldsymbol{\Phi}}{\boldsymbol{U}}\right]$ by $N+P$. This means that we have to impose some structural constraint on ${\boldsymbol{q}}$ to recover it uniquely from the subsampled data when the diffusive field is induced by both ${\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}$. In other words, by sampling in time and observing all the $N$ nodes, we can recover $2N$ unknowns ${\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}$ without any band-limiting constraints.

6 Numerical experiments

In this section, we apply the developed theory of graph sampling for reconstructing diffusive fields induced by hot spots on a metal block with a cavity. We use the partial differential equation toolbox from MATLAB to mesh the surface. The generated mesh with $N=134$ vertices is shown in Fig. 1. We observe $K=32$ vertices uniformly in time in the interval $[0,1.44]$s with step size $\Delta=0.16$s and $T=10$. The sampled vertices are also indicated in Fig. 1. We present results for the following two cases: (i) diffusive field induced by ${\boldsymbol{x}}(0)$, and (ii) diffusive field induced by ${\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}$.

As discussed in Section 4, to recover diffusive fields induced by ${\boldsymbol{x}}(0)$ with ${\boldsymbol{q}}={\bf 0}$, we do not require any band-limiting constraints. To demonstrate this, for ${\boldsymbol{x}}(0)$, we use a very sparse vector with only two non-zero entries at vertices $v_{88}$ and $v_{89}$. Since this initial field distribution is highly localized in the vertex domain, it is not bandlimited. Fig. 2(a) shows the initial field distribution at $t=0$, and Fig. 2(b) shows the evolution of the diffusive field at vertices $v_{73}$, $v_{88}$, $v_{89}$, and $v_{90}$ for different time instances. In Fig. 2(c), we can see the exact localization of the hot spots in the noiseless setting using a simple linear least squares estimator, and more importantly, without using any sparsity constraints. In Fig. 3, we consider a noisy setting in which the observations in (8) are corrupted with Gaussian noise having zero mean and variance $10^{-5}$. We show the normalized root mean squared error (RMSE), averaged over 1000 independent Monte-Carlo experiments, for different values of $K$. Although the error decreases as $K$ increases, we can see that increasing $T$ beyond a certain value does not lead to better performance. This is because ${\boldsymbol{A}}$ becomes ill-conditioned as $T$ increases.

For the case in which the diffusion field is induced due to both ${\boldsymbol{x}}(0)$ and ${\boldsymbol{q}}$, we use a sparse ${\boldsymbol{x}}(0)$ as before, and a bandlimited ${\boldsymbol{q}}$ with $P=5$. Fig. 2(d) shows the external time-invariant input, which is smooth on the surface. When ${\boldsymbol{q}}\neq{\bf 0}$, we can see in Fig. 2(e) that the field values do not decay with time as earlier. Fig. 2(f) shows the exact recovery of ${\boldsymbol{q}}$ using a simple linear least squares estimator (the reconstruction of ${\boldsymbol{x}}(0)$ is similar to Fig. 2(c), hence not shown), where we do not impose any sparsity constraints for recovering ${\boldsymbol{x}}(0)$. As before, gathering more samples in time does not lead to better performance as both ${\boldsymbol{A}}$ and ${\boldsymbol{B}}$ become ill-conditioned, and as a consequence we need to sample more vertices.

7 Concluding remarks

In this paper, we discussed the sampling and recovery of diffusive fields on graphs induced by possibly non-bandlimited sources. When the diffusion field is induced by an initial field or a time-invariant external input, we can localize and recover the sources by sampling a significantly smaller subset of nodes uniformly in time without imposing any band-limiting constraints and by using a simple least squares estimator. For diffusive fields induced due to an initial field and external input, we can exactly recover the sources from noiseless subsampled data when we constrain the external input to be bandlimited. When the observations are noiseless, the recovery is exact. In essence, for diffusion models on graphs, we can compensate for the unobserved vertices with the temporal samples at the observed vertices.