Robust Adaptive Estimation of Graph Signals Based on Welsch Loss

Abstract

This paper considers the problem of adaptive estimation of graph signals under the impulsive noise environment. The existing least mean squares (LMS) approach suffers from severe performance degradation under an impulsive environment that widely occurs in various practical applications. We present a novel adaptive estimation over graphs based on Welsch loss (WL-G) to handle the problems related to impulsive interference. The proposed WL-G algorithm can efficiently reconstruct graph signals from the observations with impulsive noises by formulating the reconstruction problem as an optimization based on Welsch loss. An analysis on the performance of the WL-G is presented to develop effective sampling strategies for graph signals. A novel graph sampling approach is also proposed and used in conjunction with the WL-G to tackle the time-varying case. The performance advantages of the proposed WL-G over the existing LMS regarding graph signal reconstruction under impulsive noise environment are demonstrated.

1. Introduction

1.1. Background and Motivation

The area of graph signal processing (GSP) has received extensive attention [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]. Graph signal processing has been found in social and economic networks, climate analysis, traffic patterns, marketing preferences, and so on [19,20,21,22,23,24,25,26]. The objective of GSP is to use the tools in DSP to the irregular domain in which the relationship between the elements are characterized via the graph. Under this framework, a signal occurring at graph nodes is handled over the graph topology. As an important GSP tool, the graph Fourier transform (GFT) has been introduced to decompose an observed graph signal into orthonormal components over the graph topology [10,12].

The graph sampling theory, one of the most important topics in GSP, aims to reconstruct the graph signals which are band-limited from the partial samples on the graph [27,28,29]. If the samples are not appropriately selected, the problem of recovering graph signals may be ill-conditioned. Therefore, optimizing the sampling set is critically important to the success of the graph signal recovery problem. The first method for sampling theory is proposed by Pesenson in [29]. On the other hand, since the adaptive algorithms are flexible [30,31,32,33,34,35], online graph signal reconstruction methods based on adaptive strategies have been proposed [36].

The major drawback of the aforementioned methods is that they undergo severe performance depression when encountering heavy-tailed impulsive noises, which are usually confronted in plenty of practical applications [37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54]. The aforementioned algorithms were proposed by employing the mean square error (MSE) method. The MSE-based algorithms may diverge when impulsive interference occurs [45]. Therefore, this urges for the development of a new algorithm to reconstruct graph signals under the impulsive noise environment. In summary, reconstructing the graph signal from the partial observations with impulsive noise is a significant problem. It is because that graph signal recovery and estimate have been widely studied with many promising applications. Applications contain power systems estimation [55,56,57,58], network time synchronization, and data registration [59,60].

1.2. Our Contributions

Our objective is to handle the problem of adaptive graph signal estimation under impulsive noise. The main aim of the paper is to meliorate a novel adaptive algorithm over graphs based on Welsch loss (WL-G). In contrast with the existing methods, which are predicated upon the MSE, the proposed WL-G converts the problem of the reconstruction for the graph signal to an optimization problem for a Welsch-loss-based cost function [61]. Unlike the MSE criterion, which undergoes severe performance depression in impulsive interference because of the property of the optimization based on $l_2$-norm [62], the Welsch loss is insensitive to the impulsive noise since the Welsch loss is a bounded nonlinear function that can eliminate large outliers. Therefore, the WL-G could effectively recover graph signals from imperfect observations under band-limited condition when impulsive noise occurring. The mean square performance is analyzed to highlight the importance of the selection of the sampling set. The results from this analysis are then exploited in the derivation of effective sampling strategies for the WL-G. We also take into account the case where the bandwidth and spectral contents are unknown and time-varying. To address this problem, we present an adaptive graph sampling (AGS) technique, which is used in conjunction with the WL-G to determine the signal support. The performance effectiveness of the proposed WL-G algorithm in impulsive noise is numerically demonstrated via various simulation examples.

The contributions of this paper are summarized as:

• We proposed a novel cost function on graph to deal with the impulsive noise environment.
• The detail analysis of the proposed algorithm is provided.
• The partial sampling strategy is proposed for WL-G algorithm.
• WL-G estimation with adaptive graph sampling is also considered to deal with the time-variant graphs.

2. Related Work

First, the papers regarding graph sampling without adaptive strategy are reviewed. Then, we overview graph sampling with adaptive strategy.

2.1. Graph Sampling without Adaptive Strategy

The graph sampling method based on Paley–Wiener was extended and developed in [19,63,64,65,66,67]. For instance, the prerequisite of individual restoration in GSP was provided in [27,66]. The authors also proposed the sampling method based on the Nyquist–Shannon theory. In [27], several greedy sampling methods on graphs were presented and the reconstruction (interpolation) performance of those schemes was also evaluated. Other effective sampling approaches in the graph spectral field were developed in [67].

2.2. Graph Sampling with Adaptive Strategy

Due to the merits of adaptive filter [68,69,70,71], the graph sampling with adaptive strategy has been proposed and developed in [36,72,73,74,75,76]. Specifically, Lorenzo et al. proposed an adaptive method to reconstruct graph signals by exploiting the use of the least mean squares (LMS) strategy [72]. The LMS approach is then extended to the distributed scheme in [74]. To speed up the adaptive estimation, the recursive least square (RLS) method is applied to the reconstruction of signals on graphs [75]. The probabilistic-sampling-based reconstruction approaches are then proposed by using the LMS and RLS strategies in [76]. In [77], to enhance convergence rate of the previous estimation algorithm of graph signals, the authors present two novel adaptive algorithms. First, the extended LMS (ELMS) was proposed, which extends the LMS algorithm via redeveloping the former vectors. To enhance the performance of ELMS, authors present the fast ELMS via using optimization for the gradient MSD. However, these graph sampling methods will suffer from severe performance deterioration when heavy-tailed impulsive noise occurs. The previous references regarding the graph sampling with adaptive strategy have been summarized in

