Vector-Valued Graph Trend Filtering with Non-Convex Penalties

TL;DR

This paper introduces a non-convex regularization approach for denoising vector-valued graph signals, demonstrating improved recovery performance and providing theoretical error bounds and an efficient algorithm.

Contribution

It extends graph trend filtering to vector-valued signals with non-convex penalties, offering better denoising and recovery capabilities along with convergence guarantees.

Findings

Non-convex regularizers outperform convex ones in denoising tasks.

The proposed ADMM algorithm converges reliably.

Numerical experiments validate improved performance on real-world data.

Abstract

This work studies the denoising of piecewise smooth graph signals that exhibit inhomogeneous levels of smoothness over a graph, where the value at each node can be vector-valued. We extend the graph trend filtering framework to denoising vector-valued graph signals with a family of non-convex regularizers, which exhibit superior recovery performance over existing convex regularizers. Using an oracle inequality, we establish the statistical error rates of first-order stationary points of the proposed non-convex method for generic graphs. Furthermore, we present an ADMM-based algorithm to solve the proposed method and establish its convergence. Numerical experiments are conducted on both synthetic and real-world data for denoising, support recovery, event detection, and semi-supervised classification.

Tables

Table 1. TABLE I : Key notations used in this paper.

Symbol	Description	Dimension
$𝚫$	oriented incidence matrix	$m\times n$
${𝚫}^{\left(k+1\right)}$	$k$th order graph difference operator	$r\times n$
$𝜷$	scalar-valued graph signal	$n$
$𝑩$	vector-valued graph signal	$n\times d$
$𝒚$	noisy observation of $𝜷$	$n$
$𝒀$	noisy observation of $𝑩$	$n\times d$
${𝚫}_{\ell \cdot }$	$\ell$-th row of $𝚫$	$n$
${𝑩}_{\cdot j}$	$j$-th column of $𝑩$	$n$
$\Vert {𝚫}^{\left(k+1\right)}\Vert$	spectral norm of ${𝚫}^{\left(k+1\right)}$	$1$

Table 2. TABLE II : Time complexity analysis of Alg.  1 .

	$k\ge 1$	$k=0$
${𝑫}^{\top }𝑫$ eigen decomposition	$O\left(r{n}^2+n^3\right)$	$O\left(m^2+n^3\right)$
$𝒁$ initialization	$O\left(rnd\right)$	$O\left(md\right)$
$𝑿$ initialization	$O\left(n^2\right)$	$O\left(n^2\right)$
$𝑩$ update	$O\left(d\left(nr+n^2\right)\right)$	$O\left(d\left(m+n^2\right)\right)$
$𝑫𝑩$ calculation	$O\left(rnd\right)$	$O\left(md\right)$
$𝒁,𝑼$ update	$O\left(rd\right)$	$O\left(rd\right)$
Total after $t$ iterations	$O\left(tdrn+td{n}^2\right)$	$O\left(tdm+td{n}^2\right)$

Table 3. TABLE III : Noisy input and reconstructed signal SNRs for eight measurements of varying input SNRs, rounded to two significant figures. Highest reconstructed signal SNR for each measurement is in bold .

	Average	#1	#2	#3	#4	#5	#6	#7	#8
Input SNR (dB)	8.7	-14	0	0	3.5	5.8	12	29	34
Vector-GTF + ${\ell }_1$	29	10	20	23	26	36	37	39	38
Scalar-GTF + ${\ell }_1$	21	0	11	13	16	18	26	41	45
Vector-GTF + SCAD	32	10	20	22	25	36	35	49	61
Scalar-GTF + SCAD	29	0	15	17	25	35	34	47	60
Vector-GTF + MCP	32	10	20	22	25	36	35	49	61
Scalar-GTF + MCP	29	0	15	22	24	30	33	49	60

Table 4. TABLE IV : Misclassification rates averaged over 10 trials, with p-values from running sampled t-tests between SCAD/MCP misclassification rates and the corresponding rates using ℓ 1 subscript ℓ 1 \ell_{1} . Cases where non-convex penalties perform better than ℓ 1 subscript ℓ 1 \ell_{1} with p-value below 0.1 0.1 0.1 are highlighted in bold , and where they perform worse are in italic .

		Heart	Wine quality	Wine	Iris	Breast	Car
# of samples ($n$)		303	1599	178	150	569	1728
# of classes ($K$)		2	6	3	3	2	4
$k=0$	${\ell }_1$	0.148	0.346	0.038	0.036	0.042	0.172
	SCAD p-value	0.148	0.353	0.038	0.033	0.042	0.149
		1.	0.06	1.	0.27	1.	0.06
	MCP p-value	0.144	0.351	0.037	0.035	0.040	0.148
		0.23	0.18	0.34	0.34	0.35	0.05
$k=1$	${\ell }_1$	0.143	0.351	0.034	0.039	0.035	0.104
	SCAD p-value	0.144	0.350	0.034	0.039	0.035	0.104
		0.30	0.43	0.34	1.	0.71	0.66
	MCP p-value	0.146	0.350	0.034	0.039	0.034	0.103
		0.05	0.44	0.34	1.	0.02	0.23

Equations

$$\mathbf{\beta }={\left[\begin{array}{c}{\beta }_1,{\beta }_2,\dots ,{\beta }_n\end{array}\right]}^{\top }\in {\mathbb{R}}^n,$$