Table 1. Graph sampling with adaptive strategy.

The Literature	Algorithm
[72]	LMS on graph
[74]	Distributed LMS on graph
[75]	RLS on graph
[76]	probabilistic LMS and RLS on graph
[77]	ELMS and FELMS on graph

.

3. Background of Graph Signal Processing

Consider an undirected graph $\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$ that is defined by two sets: the set of vertices $\mathcal{V}=\left\{1,2,\dots ,N\right\}$ and the set of weight edges $\mathcal{E}={\left\{a_{ij}\right\}}_{i,j\in \mathcal{V}}$, where $a_{ij}>0$ if nodes i and j are connected and $a_{ij}=0$ otherwise. Let $\mathbf{A}$ be the graph adjacency matrix whose ith entry is $a_{ij}$ representing the edge weight from node i to node j. For undirected graphs, $\mathbf{A}$ is symmetric. The degree matrix is a diagonal matrix $\mathbf{K}$, whose $ith$ diagonal entry is expressed as $k_i={\displaystyle \sum _{j=1}^N}{a}_{ij}$.

As the most fundamental operator in GSP, the graph Laplacian matrix takes the form $\mathbf{L}:=\mathbf{K}-\mathbf{A}$. Obviously, for undirected graphs, $\mathbf{L}$ is also positive semidefinite and symmetric. Therefore, $\mathbf{L}$ can be eigendecomposed. Given the eigendecomposition of $\mathbf{L}$ as $\mathbf{U}\Lambda {\mathbf{U}}^T$ with $\mathbf{U}=\left[{\mathit{u}}_1,\cdots {\mathit{u}}_N\right]$ and $\Lambda =\mathrm{diag}\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_N\right)$, the ascending-order set of eigenvalues $0<{\lambda }_1<{\lambda }_2<\cdots <{\lambda }_N$ denote the graph frequencies.

A signal $\mathit{x}:\mathcal{V}\to {\mathbb{R}}^N$ is defined as an $N\times 1$ vector whose ith entry $x_i$ represents the vertex value of node i. For graph signal $\mathit{x}$, the definition of the GFT is given by

$$\begin{array}{c}\hfill \mathit{s}={\mathbf{U}}^T\mathit{x}.\end{array}$$

(1)

Alternative definitions for the GFT are also available, e.g., via the use of adjacency matrix [10]. We focus on the definition dependent on the Laplacian matrix in this paper.

In many cases, the graph signal is usually band-limited and can be given by

$$\begin{array}{c}\hfill \mathit{x}=\mathbf{U}\mathit{s}\end{array}$$

(2)

where the GFT $\mathit{s}$ is sparse. The support of $\mathit{s}$ is defined as

$$\begin{array}{c}\hfill \mathcal{F}=\left\{i\in \left\{1,\dots ,N\right\}:s_i\ne 0\right\}.\end{array}$$

(3)

The cardinality of $\mathcal{F}$, i.e., $\left|\mathcal{F}\right|$, leads to the bandwidth of $\mathit{x}$. We denote ${\mathbf{D}}_{\mathcal{S}}$ as the vertex limiting operator:

$$\begin{array}{c}\hfill {\mathbf{D}}_{\mathcal{S}}=\mathrm{diag}\left\{1_{\mathcal{S}}\right\}\end{array}$$

(4)

where $1_{\mathcal{S}}$ denotes a vector and $\mathcal{S}$ denotes a subset of $\mathcal{V}$, i.e., $\mathcal{S}\subseteq \mathcal{V}$. If $i\notin \mathcal{S}$, the ith element of $1_{\mathcal{S}}$ becomes 0. For $i\in \mathcal{S}$, the ith element of $1_{\mathcal{S}}$ is 1. The band-limiting operator over the set $\mathcal{F}$ that satisfies $\mathcal{F}\subseteq \mathcal{V}$ is defined as

$$\begin{array}{c}\hfill {\mathbf{B}}_{\mathcal{F}}=\mathbf{U}{\Sigma }_{\mathcal{F}}{\mathbf{U}}_{}^T\end{array}$$

(5)

where ${\Sigma }_{\mathcal{F}}$ denotes a diagonal matrix whose ith diagonal element is zero, if $i\notin \mathcal{F}$, and one otherwise. Note that both matrices ${\mathbf{D}}_{\mathcal{S}}$ and ${\mathbf{B}}_{\mathcal{F}}$ are idempotent and self-adjoint.

4. Adaptive WL-G Estimation on Graphs

Given a graph signal ${\mathit{x}}^o={\left\{x_{oi}^{}\right\}}_{i=1}^N\in {\mathbb{R}}^N$ defined over the graph $\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$, we assume that

Assumption 1. 

(band-limited):  The signal ${\mathit{x}}^o$ is $\mathcal{F}$-bandlimited.

The observed signal is corrupted by the additive noise $\mathit{v}\left[n\right]$ including the background noise $\mathbf{\eta }\left[n\right]$ with covariance matrix ${\mathit{C}}_{\eta }$ and the impulsive noise $\mathbf{\xi }\left[n\right]$ with covariance matrix ${\mathit{C}}_{\xi }$. Therefore, the observed noisy signal has the following form at each time n

$$\begin{array}{c}\hfill \mathit{y}\left[n\right]=\mathbf{D}\left({\mathit{x}}^o+\mathit{v}\left[n\right]\right)=\mathbf{D}\mathbf{B}{\mathit{x}}^o+\mathbf{D}\mathit{v}\left[n\right]\end{array}$$

(6)

where the $\mathbf{D}$ is given by (4). For convenience, we omit the subscript in ${\mathbf{D}}_{\mathcal{S}}$. The covariance matrix of $\mathit{v}\left[n\right]$ can be denoted as ${\mathit{C}}_v$. The objective of the graph signal reconstruction under consideration is to recover ${\mathit{x}}^o$. The previous graph signal reconstruction algorithms were developed based on the mean-squared error criterion under the assumption of Gaussian noise. However, the MSE-based approaches are unstable to outliers and thus performs rather poorly under the presence of impulsive noise. To make the recovering process robust against the impulsive noise, we use the cost function based on Welsch loss [61]

$$\begin{array}{c}\hfill \mathit{f}\left(\mathit{x}\left[\mathit{n}\right]\right)=1-\sum _{i=1}^{N}exp(-\frac{e_i^2\left[n\right]}{2c^2})\end{array}$$

(7)

where c denotes a scale parameter that restraints the scale of the quadratic bowl of loss and $e_i\left[n\right]$ is the ith entry of $\mathit{e}\left[n\right]$ with $\mathit{e}\left[n\right]=\mathit{y}\left[n\right]-\mathbf{DB}\mathit{x}\left[n\right]$. The Welsch loss given by the kernel width verified to be effective in suppressing the impacts of impulsive noise or large outliers. Figure 1 depicts the cost function (7) and its first derivative to demonstrate the robustness of the arctangent function.

An estimate for ${\mathit{x}}^o$ can be obtained by

$$\begin{array}{cc}\hfill & \underset{\mathit{x}\left[n\right]}{min}E\left\{\mathit{f}\left(\mathit{x}\left[n\right]\right)\right\}\hfill \\ \hfill & s.t.\mathbf{B}\mathit{x}\left[n\right]=\mathit{x}\left[n\right]\hfill \end{array}$$

(8)

This optimization problem is robust against the impulsive noise because the objective function grounded in the Welsch loss has been proven to be insensitive to the impulsive noise [45,78,79]. In contrary to the sign algorithm, which can be derived from the $l_1$-norm-based cost function, the Welsch loss is differentiable and mathematically tractable.

A common method to address optimization problem (8) is stochastic steepest-descent method, which is expressed as [62]

$$\begin{array}{c}\hfill \mathit{x}[n+1]=\mathit{x}\left[n\right]+\mu \mathbf{B}\mathbf{D}\varphi \left(\mathit{y}\left[n\right]-\mathit{x}\left[n\right]\right)\end{array}$$

(9)