$$\Delta_{i\ell}=\left\{\begin{array}[]{rl}-1,&\ell=j;\\ 1,&\ell=k;\\ 0,&\mbox{otherwise}\end{array}\right..$$

$${\mathbf{\Delta }}^{(k+1)}=\{\begin{array}{ll}{\mathbf{\Delta }}^{(1)\top }{\mathbf{\Delta }}^{(k)}\in {\mathbb{R}}^{n\times n},& \text{ odd k }\\ {\mathbf{\Delta }}^{(1)}{\mathbf{\Delta }}^{(k)}\in {\mathbb{R}}^{m\times n},& \text{ even k }\end{array}.$$

$$\mathbf{y}={\mathbf{\beta }}^{\star }+\mathbf{ϵ},\mathbf{ϵ}\sim \mathcal{N}(0,{\sigma }^2\mathbf{I}),$$

$$\underset{\mathbf{\beta }\in {\mathbb{R}}^n}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{\beta }{\Vert }_2^2+\lambda \Vert \mathbf{\Delta }\mathbf{\beta }{\Vert }_2^2,$$

$$\underset{\mathbf{\beta }\in {\mathbb{R}}^n}{\min }\frac{1}{2}\Vert \mathbf{y}-\mathbf{\beta }{\Vert }_2^2+\lambda \Vert {\mathbf{\Delta }}^{(k+1)}\mathbf{\beta }{\Vert }_1,$$

$$f(\mathbf{\beta })=\frac{1}{2}\Vert \mathbf{y}-\mathbf{\beta }{\Vert }_2^2+g({\mathbf{\Delta }}^{(k+1)}\mathbf{\beta };\lambda ,\gamma ),\mathbf{\beta }\in {\mathbb{R}}^n,$$

$$g({\mathbf{\Delta }}^{(k+1)}\mathbf{\beta })\triangleq g({\mathbf{\Delta }}^{(k+1)}\mathbf{\beta };\lambda ,\gamma )=\sum _{\ell =1}^r\rho (({\mathbf{\Delta }}^{(k+1)}\mathbf{\beta }{)}_{\ell };\lambda ,\gamma )$$

$${\rho }_{\mathrm{SCAD}}(t;\lambda ,\gamma )$$

$${\rho }_{\mathrm{MCP}}(t;\lambda ,\gamma )$$

$$\mathbf{Y}={\mathbf{B}}^{\star }+\mathbf{E},$$

$$\underset{\mathbf{B}\in {\mathbb{R}}^{n\times d}}{\min }\sum _{j=1}^{d}f({\mathbf{B}}_{\cdot j}).$$

$$\underset{\mathbf{B}\in {\mathbb{R}}^{n\times d}}{\min }\frac{1}{2}\Vert \mathbf{Y}-\mathbf{B}{\Vert }_{\mathrm{F}}^2+h({\mathbf{\Delta }}^{(k+1)}\mathbf{B};\lambda ,\gamma ),$$

$$h({\mathbf{\Delta }}^{(k+1)}\mathbf{B};\lambda ,\gamma )=\sum _{\ell =1}^r\rho \left(\Vert ({\mathbf{\Delta }}^{(k+1)}\mathbf{B}{)}_{\ell \cdot }{\Vert }_2;\lambda ,\gamma \right).$$

$$0\in {\nabla }_{\mathbf{\beta }}f(\mathbf{\beta }){|}_{\mathbf{\beta }=\hat{\mathbf{\beta }}}.$$

$$\kappa_{T,d}({\bm{\Delta}^{(k+1)}})=\inf_{{\bm{B}}\in{\mathbb{R}}^{n\times d}}\Bigg{\{}\frac{\sqrt{|T|}\cdot\|{\bm{B}}\|_{\mathrm{F}}}{\sum_{\ell\in T}\|({\bm{\Delta}^{(k+1)}}{\bm{B}})_{\ell\cdot}\|_{2}}\Bigg{\}}.$$

$$\displaystyle\frac{\|{\widehat{{\bm{\beta}}}}-{{\bm{\beta}}^{\star}}\|_{2}^{2}}{n}\leq\inf_{{\bm{\beta}}\in{\mathbb{R}}^{n}}\Bigg{\{}\frac{\|{\bm{\beta}}-{{\bm{\beta}}^{\star}}\|_{2}^{2}+4g({({\bm{\Delta}^{(k+1)}}{\bm{\beta}})_{T^{c}}})}{n}\Bigg{\}}$$

$$+\frac{2{\sigma }^2\left[C_G+2\sqrt{2C_G\log (\frac{1}{\delta })}+\frac{8{\zeta }_k^2|T|}{{\kappa }_T^2}\log (\frac{er}{\delta })\right]}{n(1-\mu \Vert {\mathbf{\Delta }}^{(k+1)}{\Vert }^2)}$$

$$\displaystyle\frac{\|\bm{\widehat{B}}-{\bm{B}}^{\star}\|_{\mathrm{F}}^{2}}{dn}\leq\inf_{{\bm{B}}\in{\mathbb{R}}^{n\times d}}\Bigg{\{}\frac{\|{\bm{B}}-{\bm{B}}^{\star}\|_{\mathrm{F}}^{2}+4h({({\bm{\Delta}^{(k+1)}}{\bm{B}})_{T^{c}}})}{dn}\Bigg{\}}$$

$$+\frac{2{\sigma }^2\left[C_G+2\sqrt{2C_G\log (\frac{d}{\delta })}+\frac{8{\zeta }_k^2|T|}{{\kappa }_T^2}\log (\frac{edr}{\delta })\right]}{n(1-\mu \Vert {\mathbf{\Delta }}^{(k+1)}{\Vert }^2)}$$

$$\frac{\Vert \hat{\mathbf{\beta }}-{\mathbf{\beta }}^{\star }{\Vert }_2^2}{n}$$

$$+\frac{2{\sigma }^2\left[C_G^{\delta }+8{\zeta }_k^2{\kappa }_T^{-2}|T|\log (\frac{er}{\delta })\right]}{n(1-\mu \Vert {\mathbf{\Delta }}^{(k+1)}{\Vert }^2)},$$

$$\frac{\Vert \hat{\mathbf{\beta }}-{\mathbf{\beta }}^{\star }{\Vert }_2^2}{n}$$

$$\frac{\Vert \hat{\mathbf{\beta }}-{\mathbf{\beta }}^{\star }{\Vert }_2^2}{n}=O\left(\frac{\lambda \Vert {\mathbf{\Delta }}^{(k+1)}{\mathbf{\beta }}^{\star }{\Vert }_1}{n}+\frac{{\sigma }^2{C}_G}{n}\right).$$

$${\zeta }_k\le {\lambda }_{\min }({\mathbf{\Delta }}^{(2)}{)}^{-\frac{k+1}{2}},$$

$${\kappa }_T({\mathbf{\Delta }}^{(k+1)})\ge (2d_{\max }{)}^{-\frac{k+1}{2}}.$$

$$\frac{\Vert \hat{\mathbf{\beta }}-{\mathbf{\beta }}^{\star }{\Vert }_2^2}{n}$$

$$+\min \left\{\frac{g({\mathbf{\Delta }}^{(k+1)}{\mathbf{\beta }}^{\star })}{n},\frac{{\sigma }^2\Vert {\mathbf{\Delta }}^{(k+1)}{\mathbf{\beta }}^{\star }{\Vert }_0\log n}{n}\right\},$$

$$S_k({\mathbf{\beta }}^{\star })=\left\{t\in \{1,\cdots ,r\}:({\mathbf{\Delta }}^{(k+1)}{\mathbf{\beta }}^{\star }{)}_t\ne 0\right\}.$$

$$d_G(S_k(\hat{\mathbf{\beta }})|S_k({\mathbf{\beta }}^{\star }))$$

$$=\{\begin{array}{ll}\underset{t_1\in S_k({\mathbf{\beta }}^{\star })}{\max }\underset{t_2\in S_k(\hat{\mathbf{\beta }})}{\min }{d}_G(t_1,t_2),& \text{if}{S}_k({\mathbf{\beta }}^{\star })\ne \varnothing \\ \infty & \text{otherwise}\end{array}.$$

$$d_G(S_k(\hat{\mathbf{\beta }})|S_k({\mathbf{\beta }}^{\star }))=\{\begin{array}{ll}O\left(\frac{R_n}{H_r^2}\right),& k=0\\ O\left(\frac{R_n^{1/3}}{H_r^{2/3}}\right),& k=1\end{array},$$

Full text

I Introduction

Signal estimation from noisy observations is a classic problem in signal processing and has applications in signal inpainting, collaborative filtering, recommendation systems and other large-scale data completion problems. Since noise can have deleterious, cascading effects in many downstream tasks, being able to efficiently and accurately filter and reconstruct a signal is of significant importance.

With the explosive growth of information and communication, signals are generated at an unprecedented rate from various sources, including social networks, citation networks, biological networks, and physical infrastructure [2]. Unlike time series or images, these signals lie on complex, irregular structures, and require novel processing techniques, leading to the emerging field of graph signal processing [3]-[5]. This framework models the structure by a graph and generalizes concepts and tools from classical discrete signal processing to graph signal processing. The associated graph-structured data are referred to as graph signals.

In graph signal processing, a common assumption is that the graph signal is smooth with respect to the graph, that is, the signal coefficients do not vary much over local neighborhoods of the graph. However, this characterization is insufficient for many real-world signals that exhibit spatially inhomogeneous levels of smoothness over the graph. In social networks for example, within a given community or social circle, users’ profiles tend to be homogeneous, while within a different social circle they will be of different, yet still have homogeneous values. Consequently, the signal is often characterized by large variations between regions and small variations within regions such that there are localized discontinuities and patterns in the signal. As a result, it is necessary to develop representations and algorithms to process and analyze such piecewise smooth graph signals.

In this work, we study the denoising of the class of piecewise smooth graph signals (including but not limited to piecewise constant graph signals), which is complementary to the class of smooth graph signals that exhibit homogeneous levels of smoothness over the graph. The reconstruction of smooth graph signals has been well studied in previous work both within graph signal processing [4]-[10] as well as in the context of Laplacian regularization [11, 12].

The Graph Trend Filtering (GTF) framework [13], which applies total variation denoising to graph signals [14], is a particularly flexible and attractive approach that regularizes discrete graph differences using the $\ell_{1}$ norm. Although the $\ell_{1}$ norm based regularization has many attractive properties [15], the resulting estimates are biased toward zero for large coefficients. To alleviate this bias effect, non-convex penalties such as the Smoothly Clipped Absolute Deviation (SCAD) penalty [16] and the Minimax Concave Penalty (MCP) [17] have been proposed as alternatives. These penalties behave similarly to the $\ell_{1}$ norm when the signal coefficients are small, but tend to a constant when the signal coefficients are large. Notably, they possess the so-called oracle property: in the asymptotics of large dimension, they perform as well as the case where we know in advance the support of the sparse vectors [18]-[22].

In this work, we strengthen the GTF framework in [13] by considering a large family of possibly non-convex regularizers, including SCAD and MCP that exhibit superior reconstruction performance over $\ell_{1}$ minimization for the denoising of piecewise smooth graph signals. Furthermore, we extend the GTF framework to allow vector-valued signals, e.g. time series [23], on each node of the graph, which greatly broadens the applicability of GTF to applications in social networks [24], gene networks, and semi-supervised classification [25]-[27]. Through theoretical analyses and empirical performance, we demonstrate that the use of non-convex penalties improves the performance of GTF in terms of both reduced reconstruction error and improved support recovery, i.e. how accurately we can localize the discontinuities of the piecewise smooth signals. Our contributions can be summarized as follows:

• •

Theoretically, we derive the statistical error rates of the signal estimates, defined as first-order stationary points of the proposed GTF estimator. We derive the rates in terms of the noise level and the alignment of the ground truth signal with respect to the underlying graph, without making assumptions on the piecewise smoothness of the ground truth signal. The better the alignment, the more accurate the estimates. Importantly, the estimators do not need to be the global minima of the proposed non-convex problem, which are much milder requirements and important for the success of optimization. For denoising vector-valued signals, the GTF estimate is more accurate when each dimension of the signal shares similar patterns across the graph.

• •

Algorithmically, we propose an ADMM-based algorithm that is guaranteed to converge to a critical point of the proposed GTF estimator.

• •

Empirically, we demonstrate the performance improvements of the proposed GTF estimators with non-convex penalties on both synthetic and real data for signal estimation, support recovery, event detection, and semi-supervised classification.

The rest of this paper is organized as follows. Section II reviews related works and their relationships to the current paper. In Section III, we provide some background and definitions on graph signal processing and GTF. Section IV presents the proposed GTF framework with non-convex penalties and vector-valued graph signals. Section V develops its performance guarantees, and Section VI presents an efficient algorithm based on ADMM. Numerical performance of the proposed approach is examined on both synthetic and real-world data for denoising and semi-supervised classification in Section VII. Finally, we conclude in Section VIII and briefly discuss future work.

Throughout this paper, we use boldface letters $\bm{a}$ and $\bm{A}$ to represent vectors and matrices respectively. The transpose of $\bm{A}$ is denoted as $\bm{A}^{\top}$. The $\ell$-th row of a matrix $\bm{A}$ is denoted as $\bm{A}_{\ell\cdot}$, and the $j$-th column of a matrix $\bm{A}$ is denoted as $\bm{A}_{\cdot j}$. The cardinality of a set $T$ is denoted as $|T|$. For any set $T\subseteq\{1,2,...,r\}$ and ${\bm{x}}\in{\mathbb{R}}^{r}$, we denote $({\bm{x}})_{T}\in{\mathbb{R}}^{|T|}$ such that $x_{\ell}\in({\bm{x}})_{T}$ if and only if $\ell\in T$ for $\ell\in\{1,2,...,r\}$. Similarly, we define a submatrix $\bm{A}_{T\cdot}\in\mathbb{R}^{|T|\times d}$ of $\bm{A}\in\mathbb{R}^{r\times d}$ that corresponds to pulling out the rows of $\bm{A}$ indexed by $T$. The $\ell_{2}$ norm of a vector $\bm{a}$ is defined as $\|\bm{a}\|_{2}$, and the spectral norm of a matrix $\bm{A}$ is defined as $\|\bm{A}\|$. The pseudo-inverse of a matrix $\bm{A}$ is defined as $\bm{A}^{{\dagger}}$. For a function $h({\bm{x}}):{\mathbb{R}}^{p}\to{\mathbb{R}}$, we write $\nabla_{\bm{x}}h({\bm{x}})|_{{\bm{x}}={\bm{x}}^{*}}$ to denote the gradient or subdifferential of $h({\bm{x}})$, if they exist, evaluated at ${\bm{x}}={\bm{x}}^{*}$. When the intention is clear, this may be written concisely as $\nabla h({\bm{x}}^{*})$. We also follow the standard asymptotic notations. If for some constants $C,N>0$, $|f(n)|\leq C|g(n)|$ for all $n\geq N$, then $f(n)=O(g(n))$; if $g(n)=O(f(n))$, then $f(n)=\Omega(g(n))$. Finally, Table I summarizes some key notations used in this paper for convenience.

II Related Work and Connections

Estimators that adapt to spatial inhomogeneities have been well studied in the literature via regularized regression, total variation and splines [28]-[30]. Most of these methods involve locating change points or knots that denote a distinct change in the behavior of the function or the signal.

Our work is most related to the spatially adaptive GTF estimator introduced in [13] that smoothens or filters noisy signals to promote piecewise smooth behavior with respect to the underlying graph structure; see also [31]. In the same spirit as [28], the fused LASSO and univariate trend filtering framework developed in [14, 32, 33] use discrete difference operators to fit a time series signal using piecewise polynomials. The GTF framework generalizes univariate trend filtering by generalizing a path graph to arbitrarily complex graphs. Specifically, by appropriately defining the discrete difference operator, we can enforce piecewise constant, piecewise linear, and more generally piecewise polynomial behaviors over the graph structure. In comparison to previous work [13], in this paper, we have significantly expanded its scope by allowing vector-valued data over the graph nodes and a broader family of possibly non-convex penalties.

We note that while a significant portion of the relevant literature on GTF or the fused LASSO has focused on the sparsistency or support recovery conditions under which we can ensure the recovery of the location of the discontinuities or knots [34, 35], in this work, we study the asymptotic error rates of our estimator with respect to the mean squared error. Our analysis of error rates leverages techniques in [36, 37] that result in sharp error rates of total variation denoising via oracle inequalities, which we have carefully adapted to allow non-convex regularizers. The obtained error rates can be translated into bounds on support recovery or how well we can localize the boundary by leveraging techniques in [38].

Employing a graph-based regularizer that promotes similarities between the signal values at connected nodes has been investigated by many communities, such as graph signal processing, machine learning, applied mathematics, and network science. The Network LASSO proposed in [24], which is similar to the GTF framework with multi-dimensional or vector-valued data, focused on the development of efficient algorithms without any theoretical guarantees. The recent works by Jung et al. [26, 27, 39] have analyzed the performance of Network LASSO for semi-supervised learning when the graph signal is assumed to be clustered according to the labels using the network null space property and the network compatibility condition inspired by related concepts in compressed sensing [40]. In contrast, our analysis does not make assumptions on the graph signal, and the error rate is adaptive to the alignment of the signal and the graph structure used in denoising.

A well-studied generalization of the sparse linear inverse problem is when there are multiple measurement vectors (MMV), and the solutions are assumed to have a common sparsity pattern [41]-[43]. Sharing information across measurements, and thereby exploiting the conformity of the sparsity pattern, has been shown to significantly improve the performance of sparse recovery in compressive sensing and sparse coding [44]-[48]. Motivated by these works, we consider vector-valued graph signals that are regarded as multiple measurements of scalar-valued graph signals sharing discontinuity patterns.

There are a few variants of non-convex penalties that promote sparsity such as SCAD, MCP, weakly convex penalties, and $\ell_{q}$ ($0\leq q<1$) minimization [18, 49]-[52]. In this paper, we consider and develop theory for a family of non-convex penalties parametrized similarly to that in [18, 50] with SCAD and MCP as our prime examples, although it is valid for other non-convex penalties.

III Graph Signal Processing, Piecewise Smooth Signals, and Graph Trend Filtering

We consider an undirected graph $G=(\mathcal{V},\mathcal{E},\bm{A})$, where $\mathcal{V}=\{v_{1},\ldots,v_{n}\}$ is the set of nodes, $\mathcal{E}=\{e_{1},\ldots,e_{m}\}$ is the set of edges, and $\bm{A}=[A_{jk}]\in{\mathbb{R}}^{n\times n}$ is the unweighted adjacency matrix – also known as the graph shift operator [3]. The edge set $\mathcal{E}$ represents the connections of the undirected graph $G$, and the positive edge weight $A_{jk}$ measures the underlying relation between the $j$th and the $k$th node, such as a similarity, a dependency, or a communication pattern. Let a scalar-valued graph signal be defined as

[TABLE]

where $\beta_{i}$ denotes the signal coefficient at the $i$th node.

Let $\bm{\Delta}\in{\mathbb{R}}^{m\times n}$ be the oriented incidence matrix of $G$, where each row corresponds to an edge. That is, if the edge $e_{i}=(j,k)\in\mathcal{E}$ connects the $j$th node to the $k$th node ($j<k$), the entries in the $i$th row of $\bm{\Delta}$ is then given as

[TABLE]

The entries of the signal $\bm{\Delta}\bm{\beta}=[(\beta_{k}-\beta_{j})]_{(j,k)\in\mathcal{E}}$ specifies the unweighted pairwise differences of the graph signal over each edge. As a result, $\bm{\Delta}$ can be interpreted as a graph difference operator. In graph signal processing, a signal is called smooth over a graph $G$ if $\|\bm{\Delta}\bm{\beta}\|_{2}^{2}=\sum_{(j,k)\in\mathcal{E}}(\beta_{k}-\beta_{j})^{2}$ is small.

III-A Piecewise Smooth Graph Signals

In practice, the graph signal may not be necessarily smooth over the entire graph, but only locally within different pieces of the graph. To model inhomogeneous levels of smoothness over a graph, we say that a graph signal $\bm{\beta}$ is piecewise constant over a graph $G$ if many of the differences $\beta_{k}-\beta_{j}$ are zero for $(j,k)\in\mathcal{E}$. Consequently, the difference signal $\bm{\Delta}\bm{\beta}$ is sparse and $\|\bm{\Delta}\bm{\beta}\|_{0}$ is small.

We can characterize piecewise $k$th order polynomial signals on a graph, where the piecewise constant case corresponds to $k=0$, by generalizing the notion of graph difference operators. Specifically, we use the following recursive definition of the $k$th order graph difference operator ${\bm{\Delta}^{(k+1)}}$ [13]. Let $\bm{\Delta}^{(1)}=\bm{\Delta}$ for $k=0$. For $k\geq 1$, let

[TABLE]

The signal $\bm{\beta}$ is said to be a piecewise $k$th order polynomial graph signal if $\|{\bm{\Delta}^{(k+1)}}\bm{\beta}\|_{0}$ is small. To further illustrate, let us consider the piecewise linear graph signal, corresponding to $k=1$, as a signal whose value at a node can be linearly interpolated from the weighted average of the values at neighboring nodes. It is easy to see that this is the same as requiring the second-order differences $\bm{\Delta}^{\top}\bm{\Delta}\bm{\beta}$ to be sparse. Similarly, we say that a signal has a piecewise quadratic structure over a graph if the differences between the second-order differences defined for piecewise linear signals are mostly zero, that is, if $\bm{\Delta}\bm{\Delta}^{\top}\bm{\Delta}\bm{\beta}$ is sparse. Fig. 1 illustrates various orders of piecewise graph smooth signals over the Minnesota road network graph.

III-B Denoising Piecewise Smooth Graph Signals via GTF

Assume we observe a noisy signal ${\bm{y}}$ over the graph under i.i.d Gaussian noise:

[TABLE]

and seek to reconstruct ${\bm{\beta}}^{\star}$ from ${\bm{y}}$ by leveraging the graph structure. When $\bm{\beta}$ is a smooth graph signal, Laplacian smoothing [11, 12, 53]-[55] can be used, which solves the following problem:

[TABLE]

where $\lambda>0$. However, it cannot localize abrupt changes in the graph signal when the signal is piecewise smooth.

Graph trend filtering (GTF) [13] is a flexible framework for estimation on graphs that is adaptive to inhomogeneity in the level of smoothness of an observed signal across nodes. The $k$th order GTF estimate is defined as:

[TABLE]

which can be regarded as applying total variation or fused LASSO with the graph difference operator ${\bm{\Delta}^{(k+1)}}$ [14, 56]. The sparsity-promoting properties of the $\ell_{1}$ norm have been well-studied [57]. Consequently, applying the $\ell_{1}$ penalty in GTF sets many of the (higher-order) graph differences to zero while keeping a small fraction of non-zero values. GTF is then adaptive over the graph; its estimate at a node adapts to the smoothness in its localized neighborhood.

IV Vector-Valued GTF with Non-convex Penalties

In this section, we first extend GTF to allow a broader family of non-convex penalties, and then extend it to handle vector-valued signals over the graph.

IV-A (Non-)convex Penalties

The $\ell_{1}$ norm penalty considered in (4) is well-known to produce biased estimates [58], which motivates us to extend the GTF framework to a broader class of sparsity-promoting regularizers that are not necessarily convex. We wish to minimize the following generalized $k$th order GTF loss function:

[TABLE]

where

[TABLE]

is a regularizer defined as the sum of the penalty function $\rho(\cdot;\lambda,\gamma):{\mathbb{R}}\to{\mathbb{R}}$ applied element-wise to ${\bm{\Delta}^{(k+1)}}{\bm{\beta}}$. Here, $r=m$ for even $k$ and $r=n$ for odd $k$ to account for different dimensions of ${\bm{\Delta}^{(k+1)}}$; see (1). We will refer to the GTF estimator that minimizes $f(\bm{\beta})$ as scalar-GTF.

Similarly to [18, 20, 50], we consider a family of penalty functions $\rho(\cdot;\lambda,\gamma)$ that satisfies the following assumptions.

Assumption 1**.**

Assume $\rho(\cdot;\lambda,\gamma)$ satisfies the following:

1. (a)

$\rho(t;\lambda,\gamma)$* satisfies $\rho(0;\lambda,\gamma)=0$, is symmetric around [math], and is non-decreasing on the real non-negative line.* 2. (b)

For $t\geq 0$, the function $t\mapsto\frac{\rho(t;\lambda,\gamma)}{t}$ is non-increasing in $t$. Also, $\rho(t;\lambda,\gamma)$ is differentiable for all $t\neq 0$ and sub-differentiable at $t=0$, with $\lim_{t\to 0^{+}}\rho^{\prime}(t;\lambda,\gamma)=\lambda$. This upper bounds $\rho(t;\lambda,\gamma)\leq\lambda|t|$. 3. (c)

There exists $\mu>0$ such that $\rho(t;\lambda,\gamma)+\frac{\mu}{2}t^{2}$ is convex.

Many penalty functions satisfy these assumptions. Besides the $\ell_{1}$ penalty, the non-convex SCAD [16] penalty

[TABLE]

and the MCP [17]

[TABLE]

also satisfy them. We note that Assumption 1 (c) is satisfied for SCAD with $\mu\geq\mu_{\mathrm{SCAD}}=\frac{1}{\gamma-1}$ and for MCP with $\mu\geq\mu_{\mathrm{MCP}}=\frac{1}{\gamma}$. Fig. 2 illustrates the $\ell_{1}$, SCAD and MCP penalties for comparison. While the non-convexity means that in general, we may not always find the global optimum of $f({\bm{\beta}})$, it often affords us many other advantages. SCAD and MCP both taper off to a constant value, and hence apply less shrinkage for higher values. As a result, they mitigate the bias effect while promoting sparsity. Further, they are smooth and differentiable for $t\geq 0$ and are both upper bounded by the $\ell_{1}$ penalty for all $t$.

IV-B Vector-Valued GTF

In many applications, the signals on each node are in fact multi-dimensional or vector-valued, e.g. time series in social networks, multi-class labels in semi-supervised learning, feature vectors of different objects in feature selection. Therefore, it is natural to consider an extension to the graph signal denoising problem, where the graph signal on each node is a $d$-dimensional vector instead of a scalar. In this scenario, we define a vector-valued graph signal to be piecewise smooth if it is piecewise smooth in each of its $d$ dimensions, and assume their discontinuities to coincide over the same small set of edges or nodes. Further, we denote the vector-valued signal of interest as ${\bm{B}}^{\star}\in{\mathbb{R}}^{n\times d}$, such that the $i$th row of the matrix ${\bm{B}}$ corresponds to the $i$th node of the graph. The noise model for the observation matrix $\bm{Y}\in{\mathbb{R}}^{n\times d}$ is defined as

[TABLE]

where each element of $\bm{E}\in{\mathbb{R}}^{n\times d}$ is drawn i.i.d from $\mathcal{N}(0,\sigma^{2})$. A naïve approach is to estimate each column ${\bm{B}}_{\cdot j}$ of ${\bm{B}}$ separately via scalar-GTF:

[TABLE]

However, this formulation does not take full advantage of the multi-dimensionality of the graph signal. Instead, when the columns of ${\bm{B}}$ are correlated, coupling them can be beneficial such that we encourage the sharing of information across dimensions or features. For example, if one column ${\bm{B}}_{\cdot i}$ exhibits strong piecewise smoothness over the graph, and therefore has compelling evidence about the relationship between nodes, sharing that information to a related column ${\bm{B}}_{\cdot j}$ can improve the overall denoising and filtering performance. As a result, we formulate a vector-GTF problem as follows:

[TABLE]

where the new penalty function $h({\bm{\Delta}^{(k+1)}}{\bm{B}})\triangleq h({\bm{\Delta}^{(k+1)}}{\bm{B}};\lambda,\gamma):{\mathbb{R}}^{r\times d}\to{\mathbb{R}}$ is the sum of $\rho(\cdot;\lambda,\gamma)$ applied to the $\ell_{2}$ norm of each row of ${\bm{\Delta}^{(k+1)}}{\bm{B}}\in{\mathbb{R}}^{r\times d}$:

[TABLE]

By enforcing sparsity on $\left\{\|(\bm{\Delta}^{(k+1)}{\bm{B}})_{\ell\cdot}\|_{2}\right\}_{1\leq l\leq r}$, we are coupling $\bm{\Delta}^{(k+1)}{{\bm{B}}}_{\cdot j}$ to be of similar sparsity patterns across $j=1,\ldots,d$. Note the difference from (9), where elements of $(\bm{\Delta}^{(k+1)}{\bm{B}})_{\ell\cdot}$ can be set to zero or non-zero independently.

V Theoretical Guarantees

In this section, we present the error rates and support recovery guarantees of the generalized GTF estimators, namely scalar-GTF (5) and vector-GTF (10), under the AWGN noise model. Before continuing, we first define a few useful quantities. Let $C_{G}$ be the number of connected components in the graph $G$, or equivalently, the dimension of the null space of ${\bm{\Delta}^{(k+1)}}$. Further, let $r$ be the number of rows of ${\bm{\Delta}^{(k+1)}}$, and $\zeta_{k}$ be the maximum $\ell_{2}$ norm of the columns of ${\bm{\Delta}^{(k+1)\dagger}}$.

V-A Error Rates of First-order Stationary Points

Due to non-convexity, global minima of the proposed GTF estimators may not be attainable. Therefore, it is more desirable to understand the statistical performance of any first-order stationary points of the GTF estimators. We call ${\widehat{{\bm{\beta}}}}\in{\mathbb{R}}^{n}$ a stationary point of $f({\bm{\beta}})$, if it satisfies

[TABLE]

We further introduce the compatibility factor, which generalizes the notion used in [36] to allow vector-valued signals.

Definition 1** (Compatibility factor).**

Let ${\bm{\Delta}^{(k+1)}}$ be fixed. The compatibility factor $\kappa_{T,d}$ of a set $T\subseteq\{1,2,\ldots,r\}$ is defined as $\kappa_{\varnothing,d}=1$, and for nonempty set $T$,

[TABLE]

To further build intuition, consider $\sqrt{|T|}\kappa_{T,1}(\bm{\Delta})^{-1}=\sup_{{\bm{\beta}}\in{\mathbb{R}}}\{\|(\bm{\Delta})_{T}{\bm{\beta}}\|_{1}/\|{\bm{\beta}}\|_{2}\}$. This is precisely the definition of $\|(\bm{\Delta})_{T}\|_{1,2}$, an induced norm of the $|T|\times n$ submatrix of $\bm{\Delta}$. If we consider signals with fixed power $\|{\bm{\beta}}\|_{2}^{2}=1$, $\|(\bm{\Delta})_{T}{\bm{\beta}}\|_{1}$ will depend on how much the $T$ edges are connected to each other. Together with $\|(\bm{\Delta})_{T}{\bm{\beta}}\|_{1}\leq\sqrt{|T|}\|(\bm{\Delta})_{T}{\bm{\beta}}\|_{2}$, $\kappa_{T,1}(\bm{\Delta})$ can be related to the restricted eigenvalue condition, which is often used to bound the performance of LASSO [40]. With slight abuse of notation, we write $\kappa_{T}:=\kappa_{T,d}$.

We have the following oracle inequality that is applicable to the stationary points of the GTF estimators. The proof in Appendix A follows a construction that is similar to Theorem 2 in [36]. The oracle inequality holds for any ${\widehat{{\bm{\beta}}}}$ that satisfies the first-order optimality condition, allowing the use of non-convex penalties. This mild condition on ${\widehat{{\bm{\beta}}}}$ is a key difference from [13, Theorem 3] and [1, Theorem 1] that are applicable to the global minimizer, which is difficult to guarantee when using non-convex penalties. We also stress that although GTF was motivated by piecewise smooth graph signals, Theorem 1 holds for any graph $G$ and graph signal ${{\bm{\beta}}^{\star}}$.

Theorem 1** (Oracle inequality of GTF stationary points).**

Assume $\mu<1/\|{\bm{\Delta}^{(k+1)}}\|^{2}$. Fix $\delta\in(0,1)$. For scalar-GTF (5), let ${\widehat{{\bm{\beta}}}}$ be a stationary point. Set $\lambda=\sigma\zeta_{k}\sqrt{2\log\left(\frac{er}{\delta}\right)}$, then

[TABLE]

with probability at least $1-2\delta$ for any $T\subseteq\{1,2,...,r\}$. Similarly, for vector-GTF (10), let $\bm{\widehat{B}}$ be a stationary point. Set $\lambda=\sigma\zeta_{k}\sqrt{2d\log(\frac{edr}{\delta})}$, then

[TABLE]

with probability at least $1-2\delta$ for any $T\subseteq\{1,2,...,r\}$.

Remark 1**.**

Recall that $\mu$ is defined in Assumption 1 (c), which characterizes how “non-convex” the regularizer is, and dictates the inflection point in Fig. 2. The assumption $\mu<1/\|{\bm{\Delta}^{(k+1)}}\|^{2}$ in Theorem 1 therefore implicitly constrains the level of non-convexity of the regularizer. Take MCP in (7) for example: since $\mu\geq 1/\gamma$, we can guarantee the existence of a valid $\mu$ such that $\mu<1/\|{\bm{\Delta}^{(k+1)}}\|^{2}$ as long as we set $\gamma>\|{\bm{\Delta}^{(k+1)}}\|^{2}$.

Theorem 1 allows one to select ${\bm{\beta}}$ and $T$ to optimize the error bounds on the right hand side of (12) and (13). For example, pick ${\bm{\beta}}={{\bm{\beta}}^{\star}}$ in (12) (hence an “oracle”) to have

[TABLE]

where $C_{G}^{\delta}=C_{G}+2\sqrt{2C_{G}\log(\frac{1}{\delta})}$.

• •

By setting $T$ as an empty set, we have

[TABLE]

which suggest that the reconstruction accuracy improves when the ground truth $\bm{\beta}^{\star}$ is better aligned with the graph structure, and consequently the value of $g({\bm{\Delta}^{(k+1)}}{\bm{\beta}}^{\star})$ is small.

• •

On the other hand, by setting $T$ as the support of ${\bm{\Delta}^{(k+1)}}{{\bm{\beta}}^{\star}}$, we achieve

[TABLE]

which grows linearly as we increase the sparsity level $\|{\bm{\Delta}^{(k+1)}}{{\bm{\beta}}^{\star}}\|_{0}$.

Similar discussions can be conducted for vector-GTF by choosing ${\bm{B}}={\bm{B}}^{\star}$ in (13). More importantly, we can directly compare the performance of vector-GTF with scalar-GTF, which was formulated for vector-valued graph signals in (9). The error bound of vector-GTF pays a small price in the order of $\log d$, but is tighter than scalar-GTF if $h(({\bm{\Delta}^{(k+1)}}{\bm{B}}^{\star})_{T^{c}})\ll\sum_{j=1}^{d}g(({\bm{\Delta}^{(k+1)}}{\bm{B}}^{\star}_{\cdot j})_{T^{c}})$. This suggests that vector-GTF is much more advantageous when the support sets of ${\bm{\Delta}^{(k+1)}}{\bm{B}}^{\star}_{\cdot j}$ for $j=1,\ldots,d$ overlap, i.e. when the local discontinuities and patterns in ${\bm{B}}^{\star}_{\cdot j}$ are shared.

V-B Comparison with Scalar-GTF using $\ell_{1}$ Regularization

We compare our error bound for scalar-GTF that is on $\|{\widehat{{\bm{\beta}}}}-{{\bm{\beta}}^{\star}}\|_{2}^{2}/n$ with [13, Theorem 3], which is obtained for GTF with the $\ell_{1}$ penalty, reproduced below for convenience.

Theorem 2** (Basic error bound of $\ell_{1}$ GTF minimizer).**

If $\lambda=\Theta(\sigma\zeta_{k}\sqrt{\log r})$, then ${\widehat{{\bm{\beta}}}}$, the minimizer of (4), satisfies

[TABLE]

The above bound is comparable to our bound in the special case of setting $T$ to an empty set, i.e. (15). The first term of the bound in (15) is upper bounded by that of Theorem 2. The non-convex regularization yields especially tighter bounds when ${\bm{\Delta}^{(k+1)}}{\bm{\beta}}^{\star}$ contains large coefficients, so that $g({\bm{\Delta}^{(k+1)}}{\bm{\beta}}^{\star})\ll\lambda\|{\bm{\Delta}^{(k+1)}}{\bm{\beta}}^{\star}\|_{1}$. On the other hand, the second term of (15) contains $1-\mu\|{\bm{\Delta}^{(k+1)}}\|^{2}$ in the denominator, which makes it an upper bound of the second term in Theorem 2. This gap can be brought down by choosing a larger $\gamma$, which allows one to pick a smaller $\mu$, as mentioned in Remark 1. However, as $\gamma\to\infty$, non-convex SCAD and MCP also tends to $\ell_{1}$, which erases the improvement from using non-convex regularizers in the first term of the bound. This indicates a trade-off in the overall error bound based on $\gamma$, or the “non-convexity” of the regularizers chosen for scalar-GTF.

To sum up, despite being non-convex, we can guarantee that any stationary point of the proposed GTF estimator possesses strong statistical guarantees.

V-C Error Rates for Erdős-Rényi Graphs

We next specialize Theorem 1 to the Erdős-Rényi random graphs using spectral graph theory [59]. Let $d_{\max}$ and $d_{0}$ respectively be the maximum and expected degree of the graph. It is known that for any graph it holds [13]

[TABLE]

where $\lambda_{\min}(\bm{\Delta}^{(2)})$ is the smallest non-zero eigenvalue of the graph Laplacian matrix $\bm{\Delta}^{(2)}$. Moreover, we have $\|{\bm{\Delta}^{(k+1)}}\|^{2}=(\lambda_{\max}(\bm{\Delta}^{(2)}))^{k+1}$, and $d_{\max}+1\leq\lambda_{\max}(\bm{\Delta}^{(2)})\leq 2d_{\max}$ [60]. Next, we present a simple lower bound on $\kappa_{T}$, which is proved in Appendix -B.

Proposition 1** (Bound on $\kappa_{T}$).**

$\kappa_{T}$* is bounded for any $T$ and $d$ as*

[TABLE]

For an Erdős-Rényi random graph, if $d_{0}=\Omega\left(\log(n)\right)$, we have $d_{\max}=O(d_{0})$ almost surely [61, Corollary 8.2] and $C_{G}=1$. Furthermore, $\lambda_{\min}(\bm{\Delta}^{(2)})=\Omega(d_{0}-\sqrt{d_{0}})$ [13, 59, 62], and $r=n$ for odd $k$ and $r=O(nd_{0})$ for even $k$. Therefore, with probability at least $1-n^{-10}$, we have

[TABLE]

where $g({{\bm{\Delta}^{(k+1)}}{{\bm{\beta}}^{\star}}})\lesssim\frac{\sigma\|{\bm{\Delta}^{(k+1)}}{{\bm{\beta}}^{\star}}\|_{1}\sqrt{\log n}}{d_{0}^{(k+1)/2}}$ by plugging in $g({{\bm{\Delta}^{(k+1)}}{{\bm{\beta}}^{\star}}})\leq\lambda\|{\bm{\Delta}^{(k+1)}}{{\bm{\beta}}^{\star}}\|_{1}$.

These results are also applicable to $d_{0}$-regular Ramanujan graphs [62].

V-D Support Recovery

An alternative yet important metric for gauging the success of the proposed GTF estimators is support recovery, which aims to localize the discontinuities in the piecewise smooth graph signals, i.e. the support set of ${\bm{\Delta}^{(k+1)}}{\bm{\beta}}^{\star}$, that is

[TABLE]

In particular, for odd $k$, the discontinuities correspond to graph nodes; and for even $k$, they correspond to the edges. Let ${\widehat{{\bm{\beta}}}}$ be the GTF estimate of the graph signal. The quality of the support recovery can be measured using the graph screening distance [38]. For any $t_{1}\in S_{k}({\bm{\beta}}^{\star})$ and $t_{2}\in S_{k}({\widehat{{\bm{\beta}}}})$, let $d_{G}(t_{1},t_{2})$ denote the length of the shortest path between them. The distance of $S_{k}({\widehat{{\bm{\beta}}}})$ from $S_{k}({\bm{\beta}}^{\star})$ is then defined as

[TABLE]

Interestingly, Lin et.al. [38] showed recently that under mild assumptions, one can translate the error bound into a support recovery guarantee. Specifically, letting $R_{n}$ be the RHS of (12) that bounds the error $\|{\widehat{{\bm{\beta}}}}-{{\bm{\beta}}^{\star}}\|_{2}^{2}/n$ in Theorem 1, we have

[TABLE]

where $H_{r}$ quantifies the minimum level of discontinuity, defined as the minimum absolute value of the non-zero values of ${\bm{\Delta}^{(k+1)}}{\bm{\beta}}^{\star}$, i.e.

[TABLE]

Consequently, this leads to support recovery guarantees of the proposed GTF estimators. Numerical experiment in Section VII-A verifies the superior performance of the non-convex regularizers over the $\ell_{1}$ regularizer for support recovery.

VI ADMM Algorithm and its Convergence

There are many algorithmic approaches to optimize the vector-GTF formulation in (10), since scalar-GTF (5) can be regarded as a special case with $d=1$. In this section, we illustrate the approach adopted in this paper, which is the Alternating Direction Method of Multipliers (ADMM) framework for solving separable optimization problems [63].

Via a change of variable as $\bm{Z}={\bm{\Delta}^{(k+1)}}{\bm{B}}$, we can transform (10) to

[TABLE]

Its corresponding Lagrangian can be written as:

[TABLE]

where $\bm{U}\in{\mathbb{R}}^{r\times d}$ is the Lagrangian multiplier, and $\tau$ is the parameter. Alg. 1 shows the ADMM updates based on the Lagrangian in (21). Recall the proximal operator is defined as $\texttt{Prox}_{f}(\bm{v};\alpha)=\operatorname{argmin}_{\bm{x}}\frac{1}{2}\|{\bm{x}}-\bm{v}\|_{2}^{2}+\alpha f({\bm{x}})$ for a function $f(\cdot)$. $\ell_{1}$, SCAD and MCP all admit closed-form solutions of Prox, which are simple thresholding operations [64]. Furthermore, we have the following convergence guarantee for Alg. 1, whose proof is provided in Appendix -C. Theorem 3 implies that the output of Alg. 1 satisfies Theorem 1.

Theorem 3**.**

Let $\tau\geq\mu$, then Alg. 1 converges to a stationary point of (10).

In addition, we provide a detailed time complexity analysis of Alg. 1 in Table II. Note that since $\bm{\Delta}$ is a sparse matrix with exactly $2m$ non-zero entries, Alg. 1 can run much faster when $k=0$. As a preprocessing step for each ${\bm{D}}$, we compute $\bm{V}\in{\mathbb{R}}^{n\times n}$ and $\bm{S}\in{\mathbb{R}}^{n\times n}$, the eigenvectors and eigenvalues of ${\bm{D}}^{\top}{\bm{D}}$, exactly once. $\bm{X}=\bm{V}(\bm{I}+\tau\bm{S})^{-1}\bm{V}^{\top}$ can then be initialized very efficiently for all experiments that use ${\bm{D}}$.

VII Numerical Experiments

For the following experiments, we fixed $\gamma=3.7$ for SCAD, $\gamma=1.4$ for MCP. Further, the graphs we use in the following experiments satisfy Assumption 2 for this choice of $\gamma$. Unless explicitly mentioned, we tuned $\lambda$ and $\frac{\tau}{\lambda}$ for each experiment using the Hyperopt toolbox [65]. To meet the convergence criteria in Theorem 3, we enforced $\tau\geq 1/\gamma$. SCAD/MCP were warm-started with the GTF estimate with $\ell_{1}$ penalty. Python packages PyGSP [66] and NetworkX [67] were used to construct and plot graphs. The input signal SNR was calculated as $10\log_{10}(\|{\bm{B}}^{\star}\|_{\mathrm{F}}/\sigma^{2}nd)$, while the reconstructed signal SNR was calculated as $10\log_{10}(\|{\bm{B}}^{\star}\|_{\mathrm{F}}/\|\bm{\widehat{B}}-{\bm{B}}^{\star}\|_{\mathrm{F}})$, where $\bm{\widehat{B}}$ was the reconstruction. Computation time was measured with MacBook Pro 2017 with an 2.9 GHz Intel Core i7 and 16GB RAM. Our code is available at https://github.com/HarlinLee/nonconvex-GTF-public.

VII-A Denoising via GTF with Non-convex Regularizers

We first highlight via synthetic examples two important advantages that non-convex regularizers provide over the $\ell_{1}$ penalty.

• •

Bias Reduction: We demonstrate the reduction in signal bias in Fig. 3 for the graph signal defined over a $12\times 12$ 2D-grid graph, using both the $\ell_{1}$ penalty and the MCP penalty. Clearly, the MCP estimate (orange) has less bias than the $\ell_{1}$ estimate (blue), and can recover the ground truth surface (purple) more closely.

• •

Support Recovery: We illustrate the improved support recovery performance of non-convex penalties on localizing the boundaries for a piecewise constant signal on the Minnesota road graph, shown in Fig. 5. Particularly, we look at how well our estimator localizes the support of ${\bm{\Delta}^{(k+1)}}{{\bm{\beta}}^{\star}}$, that is, the discontinuity of the piecewise constant graph signal by looking at how well we can classify an edge as connecting two nodes in the same piece or being a cut edge across two pieces. By sweeping the regularization parameter $\lambda$, we obtain the ROC curve in Fig. 4, i.e. the true positive rate versus the false positive rate of classifying a cut edge correctly, and see that scalar-GTF with MCP and SCAD consistently outperforms the scalar-GTF with $\ell_{1}$ penalty.

Then, we compare the performance of GTF using non-convex regularizers such as SCAD and MCP with that using the $\ell_{1}$ norm more rigorously. For the ground truth signal ${{\bm{\beta}}^{\star}}$, we construct a piecewise constant signal on a $20\times 20$ 2D-grid graph and the Minnesota road graph [66] as shown in the left panel of Fig. 5, and add different levels of noise following (2). We recover the signal by scalar-GTF with Alg. 1, and plot the SNR of the reconstructed signal versus the SNR of the input signal in solid lines in the middle panel of Fig. 5, averaged over $10$ and $20$ realizations, respectively. SCAD/MCP consistently outperforms $\ell_{1}$ in denoising graph signals defined over both regular and irregular structures.

VII-B Denoising Vector-valued Signals via GTF

We compare the performance of vector-GTF in (10) with (9), which applies scalar-GTF to each column of the vector-valued graph signal. The convex $\ell_{1}$ norm, and the non-convex SCAD and MCP are employed. We reuse the same ground truth graph signals over the 2D-grid graph and the Minnesota road graph constructed in Section VII-A in Fig. 5. $d$ independent noisy realizations of the graph signal are concatenated to construct a noisy vector-valued graph signal with dimension $d=10$ on the 2D-grid graph and with $d=20$ on the Minnesota road graph. We recover the vector-valued signal by minimizing vector-GTF (10) with Alg. 1.

The middle panel of Fig. 5 plots the average SNR of the reconstructed signal versus the average SNR of the input signal in dotted lines. We emphasize that the performance of (9) is the same as applying scalar-GTF to each realization, which is shown in the middle panel of Fig. 5 in solid lines. As before, SCAD/MCP consistently outperforms $\ell_{1}$ in denoising signals over both regular and irregular graphs. Furthermore, as expected, due to the sharing of information across realizations, vector-GTF consistently outperforms scalar-GTF, especially in the low SNR regime.

The right panel of Fig. 5 plots the computation time versus the gain in SNR from denoising via vector-GTF. 10 trials are performed for each regularizer with the input signal SNR fixed at 20dB. Parameter tuning and eigen decomposition of $\bm{\Delta}^{(2)}$ are preprocessing steps, and hence they are not included in the time measurement; but for reference, the eigen decomposition took 0.025 and 2.5 seconds for 2D-grid and Minnesota graphs, respectively. Since GTF with non-convex regularizers are warm-started by the $\ell_{1}$ estimate, the runtime for $\ell_{1}$ GTF is added to the SCAD/MCP runtime. Overall, running vector-GTF with SCAD/MCP after once with $\ell_{1}$ takes more time, but with large benefits in the denoising performance. Even with the additional computation time, Vector-GTF runs reasonably fast; with the Minnesota road network, where $n=2642$ and $m=3304$, computation takes less than 25 seconds.

We further investigate the benefit of sharing information across measurements or realizations in the following experiment, using the same ground truth signal on the 2D-grid graph. We stack eight noisy realizations of this same piecewise constant signal to build a vector-valued signal. We construct these noisy measurements by scaling each one of them differently and randomly such that each will have SNR $\sim\log_{10}\mbox{Uniform}[-10,30]$dB under (2). This has the effect of rendering some measurements more informative than others, and potentially allowing vector-GTF to reap the benefits of sharing information across measurements. We recover the 8-dimensional graph signal via Alg. 1 using $\ell_{1}$, SCAD, and MCP regularizers, and in Table III, report the input signal and reconstructed signal SNRs for each measurement in addition to the average SNRs. $\lambda$ is fixed at $0.5\sigma^{2}$.

First of all, notice that as before, using SCAD/MCP generally achieves results with higher SNR than using $\ell_{1}$, and that on average, minimizing (10) outperforms minimizing (9). The effect of sharing information across measurements is most apparent in low SNR settings, when information about the boundaries of the graph signal can be borrowed from higher SNR signals to improve the estimation. On the other hand, sharing information with noisier signals does not help denoising signals with high input SNR. However, it is worth noting that, unlike $\ell_{1}$, SCAD/MCP does not see decrease in its performance in the high SNR settings.

VII-C Event Detection with NYC Taxi Data

To further illustrate graph trend filtering on a real-world dataset, we consider the road network of Manhattan where the nodes correspond to junctions [68]. We map the pickups and dropoffs of the NYC taxi trip dataset [69, 70] to the nearest road junctions, and define the total count at that junction to be the signal value on the corresponding graph node. The signal of interest, plotted on the top left panel of Fig. 6, is the difference between the event graph signal on the day of NYC Gay Pride parade, 12-2pm on June 26, 2011, and the seasonal average graph signal at the same time during the 8 nearest Sundays. During the event, no pickups and dropoffs could occur in the areas shown in the top right panel of Fig. 6 [71]. We denoise the signal via GTF using both $\ell_{1}$ and MCP, where we chose $\lambda$ such that $\|\bm{\Delta}{\widehat{{\bm{\beta}}}}\|_{0}\approx 200$. Once again, we observe the GTF estimate with MCP produces sharper traces around the parade route, indicating better capabilities of event detection and localization.

VII-D Semi-supervised Classification

Graph-based learning provides a flexible and attractive way to model data in semi-supervised classification problems when vast amounts of unlabeled data are available compared to labeled data, and labels are expensive to acquire [11, 12, 55]. One can construct a nearest-neighbor graph based on the similarities between each pair of samples, and hope to propagate the label information from labeled samples to unlabeled ones. We move beyond our original problem in (5) to a $K$-class classification problem in a semi-supervised learning setting, where for a given dataset with $n$ samples, we observe a subset of the one-hot encoded class labels, $\bm{Y}\in{\mathbb{R}}^{n\times K}$, such that $Y_{ij}=1$ if $i$th sample has been observed to be in $j$th class, and $Y_{ij}=0$ otherwise. A diagonal indicator matrix $\bm{M}\in{\mathbb{R}}^{n\times n}$ denotes samples whose class labels have been observed. Then, we can define the modified absorption problem [12, 13, 55] using a variation of GTF to estimate the unknown class probabilities ${\bm{B}}\in\mathbb{R}^{n\times K}$:

[TABLE]

where $\bm{R}\in{\mathbb{R}}^{n\times K}$ (set to be uniform in the experiment) is a fixed prior belief, and $\epsilon>0$ determines how much emphasis to be given to the prior belief. The labels $\widetilde{\bm{Y}}$ can be estimated using $\widetilde{{\bm{B}}}$ such that $\widetilde{Y}_{ij}=1$ if and only if $j=\arg\max_{1\leq\ell\leq K}\widetilde{B}_{i\ell}$, and otherwise $\widetilde{Y}_{ij}=0$. Note that this can be completely separated into $K$ scalar-GTF problems, one corresponding to each class.

We applied the algorithm in $\eqref{eq:ssl}$ to 6 popular UCI classification datasets [72] with $\epsilon=0.01$. For each dataset, we normalized each feature to have zero mean and unit variance, and constructed a 5-nearest-neighbor graph of the samples based on the Euclidean distance between their features, with edge weights from the Gaussian radial basis kernel. We observed the labels of 20% of samples in each class randomly. Table IV shows the misclassification rates averaged over $10$ repetitions, which demonstrates that the performance using non-convex penalties such as SCAD/MCP are at least competitive with, and often better than, those with the $\ell_{1}$ penalty.

VIII Conclusions

We presented a framework for denoising piecewise smooth signals on graphs that generalizes the graph trend filtering framework to handle vector-valued signals using a family of non-convex regularizers. We provided theoretical guarantees on the error rates of our framework, and derived a general ADMM-based algorithm to solve this generalized graph trend filtering problem. Furthermore, we demonstrated the superior performance of these non-convex regularizers in terms of reconstruction error, bias reduction, and support recovery on both synthetic and real-world data. In particular, its performance on semi-supervised classification is investigated. In future work, we plan to further study this approach when the graph signals are observed by indirect and incomplete measurements.

-A Proof of Theorem 1

Proof.

We denote ${\bm{D}}={\bm{\Delta}^{(k+1)}}$. Define $R$ as the row space of ${\bm{D}}$, and $R^{\bot}$ the null space. Let $\mathcal{P}_{\mathrm{R}}={\bm{D}^{\dagger}}{\bm{D}}$, the projection onto $R$, and $\|{\bm{x}}\|_{\mathrm{R}}=\|\mathcal{P}_{\mathrm{R}}{\bm{x}}\|_{2}$. Additionally, $\mathcal{P}_{\mathrm{R}^{\bot}}=\bm{I}-{\bm{D}^{\dagger}}{\bm{D}}$, the projection onto $R^{\bot}$. Since ${\widehat{{\bm{\beta}}}}$ is a stationary point of $f({\bm{\beta}})$, it follows that

[TABLE]

By the chain rule, $\nabla_{\bm{\beta}}g({{\bm{D}}{\bm{\beta}}})|_{{\bm{\beta}}={\widehat{{\bm{\beta}}}}}=\{{\bm{D}}^{\top}\bm{z}:\bm{z}\in\nabla_{\bm{x}}g({\bm{x}})|_{\bm{x}={\bm{D}}{\widehat{{\bm{\beta}}}}}\}$. Then by (23), there exists $\bm{z}\in\nabla_{\bm{x}}g({\bm{x}})|_{\bm{x}={\bm{D}}{\widehat{{\bm{\beta}}}}}$, such that

[TABLE]

In particular, $\forall{\bm{\beta}}\in{\mathbb{R}}^{n}$, we have

[TABLE]

and, specializing to ${\widehat{{\bm{\beta}}}}$,

[TABLE]

Subtract (25) from (24), and use the definition of subgradient to get $\forall{\bm{\beta}}\in{\mathbb{R}}^{n}$,

[TABLE]

By the measurement model ${\bm{y}}={{\bm{\beta}}^{\star}}+\bm{\epsilon}$ and the polarization equality, i.e. $2\bm{a}^{\top}\bm{b}=\|\bm{a}\|_{2}^{2}+\|\bm{b}\|_{2}^{2}-\|\bm{a}-\bm{b}\|_{2}^{2}$, the left-hand side of (-A) can be rewritten as

[TABLE]

Combining (-A) and (27) gives us $\forall{\bm{\beta}}\in{\mathbb{R}}^{n}$

[TABLE]

Let us first consider $\bm{\epsilon}^{\top}({\widehat{{\bm{\beta}}}}-{\bm{\beta}})$. From Hölder’s inequality,

[TABLE]

By standard tail bounds for independent Gaussian random variables, we have with probability at least $1-\delta$,

[TABLE]

Additionally, recognize that $\|\bm{\epsilon}\|_{\mathrm{R}^{\bot}}^{2}$ is a chi-squared random variable with $C_{G}$ degrees of freedom. We can then invoke the one-sided tail bound for chi-squared random variables (c.f. [73, Example 2.5]) such that for any $0\leq t\leq 1$,

[TABLE]

Consequently, with probability at least $1-\delta$,

[TABLE]

The inequalities (31) and (30) hold simultaneously with probability at least $1-2\delta$. Then, using $\lambda\|\bm{x}\|_{1}\leq g(\bm{x})+\frac{\mu}{2}\|\bm{x}\|_{2}^{2}$ and $\lambda=\sigma\zeta_{k}\sqrt{2\log(\frac{er}{\delta})}\geq\|({\bm{D}}^{\dagger})^{\top}\bm{\epsilon}\|_{\infty}$, we can bound (29) further as

[TABLE]

Together with $\|{\bm{D}}({\widehat{{\bm{\beta}}}}-{\bm{\beta}})\|_{2}^{2}\leq\|{\bm{D}}\|^{2}\|({\widehat{{\bm{\beta}}}}-{\bm{\beta}})\|_{2}^{2}$, we can upper bound (-A) as

[TABLE]

Note that for any set $T$, $g({{\bm{x}}})=\sum_{i\in T}\rho({x_{i}})+\sum_{j\in T^{c}}\rho({x_{j}})=g({({\bm{x}})_{T}})+g({({\bm{x}})_{T^{c}}})$. Therefore, using the triangle inequality and subadditivity and symmetry of $\rho$,

[TABLE]

We bound (33) further by the compatibility factor,

[TABLE]

Now combining (32), (33), and (34), we then have

[TABLE]

Apply Young’s inequality, which is $2ab\leq a^{2}/\epsilon+\epsilon b^{2}$ for $\epsilon>0$, with $a=\|\mathcal{P}_{\mathrm{R}^{\bot}}\bm{\epsilon}\|_{2}+2\lambda\sqrt{|T|}\kappa_{T}^{-1}$, $b=\|{\widehat{{\bm{\beta}}}}-{\bm{\beta}}\|_{2}$, and $\epsilon=1-\mu\|{\bm{D}}\|^{2}>0$, we have

[TABLE]

Therefore,

[TABLE]

Cancel $\|{\widehat{{\bm{\beta}}}}-{\bm{\beta}}\|_{2}^{2}$ on both sides, apply the infimum over ${\bm{\beta}}$ and plug in the bounds (31) to get

[TABLE]

The proof extends for the vector-GTF (10) in a similar manner. We need to replace (29) by

[TABLE]

where \|\mathcal{P}_{\mathrm{R}^{\bot}}\bm{E}\|_{\mathrm{F}}^{2}\leq d\sigma^{2}\Big{(}C_{G}+2\sqrt{2C_{G}\log(d/\delta)}\Big{)} with probability at least $1-\delta$. Similarly, for (34), we use the generalized definition of the compatibility factor $\kappa_{T}$, given as

[TABLE]

which will lead to the claimed bound in the theorem. ∎

-B Proof of Proposition 1

Proof.

By Cauchy-Schwartz inequality, we have

[TABLE]

and note that given two matrices $\bm{U}$ and $\bm{V}$, $(\bm{U}\bm{V})_{T}=(\bm{U})_{T}\bm{V}$ where $T$ is a subset of rows indices. We also use the fact that $\|\bm{U}\bm{V}\|_{\mathrm{F}}\leq\|\bm{U}\|\|\bm{V}\|_{\mathrm{F}}$. We consider two cases:

• •

For even $k$, we have

[TABLE]

Note that $(\bm{\Delta})_{T}$ is equivalent to the incidence matrix of a subgraph with only $T$ edges, which allows us to bound,

[TABLE]

where $d_{i}$ is the degree of node $i$.

• •

For odd $k$, we have

[TABLE]

where $\bm{\Delta}^{(2)}_{T\times T}\in{\mathbb{R}}^{|T|\times|T|}$ is the principal submatrix of $\bm{\Delta}^{(2)}$ indexed by $T$. By Cauchy’s interlacing theorem, the maximum eigenvalue of the submatrix is upper bounded, so

[TABLE]

Therefore, for all $k$, $\|({\bm{\Delta}^{(k+1)}}{\bm{B}})_{T}\|_{\mathrm{F}}\leq\sqrt{2d_{\max}}\|\bm{\Delta}^{(k)}{\bm{B}}\|_{\mathrm{F}}$. To conclude the proof, note

[TABLE]

∎

-C Proof of Theorem 3

We show the convergence of Alg. 1 by proving a modified version of [22, Proposition 1]. The superscript $(m)$ denotes the values of ${\bm{B}},{\bm{Z}},{\bm{U}}$ at the $m$th iteration of the loop inside Alg. 1.

Proposition 2** (Convergence to a feasible solution).**

If $\tau\geq\mu$, then the primal residual $r^{(m)}=\|{\bm{\Delta}^{(k+1)}}{\bm{B}}^{(m)}-{\bm{Z}}^{(m)}\|_{\mathrm{F}}$ and the dual residual $s^{(m+1)}=\|\tau({\bm{\Delta}^{(k+1)}})^{\top}({\bm{Z}}^{(m+1)}-{\bm{Z}}^{(m)})\|_{\mathrm{F}}$ of Alg. 1 satisfy that $\lim_{m\to\infty}r^{(m)}=0$ and $\lim_{m\to\infty}s^{(m)}=0$.

Proof.

Denote $\bm{D}={\bm{\Delta}^{(k+1)}}$, and $\rho_{\lambda}(\cdot)=\rho(\cdot;\lambda,\gamma)$. Recall from Assumption 1 (c) that there exists $\mu>0$ such that $\rho_{\lambda}(\|{\bm{x}}\|_{2})+\frac{\mu}{2}\|{\bm{x}}\|_{2}^{2}$ is convex. Now consider the Lagrangian $\mathcal{L}({\bm{B}},{\bm{Z}},{\bm{U}})$ with regard to the $\ell$-th row ${\bm{z}}_{\ell}$ of ${\bm{Z}}=[{\bm{z}}_{1}^{\top},...{\bm{z}}_{r}^{\top}]^{\top}$, assuming all other variables are fixed:

[TABLE]

where $\bm{c}_{1}$ and $c_{2}$ represent terms of $\mathcal{L}({\bm{B}},{\bm{Z}},{\bm{U}})$ that do not depend on ${\bm{z}}_{\ell}$. With our choice of $\tau\geq\mu$, then $\mathcal{L}({\bm{B}},{\bm{Z}},{\bm{U}})$ is convex with regard to each of ${\bm{B}}$, ${\bm{U}}$, and for each row of ${\bm{Z}}$, allowing us to apply [74, Theorem 5.1]. Therefore, Alg. 1 converges to limit points ${\bm{B}}^{\star},{\bm{Z}}^{\star},{\bm{U}}^{\star}$.

Then it follows that the dual residual $\lim_{m\to\infty}s^{(m)}=\|\tau{\bm{D}}^{\top}({\bm{Z}}^{\star}-{\bm{Z}}^{\star})\|_{\mathrm{F}}=0$. For the primal residual, notice that the ${\bm{U}}$ update step in line 10 of Alg. 1 also shows that $\forall m,t\geq 0$,

[TABLE]

Fixing $t$ and setting $m\to\infty$, we have

[TABLE]

holds $\forall t\geq 0$. Hence, ${\bm{D}}{\bm{B}}^{\star}-{\bm{Z}}^{\star}=\bm{0}$, and therefore $\lim_{m\to\infty}r^{(m)}=\|{\bm{D}}{\bm{B}}^{\star}-{\bm{Z}}^{\star}\|_{\mathrm{F}}=0$. ∎

This proposition shows that the algorithm in the limit achieves primal and dual feasibility, and that the Augmented Lagrangian in (21) with ${\bm{Z}}^{\star}$ and ${\bm{U}}^{\star}$ becomes the original GTF formulation in (10). ${\bm{B}}$ that is produced by every iteration of Alg. 1 is a stationary point of (21) with fixed ${\bm{Z}}$ and ${\bm{U}}$. As a result, ${\bm{B}}^{\star}$ is a stationary point of (10).