where $\mu >0$ denotes the convergence rate, $\varphi \left(\mathit{y}\left[n\right]-\mathit{x}\left[n\right]\right)$ denotes a vector whose ith element is given by $\varphi \left(y_i\left[n\right]-x_i\left[n\right]\right)=exp(-\frac{(y_i\left[n\right]-x_i{\left[n\right])}^{\mathbf{2}}}{2c^2}\left(y_i\right[n]-x_i\left[n\right])$ and the second equation is relied on the fact that D is an idempotent operator. To satisfy the constraint in (8), the following equation should hold

$$\mathit{x}[n+1]=\mathbf{B}\mathit{x}[n+1]$$

(10)

Using (9), we have

$$\mathbf{B}\mathit{x}[n+1]=\mathbf{B}\{\mathit{x}\left[n\right]+\mu \mathbf{B}\mathbf{D}\varphi \left(\mathit{y}\left[n\right]-\mathit{x}\left[n\right]\right)\}=\mathbf{B}\mathit{x}\left[n\right]+\mu \mathbf{B}\mathbf{D}\varphi \left(\mathit{y}\left[n\right]-\mathit{x}\left[n\right]\right)$$

(11)

According to (10), to make $\mathit{x}[n+1]=\mathbf{B}\mathit{x}[n+1]$ hold, $\mathit{x}\left[n\right]$ should equal to $\mathbf{B}\mathit{x}\left[n\right]$. Then, we have the following equations

$$\mathit{x}[n-1]=\mathbf{B}\mathit{x}[n-1],\mathit{x}[n-2]=\mathbf{B}\mathit{x}[n-2],\cdots ,\mathit{x}\left[1\right]=\mathbf{B}\mathit{x}\left[1\right]$$

(12)

To make (12) hold, we have

$$\mathit{x}\left[0\right]=\mathbf{B}\mathit{x}\left[0\right]$$

(13)

In the other words, when the initialization of $\mathit{x}\left[n\right]$ satisfies $\mathit{x}\left[0\right]=\mathbf{B}\mathit{x}\left[0\right]$, the algorithm (9) satisfies the constraint in (8). The algorithm (9) is referred to as the WL-G algorithm. It is important to note that the properties of the proposed WL-G algorithm is strongly dependent on the sampling matrix $\mathbf{D}$ and the band-limiting operator matrix $\mathbf{B}$. Therefore, to optimize the WL-G method, a mean square analysis is presented in the next part to characterize the dependence of the WL-G performance on the choice of sampling matrix D. Using these outcomes, the effective sampling strategies for the WL-G are developed in Section 6. An adaptive graph sampling method is then presented in Section 7.

Computational Complexity

The computational complexity of the LMS, LMP, and WL-G algorithms is provided in

Table 2. Graph sampling with adaptive strategy.

Algorithm	Multiplications	Additions	p-Norm	Exponent
LMS on graph	$N(N+1)$	$2N^2$	-	-
LMP on graph	$N(N+1)$	$2N^2$	N	-
WL-G	$2N^2$	$2N^2+2N$	-	N

. Compared with LMS, the additional complexity of the WL-G arises from the calculation of function $\varphi$. Significantly, in spite of growth complexity, the WL-G surpasses the LMS and LMP algorithm.

5. Mean Square Analysis

Rewriting the WL-G algorithm as follows

$$\begin{array}{c}\hfill x_i[n+1]=x_i\left[n\right]+\mu \sum _{j=1}^N{b}_{ij}{d}_j\varphi \left(y_j\left[n\right]-x_j\left[n\right]\right),i=1,2,\cdots ,N\end{array}$$

(14)

We then define the error vector $\tilde{\mathit{x}}\left[n\right]=\mathit{x}\left[n\right]-{\mathit{x}}^o$ whose ith entry is ${\tilde{x}}_i\left[n\right]$.

Lemma 1. 

Using (14), the recursion of error vector $\tilde{\mathit{x}}\left[n\right]$ is given by

$$\begin{array}{c}\hfill \tilde{\mathit{x}}[n+1]=\left(\mathbf{I}-\mathbf{B}\mathbf{D}\mathcal{D}\mathbf{B}\right)\tilde{\mathit{x}}\left[n\right]+\mu \mathbf{B}\mathbf{D}{\mathit{v}}_n^g\left(\mathit{x}\left[n\right]\right).\end{array}$$

(15)

where${\mathit{v}}_n^g\mathit{x}\left[n\right]$is defined in Appendix A.

Proof. 

See Appendix A. □

Multiplying both sides of (15) by ${\mathbf{U}}^T$ and using (5), we obtain

$$\begin{array}{c}\hfill \tilde{\mathit{s}}[n+1]=\left(\mathbf{I}-\Sigma {\mathbf{U}}^T\mathbf{D}\mathcal{D}\mathbf{U}\Sigma \right)\tilde{\mathit{s}}\left[n\right]+\mu \Sigma {\mathbf{U}}^T\mathbf{D}{\mathit{v}}_n^g\left(\mathit{x}\left[n\right]\right)\end{array}$$

(16)

where $\tilde{\mathit{s}}\left[n\right]={\mathbf{U}}^T\tilde{\mathit{x}}\left[n\right]$ is the GFT of $\tilde{\mathit{x}}\left[n\right]$. We only consider $\widehat{\mathit{s}}\left[n\right]=\left[\tilde{\mathit{s}}\left[n\right],i\in \mathcal{F}\right]\in {\mathbb{R}}^{\left|\mathcal{F}\right|}$. The error recursion (16) becomes

$$\begin{array}{c}\hfill \widehat{\mathit{s}}[n+1]=\left(\mathbf{I}-\Sigma {\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}\right)\widehat{\mathit{s}}\left[n\right]+\mu {\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}{\mathit{v}}_n^g\left(\mathit{x}\left[n\right]\right)\end{array}$$

(17)

where $\widehat{\mathit{s}}\left[n\right]={\mathbf{U}}_{\mathcal{F}}^T\tilde{\mathit{x}}\left[n\right]$.

Using (17), the recursion of $E\left\{{∥\widehat{\mathit{s}}[n+1]∥}_{\mathit{c}}^2\right\}$ can be derived as follows.

Lemma 2. 

$E\left\{{∥\widehat{\mathit{s}}[n+1]∥}_{\mathit{c}}^2\right\}$can be calculated as

$$\begin{array}{c}\hfill E\left\{{∥\widehat{\mathit{s}}[n+1]∥}_{\mathit{c}}^2\right\}=E\left\{{∥\widehat{\mathit{s}}\left[n\right]∥}_{\mathbf{Q}\mathit{c}}^2\right\}+{\mu }^2\mathit{vec}{\left(\mathbf{G}\right)}^T\mathit{c}\end{array}$$

(18)

where$\mathbf{D}$and$\mathbf{Q}$are given in Appendix B.

Proof. 

See Appendix B. □

According to Lemma 2, we obtain two theorems regarding steady-state mean square deviation (MSD) and convergence stability as follows.

Theorem 1.

Assuming that the data model (6) and Assumption 1 hold, the steady state of the WL-G algorithm is given by

$$\begin{array}{c}\hfill \mathrm{MSD}={\mu }^{2}vec{\left(\mathbf{G}\right)}^T{\left(\mathbf{I}-\mathbf{Q}\right)}^{-1}vec\left(\mathbf{I}\right).\end{array}$$

(19)

Proof. 

In steady-state, and assuming that the matrix $\mathbf{I}-\mathbf{Q}$ is invertible, from the recursive expression (18), we obtain

$$\begin{array}{c}\hfill \underset{n\to \infty }{lim}E{∥\widehat{\mathit{s}}\left[n\right]∥}_{\left(\mathbf{I}-\mathbf{Q}\right)\mathit{c}}^2={\mu }^2\mathit{vec}{\left(\mathbf{G}\right)}^T\mathit{c}.\end{array}$$

(20)

Let $\mathit{c}={\left(\mathbf{I}-\mathbf{Q}\right)}^{-1}\mathit{vec}\left(\mathbf{I}\right)$, we obtain

$$\begin{array}{c}\hfill \mathrm{MSD}={\mu }^{2}vec{\left(\mathbf{G}\right)}^T{\left(\mathbf{I}-\mathbf{Q}\right)}^{-1}vec\left(\mathbf{I}\right).\end{array}$$

(21)

The theorem is proved. □

Theorem 2.

Assume that Assumption 1 and (6) hold, the WL-G can converge if the step-size μ

$$\begin{array}{c}\hfill 0<\mu <\frac{2}{{\lambda }_{max}\left({\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}\right)}.\end{array}$$

(22)

Proof. 

Iterating recursion (18) starting from $n=0$, we find that

$$\begin{array}{c}\hfill E{∥\widehat{\mathit{s}}[n+1]∥}_{\mathit{c}}^2=E{∥\widehat{\mathit{s}}\left[0\right]∥}_{{\mathbf{Q}}^{n+1}\mathit{c}}^2+{\mu }^2\mathit{vec}{\left(\mathbf{G}\right)}^T\sum _{i=0}^n{\mathbf{Q}}^i\mathit{c}\end{array}$$

(23)

with initial condition $\widehat{s}\left[0\right]$. Note $\mathbf{Q}$ is stable when $\mathbf{I}-{\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}$ is stable. Thus, we require $0<\left|1-\mu {\lambda }_{max}\left({\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}\right)\right|<1$. After some algebraic manipulations, (22) can be obtained. □

6. Sampling Strategy

According to (18), (19) and (22), the performance of the WL-G algorithm relies on the vertex limiting operator $\mathbf{D}$. Therefore, sampling signals defined on graphs is not only about choosing the number of samples but also about (if possible) having an appropriate strategy to optimally determine where to sample because the sampling locations has a great influence on the performance of the WL-G.

The objective is to determine the optimal sampling set (i.e., the vertex limiting operator $\mathbf{D}$), which minimizes the value of MSD for sampling strategy. Assuming that the matrix $\mathbf{Q}$ can be eigendecomposed by $\mathbf{Q}=\mathbf{V}\Pi {\mathbf{V}}^T$, the MSD in (19) can be rewritten as

$$\begin{array}{cc}\hfill & \mathrm{MSD}={\mu }^2\mathit{vec}{\left(\mathbf{T}\right)}^T{\left(\mathbf{I}-\mathbf{V}\Pi {\mathbf{V}}^T\right)}^{-1}\mathit{vec}\left(\mathbf{I}\right)\hfill \\ \hfill & ={\mu }^2\mathit{vec}{\left(\mathbf{T}\right)}^T\mathbf{V}{\left(\mathbf{I}-\Pi \right)}^{-1}{\mathbf{V}}^T\mathit{vec}\left(\mathbf{I}\right)\hfill \\ \hfill & ={\mu }^2\sum _{i=1}^{{\left|\mathcal{F}\right|}^2}\frac{p_i{q}_i}{1-{\lambda }_i\left(\mathbf{Q}\right)}.\hfill \end{array}$$

(24)

where $q_i$ and $p_i$ are the ith terms of the vectors ${\mathbf{V}}^T\mathit{vec}\left(\mathbf{I}\right)$ and ${\mathbf{V}}^T\mathit{vec}\left(\mathbf{T}\right)$, respectively. In order to get as a low MSD value as possible, the matrix $\mathbf{Q}$ should be selected such that its eigenvalues

$${\lambda }_i\left(\mathbf{Q}\right)=\left(1-\mu {\lambda }_k\left({\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}\right)\right)\left(1-\mu {\lambda }_l\left({\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}\right)\right)$$

are as far as from 1, where $k,l=1,2,\cdots ,\left|\mathcal{F}\right|$. In other words, the eigenvalues of the matrix ${\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}$ should differ from 0 as much as possible. As a result, we can use the greedy approximation method to obtain an approximate minimization of (24). The main idea is to iteratively select the samples from the graph that maximize the pseudodeterminant of the matrix ${\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}$, denoted by $\mathrm{Det}\left({{\mathbf{U}}_{\mathcal{F}}^T\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}}_{}\right)$. Using the greedy method results in the sampling strategy (called the greedy determinant–maximization method) as given in

Table 3. Maximizing $\mathrm{Det}\left({\mathbf{U}}_{\mathcal{F}}^{\mathrm{T}}\mathbf{D}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}\right)$.

Inputdata:M

Outputdata:$\mathcal{S}$

Initialization:$\mathcal{S}\equiv ⌀$

Function:

while $\left|\mathcal{S}\right|<M$

$\hspace{1em}\hspace{1em}s=\underset{j}{argmax}\mathrm{Det}\left({{\mathbf{U}}_{\mathcal{F}}^T{\mathbf{D}}_{\mathcal{S}\cup \left\{j\right\}}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}}_{}\right)$

$\hspace{1em}\hspace{1em}\mathcal{S}=\mathcal{S}\cup \left\{s\right\}$

end

. Since the computational steps of the greedy determinant–maximization sampling strategy have high computation complexity, it is desirable to reduce the computational complexity while at the same not significantly compromising the reconstruction performance. Motivated by the partial method in [80,81], we propose the partial greedy approach as summarized in

Table 4. Partial greedy determinant–maximization sampling strategy.

Inputdata:M

Outputdata:$\mathcal{S}$

Initialization:$\mathcal{S}\equiv ⌀$

Function:

for $i=1:M$

$\hspace{1em}\hspace{1em}{\mathcal{K}}_i=\{1,\dots ,N\}$

$\hspace{1em}\hspace{1em}{\mathcal{T}}_i=\mathrm{randit}(1,(1-p)N,[1,N\left]\right)$

$\hspace{1em}\hspace{1em}s=\underset{j\in {\mathcal{K}}_i\backslash {\mathcal{T}}_i}{argmax}\mathrm{Det}{\left({{\mathbf{U}}_{\mathcal{F}}^T{\mathbf{D}}_{\mathcal{S}\cup \left\{j\right\}}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}}_{}\right)}_{+}$

$\hspace{1em}\hspace{1em}\mathcal{S}=\mathcal{S}\cup \left\{s\right\}$

end

, where p is the partial probability and $\mathrm{randit}(1,(1-p)N,[1,N\left]\right)$ denotes $(1-p)N$ random integer sets between 1 and N.

In

Table 5. Framework of proposed WL-G algorithm.

Problem: Recovering the band-limited graph signal ${\mathit{x}}^o$ from partial observations $\mathit{y}\left[n\right]=\mathbf{D}\mathbf{B}{\mathit{x}}^o+\mathbf{D}\mathit{v}\left[n\right]$ with impulsive noise.

Inputdata:M, $\mathit{y}\left[n\right]$,$\mathbb{N}$

Outputdata:$\mathit{x}\left[n\right]$

Initialization:$\mathcal{S}\equiv ⌀$, $\mathit{x}\left[0\right]$

Function:

while $\left|\mathcal{S}\right|<M$

$\hspace{1em}\hspace{1em}s=\underset{j}{argmax}\mathrm{Det}\left({{\mathbf{U}}_{\mathcal{F}}^T{\mathbf{D}}_{\mathcal{S}\cup \left\{j\right\}}\mathcal{D}{\mathbf{U}}_{\mathcal{F}}^{}}_{}\right)$

$\hspace{1em}\hspace{1em}\mathcal{S}=\mathcal{S}\cup \left\{s\right\}$

end

Using $\mathcal{S}$ to calculate $\mathbf{D}$

while $n<\mathbb{N}$

$\hspace{1em}\hspace{1em}\mathit{x}[n+1]=\mathit{x}\left[n\right]+\mu \mathbf{B}\mathbf{D}\varphi \left(\mathit{y}\left[n\right]-\mathit{x}\left[n\right]\right)$

end

, we provide the framework of proposed WL-G estimation algorithm. $\mathbb{N}$ in Table 5 denotes the iteration number.

7. Wl-G Estimation with AGS

The band-limiting operator $\mathbf{B}$ is assumed in (9) to be known in advance. However, since the graph may vary over time, the prior knowledge of $\mathbf{B}$ is sometimes unrealistic in certain applications. To overcome the lack of prior knowledge of $\mathbf{B}$, we present an AGS technique for the WL-G algorithm.

Using (1), the graph signal observation model can be rewritten as

$$\begin{array}{c}\hfill \mathit{y}\left[n\right]=\mathbf{D}\mathbf{U}{\mathit{s}}^o+\mathbf{D}\mathit{v}\left[n\right].\end{array}$$

(25)

As a result, estimating ${\mathit{x}}^o$ is equivalent to estimating ${\mathit{s}}^o$. Motivating by the fact that the support identification for ${\mathit{s}}^o$ is acutely associated with the set of sampling, the optimal problem (8) could be recast as

$$\begin{array}{cc}\hfill & \underset{\mathit{s},\mathbf{D}}{min}E\left\{\mathit{f}\left(\mathit{s}\right)\right\}+\lambda g\left(\mathit{s}\right)\hfill \end{array}$$

(26)

where $g(·)$ is a $l_0-$norm or $l_1-$ norm, and $\lambda$ denotes a constant that regulates the sparsity level of $\mathit{s}$. The minimization problem (26) is a nonconvex program; thus, it is generally challenging to be solved. The ISTA [82] is adopted to solve (26). Provided $\mathbf{D}\left[n\right]$, $\mathit{s}$ is updated through the improved ISTA algorithm as

$$\begin{array}{cc}\hfill \mathit{s}[n+1]& =T_{\gamma }(\mathit{s}\left[n\right]+\mu \frac{\partial \mathit{f}\left(\mathit{s}\right)}{\partial \mathit{s}\left[n\right]})\hfill \\ \hfill & =T_{\gamma }(\mathit{s}\left[n\right]+\mu {\mathbf{U}}^T\mathbf{D}\left[n\right]\varphi \left(\mathit{y}\left[n\right]-\mathbf{U}\mathit{s}\left[n\right]\right))\hfill \end{array}$$

(27)

where $T_{\gamma }$ is a thresholding function, $\mu$ denotes a step-size. Here, the thresholding function is selected as a hard threshold given by

$$\begin{array}{c}\hfill T_{\gamma }\left(s_m\right)=\left\{\begin{array}{cc}{s}_m& \left|s_m\right|>\gamma \\ 0& \left|s_m\right|<\gamma .\end{array}\right.\end{array}$$

(28)

where $T_{\gamma }\left(s_m\right)$ is the mth element of the thresholding function ${\mathit{T}}_{\gamma }\left(\mathit{s}\right)$.

8. Simulation

This part carries out simulation examples to demonstrate the advantages of the proposed algorithm on graphs and verify the correctness of the performance analysis. A graph signal with 50 nodes is considered. The graph is depicted in Figure 2.

8.1. On the Theoretical Results

We now verify the precision of the mean square analytical theory provided in Theorem 1 by comparing the analytical MSD values with the MSD obtained via simulations for both the cases of nonimpulsive noise and impulsive noise. $\left|\mathcal{F}\right|$ is set to 10.

Example 1. 

(No impulsive Noise):In this scenario, the observation noise only contains the background noise. The background noise $\eta \left[n\right]$ is the zero-mean Gaussian process with diagonal covariance matrix,

$${\mathit{C}}_{\eta }=\left[\begin{array}{cccc}{c}_{\eta ,1}^2& & & \\ & c_{\eta ,2}^2& & \\ & & \ddots & \\ & & & c_{\eta ,N}^2\end{array}\right]$$

(29)

where $c_{\eta ,i}^2$ is generated uniformly randomly between 0 and 0.01. Figure 3 shows the simulated MSD of the WL-G compared with the theoretical MSD value given by Theorem 1 for various values of the step sizes $\mu =0.001,0.003,0.006,0.008$. Here, $\left|\mathcal{S}\right|$ is set to 10. The steady-state MSD of WL-G gained from the experiments closely match the analytical MSD value given by Theorem 1. This thus confirms the correctness of the mean square analysis given by Theorem 1. In addition, it is observed that, although providing a higher learning rate (convergence rate) of the WL-G algorithm, a larger μ results in an inferior property. In contrast, a smaller value of step size results in a better MSD performance of the WL-G algorithm but with a slower learning rate.

Example 2. 

(Impulsive Noise):This example considers the case of impulsive noise where the impulsive noise $\xi \left[n\right]$ is generated from the Bernoulli–Gaussian (BG) process, which is often used in the performance analysis. By denoting ${\xi }_i\left[n\right]$ as the ith element of $\mathbf{\xi }\left[n\right]$, we have ${\xi }_i\left[n\right]=b_i\left[n\right]g_i\left[n\right]$, where $b_i\left[n\right]$ is a Bernoulli process whose probabilities are expressed by $P\left[b_i\left[n\right]=1\right]=P_i$, and $g_i\left[n\right]$ is a zero-mean white Gaussian process with variance $\kappa c_{\eta ,i}^2$. Note that $p_{r,i}$ stands for the probability of occurrence of impulsive noise sample. In the simulation, the probability $p_{r,i}$ is set to 0.05, and the constant κ is set to 10,000. In addition to the impulsive noise $\xi \left[n\right]$, the background noise $\eta \left[n\right]$ is generated in the same manner as in Example 1. Figure 4 compares the transient behavior of the MSD of the WL-G algorithm with the theoretical steady-state MSD value given in Theorem 1 for the step sizes of $\mu =0.0008,0.001,0.003,0.006$. A similar observation as in Example 1 is made here, where the theoretical steady-state MSD value derived in Theorem 1 closely agrees with the simulated steady-state MSD. This once again verifies the correctness of Theorem 1.

8.2. On the Performance of The WL-G Algorithm

Example 3. 

(Effect of Cardinality $\left|\mathcal{S}\right|$):This example examines how the behavior of the WL-G algorithm is affected by the number of samples in the observation set $\mathcal{S}$. Figure 5 plots the transient MSD of the WL-G algorithm for $\left|\mathcal{S}\right|=10,20,30$, and 40. Here, the simulation setup is the same as in Example 2 except that the step size is set to $0.06$. We can see from Figure 5 that, as expected, growing the samples improve the convergence speed.

Example 4. 

(WL-G Versus LMS):In this example, the behavior of WL-G algorithm is compared with those of the LMS algorithm [72] and the LMP algorithm on graph. Here, the cardinality $\left|\mathcal{S}\right|$ is set to 10. The probabilities $P_i$ are set to 0.05, 0.01, and 0.005. Other parameters remain the same as those in Example 3.

Figure 6 compares the WL-G, LMP, and LMS on graph in the impulsive noise background. It is observed that the LMS algorithms yields a poor MSD performance due to the impulsive noise. It is as expected as the LMS algorithm is in virtue of MSE criterion which is unstable to outliers. In contrast, the WL-G and the LMP algorithm on graph can effectively cope with the impulsive noise by producing a much more reliable MSD performance.

The MSD curves of the WL-G and LMS algorithm in the Guassian noise are depicted in Figure 7. As can be seen from Figure 7, we obtain that the WL-G algorithm exhibits nearly same performance as the LMS algorithm.

Figure 8, Figure 9 and Figure 10 shows the measured graph signal and the reconstructed graph signals obtained by the LMS and WL-G algorithms. The WL-G algorithm is observed to produce a reconstructed signal that is almost identical to the ground-truth signal in Figure 2. On the other hand, the LMS algorithm results in an unsatisfactory reconstructed signal where the signal values at most of the nodes are far different from the corresponding true values. These observations demonstrate one more time the performance advantages of the proposed WL-G algorithm over the LMS.

Example 5. 

(Effect of scale parameter $\mathit{c}$):This example examines how the performance of the WL-G algorithm is affected by the scale parameter c. Figure 8 and Figure 9 depict the MSD curves of the WL-G algorithm with different scale parameters c under impulsive noise, where $p_{r,i}$ are set to 0.05 and 0.5 in Figure 11 and Figure 12, respectively. Smaller c leads to smaller steady-state error but lower convergence rate. Therefore, we can select a suitable value for c according to the specific requirement.

Example 6. 

(WL-G (greedy) Versus WL-G (partial greedy)):We conduct the simulation to test the partial greedy approach, where p is set to $0.8$. The simulated result is depicted in Figure 13. From the results of the simulation, we can determine that the algorithm in Table 4 has the similar performance to the algorithm in Table 3.

8.3. On the Performance of WL-G Algorithm with Adaptive Graph Sampling

In this section, we present a performance evaluation for the proposed WL-G algorithm with adaptive graph sampling. A time-varying graph signal with $N=50$ nodes is considered, where the spectral content of the signal switches between the first 5, 15, and 10 eigenvectors. The graph topology is the same as in Figure 3. The elements of the GFT ${\mathit{s}}^o$ inside the support are set to 1. The observation noise consists of both background noise and impulsive noise. The background noise has covariance matrix ${\mathit{C}}_{\eta }=c_{\eta }^2\mathit{I}$ where $c_{\eta }^2=2\times {10}^{-4}$. The impulsive noise follows the Bernoulli–Gaussian (BG) process as described in Section 8.1, where the probability $p_{r,i}$ is set to 0.05 and the constant $\kappa$ associated with the noise level is set to $\kappa =1000$. The step size $\mu$, the sparsity parameter $\lambda$, and the hard threshold value $\gamma$ are set to 0.9, 0.056, and 0.05, respectively. The normalized mean-square deviation (NMSD), i.e., $\mathrm{NMSD}\left[n\right]=\left|\right|\mathit{s}\left[n\right]-{\mathit{s}}_0{\left|\right|}^2/\left|\right|{\mathit{s}}_0{\left|\right|}^2$, is adopted as the measurement to evaluate the performance of the proposed algorithm.

Figure 14 shows the performance of the WL-G in comparison with that of the LMS algorithm with adaptive sampling. Similar to other simulation examples, the proposed WL-G algorithm significantly outperforms the LMS algorithm. The LMS algorithms produces an unreliable NMSD performance. In contrast, the WL-G algorithm is capable of effectively tracking the time-varying scenarios. Specifically, the sudden increases in the NMSD performance of the WL-G algorithm correspond to the changes in the spectral content of the signal. However, by adapting the sampling set, the WL-G algorithm is able to quickly converge to the steady state conditions. Figure 15 depicts the $\mathcal{F}$ curves of WL-G with adaptive sampling in comparison to LMS with adaptive sampling.

Figure 16, Figure 17 and Figure 18 report the samples chosen by the proposed WL-G algorithm at iterations $n=86$, $n=150$, and $n=279$.

9. Discussion

9.1. Discussion about Adaptive WL-G Estimation on Graphs

Comparing the Equation (9) with the LMS algorithm in [72], we obtain that the update equation of the WL-G algorithm contains just an extra scaling factor $exp(-\frac{({\mathit{y}}_{\mathit{i}}\left[\mathit{n}\right]-{\mathit{x}}_{\mathit{i}}{\left[\mathit{n}\right])}^{\mathbf{2}}}{2c^2})$ The outlier rejection property of theWelsch loss can be reflected from this factor [45]. Therefore, if the desired graph signal has impulsive characteristics or strong outliers, the WL-G is more stable than LMS algorithm in [72]. It is worth noting that the computational complexity of the WL-G algorithm is slightly more than that of the LMS algorithm. To implement the WL-G algorithm compared with LMS, we need only a few extra multiplications and additions.

9.1.1. Discussion about WL-G Estimation with AGS

The WL-G estimation on graphs in Table 5 supposes complement awareness of the support. Nevertheless, this assumption is fantastic in the practice, because the signal, the signal model, and the graph topology might be time-variant. The WL-G estimation with AGS is proposed. The signal support is tracked and estimated via the WL-G estimation, which is also fit to the AGS strategy.

9.1.2. Discussion about Simulation Results

Figure 3 and Figure 4 verify the precision of the mean square analytical theory provided in Theorem 1. From Figure 3 and Figure 4, we can determine that notional outcomes computed by (19) can fit fully with the simulated outcomes. Figure 6 verifies advantages belongs to WL-G algorithm on graph under the impulsive noise environment. From Figure 6, we get the WL-G can effectively cope with the impulsive noise by producing a much more reliable MSD performance. Figure 14 verifies the merits of WL-G algorithm with adaptive graph sampling. The proposed algorithm can get lower MSD and faster convergence rate. In summary, Figure 3 and Figure 4 in simulation section validate the precision of theory analysis in Theorem 1 while Figure 6 and Figure 14 verifies the advantages of proposed algorithms.

10. Conclusions

We propose the WL-G algorithm for adaptive graph signal reconstruction with impulsive noise in this work. Different to existing LMS methods, which are based on the least-squares criterion, the proposed WL-G leverages the use of Welsch loss to formulate its cost function. The framework of the WL-G algorithm is given in Table 5. The proposed algorithm can leverage the graph signal’s basal framework to recover signals from partial observations with impulsive noises under a band-limited supposition. The relationship between the sampling strategy and the performance of the WL-G algorithm was revealed via the theoretical analysis on the performance of the WL-G algorithm. Effective sampling strategies were then developed relied on analysis. An adaptive graph sampling technique was developed and used in conjunction with the WL-G algorithm to determine the support in the GFD while allowing the graph sampling strategy to be adapted in an online manner. Extensive simulation is carried out to verify the analysis and merits of the proposed algorithms. Specifically, Figure 3 and Figure 4 show that our theoretical analysis is correct, since theoretical results match well with the simulated results. Figure 6 show that the WL-G is sturdy against the impulsive interference. Figure 14 shows the proposed WL-G algorithm with the adaptive sampling strategy.