A sharp oracle inequality for Graph-Slope

TL;DR

This paper introduces Graph-Slope, a new estimator for signals on graphs, providing a sharp prediction error bound, an efficient algorithm, and experimental validation, advancing the analysis of graph-based denoising methods.

Contribution

It presents a sharp oracle inequality for Graph-Slope, improves results over existing Total Variation denoisers, and offers an efficient dual-based algorithm with experimental support.

Findings

Sharp oracle inequality for Graph-Slope prediction error

Improved performance over Total Variation denoiser

Efficient dual-based algorithm demonstrated in experiments

Abstract

Following recent success on the analysis of the Slope estimator, we provide a sharp oracle inequality in term of prediction error for Graph-Slope, a generalization of Slope to signals observed over a graph. In addition to improving upon best results obtained so far for the Total Variation denoiser (also referred to as Graph-Lasso or Generalized Lasso), we propose an efficient algorithm to compute Graph-Slope. The proposed algorithm is obtained by applying the forward-backward method to the dual formulation of the Graph-Slope optimization problem. We also provide experiments showing the interest of the method.

Equations

$$\text{pen }:\beta \to \lambda \Vert {\mathbf{D}}^{\top }\beta {\Vert }_1,$$

$$\Vert \theta {\Vert }_{[\lambda ]}=\sum _{j=1}^p{\lambda }_j|\theta {|}_j^{\downarrow },$$

$$\underset{b\in {\mathbb{R}}^p}{\min }\frac{1}{2n}\Vert y-X\beta {\Vert }^2+\Vert \beta {\Vert }_{[\lambda ]}.$$

$$({\mathbf{D}}^{\top }{)}_{e,v}=\{\begin{array}{ll}+1,& \text{if }v=\min (i,j)\\ -1,& \text{if }v=\max (i,j)\\ 0,& \text{otherwise}\end{array},$$

$$\Vert v{\Vert }^{\ast }=\underset{\Vert \beta \Vert ⩽1}{\sup }⟨v,\beta ⟩.$$

$$y={\beta }^{\star }+\epsilon ,$$

$$\hat{\beta }:={\hat{\beta }}^{\mathrm{GS}}\in \underset{\beta \in {\mathbb{R}}^p}{\mathrm{argmin}}\frac{1}{2n}\Vert y-\beta {\Vert }^2+\Vert {\mathbf{D}}^{\top }\beta {\Vert }_{[\lambda ]},$$

$$\Vert {\mathbf{D}}^{\top }\beta {\Vert }_{[\lambda ]}=\sum _{j=1}^p{\lambda }_j|{\mathbf{D}}^{\top }\beta {|}_j^{\downarrow },$$

$$\Vert \theta {\Vert }_{[\lambda ]}=\underset{\tau \in {\mathfrak{S}}_p}{\max }\sum _{j=1}^p{\lambda }_j|{\theta }_{\tau (j)}|=\sum _{j=1}^p{\lambda }_j|\theta {|}_j^{\downarrow }.$$

$${\hat{\beta }}^{\mathrm{GL}}\in \underset{\beta \in {\mathbb{R}}^p}{\mathrm{argmin}}\frac{1}{2n}\Vert y-\beta {\Vert }^2+{\lambda }_1\Vert {\mathbf{D}}^{\top }\beta {\Vert }_1,$$

$$\Lambda(\lambda,s)=\Big{(}\sum_{j=1}^{s}\lambda_{j}^{2}\Big{)}^{1/2}.$$

$$\frac{1}{n}\Vert {\mathbf{D}}^{†}\epsilon {\Vert }_{[\lambda ]}^{\ast }\le 1/2$$

$$\frac{1}{n}\Vert \hat{\beta }-{\beta }^{\star }{\Vert }^2\le \underset{s\in [p]}{\inf }\left[\underset{\begin{array}{c}\beta \in {\mathbb{R}}^n\\ \Vert {\mathbf{D}}^{\top }\beta {\Vert }_0\le s\end{array}}{\inf }\frac{1}{n}\Vert \beta -{\beta }^{\star }{\Vert }^2+\frac{1}{2n}{\left(\frac{3n\Lambda (\lambda ,s)}{2\kappa (s)}+\frac{\sigma +2\sigma \sqrt{2\log (1/\delta )}}{\sqrt{n}}\right)}^2\right],$$

$$\kappa (s)\triangleq \underset{v\in {\mathbb{R}}^n:3\Lambda (\lambda ,s)\Vert {\mathbf{D}}^{\top }v{\Vert }_2>{\sum }_{j=s+1}^p{\lambda }_j|{\mathbf{D}}^{\top }v{|}_j^{\downarrow }}{\inf }\left(\frac{\Vert v\Vert }{\Vert {\mathbf{D}}^{\top }v{\Vert }_2}\right).$$

$$f(e)=\underset{v\in {\mathbb{R}}^n:\Vert v\Vert =1}{\sup }\left[e^{\top }v+n\Lambda (\lambda ,s)\Vert {\mathbf{D}}^{\top }v{\Vert }_2-n\sum _{j=s+1}^p{\lambda }_j|{\mathbf{D}}^{\top }v{|}_j^{\downarrow }\right].$$

$$\frac{1}{2}(\Vert \hat{\beta }-$$

$$\le$$

$$=$$

$$\le$$

$${\epsilon }^{\top }v={\epsilon }^{\top }\Pi v+{\epsilon }^{\top }({\mathrm{Id}}_n-\Pi )v\le \Vert \Pi \epsilon {\Vert }_2+{\epsilon }^{\top }({\mathrm{Id}}_n-\Pi )v.$$

$$g(e)=\underset{v\in {\mathbb{R}}^n:\Vert v\Vert =1}{\sup }\left[e^{\top }({\mathrm{Id}}_n-\Pi )v+n\Lambda (\lambda ,s)\Vert {\mathbf{D}}^{\top }v{\Vert }_2-n\sum _{j=s+1}^p{\lambda }_j|{\mathbf{D}}^{\top }v{|}_j^{\downarrow }\right].$$

$$\mathbb{P}\left(g(\epsilon )\le \mathrm{Med}[g(\epsilon )]+\sigma \sqrt{2\log (1/\delta )}\right)\ge 1-\delta ,$$

$$\mathrm{Med}[g(\epsilon )]\le (3n/2)\Lambda (\lambda ,s)/\kappa (s).$$

$$\mathbb{P}(g(\epsilon )\le (3n/2)\Lambda (\lambda ,s)/\kappa (s))\ge 1/2.$$

$${\epsilon }^{\top }({\mathrm{Id}}_n-\Pi )v$$

$$\le \Vert (({\mathbf{D}}^{\top }{)}^{†}{)}^{\top }\epsilon {\Vert }_{[\lambda ]}^{\ast }\Vert {\mathbf{D}}^{\top }v{\Vert }_{[\lambda ]}=\Vert {\mathbf{D}}^{†}\epsilon {\Vert }_{[\lambda ]}^{\ast }\Vert {\mathbf{D}}^{\top }v{\Vert }_{[\lambda ]},$$

$$g(\epsilon )$$

$$\le \underset{v\in {\mathbb{R}}^n:\Vert v\Vert =1}{\sup }\left[\frac{n}{2}\sum _{j=1}^s{\lambda }_j|{\mathbf{D}}^{\top }v{|}_j^{\downarrow }+n\Lambda (\lambda ,s)\Vert {\mathbf{D}}^{\top }v{\Vert }_2-\frac{n}{2}\sum _{j=s+1}^p{\lambda }_j|{\mathbf{D}}^{\top }v{|}_j^{\downarrow }\right]$$

$$\le \frac{n}{2}\underset{v\in {\mathbb{R}}^n:\Vert v\Vert =1}{\sup }\left[3\Lambda (\lambda ,s)\Vert {\mathbf{D}}^{\top }v{\Vert }_2-\sum _{j=s+1}^p{\lambda }_j|{\mathbf{D}}^{\top }v{|}_j^{\downarrow }\right].$$

$$3\Lambda (\lambda ,s)\Vert {\mathbf{D}}^{\top }v{\Vert }_2-\sum _{j=s+1}^p{\lambda }_j|{\mathbf{D}}^{\top }v{|}_j^{\downarrow }\le 3\Lambda (\lambda ,s)\Vert v\Vert /\kappa (s)=3\Lambda (\lambda ,s)/\kappa (s).$$

$$3\Lambda (\lambda ,s)\Vert {\mathbf{D}}^{\top }v{\Vert }_2-\sum _{j=s+1}^p{\lambda }_j|{\mathbf{D}}^{\top }v{|}_j^{\downarrow }\le 0\le 3\Lambda (\lambda ,s)/\kappa (s).$$

$$g_j=e_j^{\top }{\mathbf{D}}^{†}\epsilon ,\text{ for all }j=1,\dots ,p.$$

$$(1/n)\Vert {\mathbf{D}}^{†}\epsilon {\Vert }_{[\lambda ]}^{\ast }\le \underset{j=1,\dots ,p}{\max }\left(|g{|}_j^{\downarrow }/(n{\lambda }_j)\right).$$

Full text

62G08,

62J07,

denoising,

graph signal regularization,

oracle inequality,

convex optimization,

and

t1J. Salmon and S. Vaiter were supported by the CNRS PE1 “ABS” grant.

1. 1 Introduction

2. 1.1 Model and notation

3. 1.2 The Graph-Slope estimator

4. 2 Theoretical guarantees: sharp oracle inequality

5. 3 Numerical experiments

6. 3.1 Algorithm for Graph-Slope

7. 3.2 Synthetic experiments

8. 3.3 Example on real data: Paris roads network

9. A Preliminary lemmas

1 Introduction

Many inference problems of interest involve signals defined on discrete graphs. This includes for instance two-dimensional imaging but also more advanced hyper-spectral imaging scenarios where the signal lives on a regular grid. Two types of structure arise naturally in such examples: The first type of structures comes from regularity or smoothness of the signal, which led to the development of wavelet methods. The second type of structure involves signals with few sharp discontinuities. For instance in one dimension, piecewise constant signals appear when transition states are present, the graph being a 1D path. In imaging, where the underlying graph is a regular 2D grid, occlusions create piece-wise smooth signals rather than smooth ones.

This paper studies regularizers for signals with sharp discontinuities. A popular choice in imaging is the Total Variation (TV) regularization [24]. For 1D signals, TV regularization has also long been used in statistics [18]. If an additional $\ell_{1}$ regularization is added, this is sometimes referred to as the fused Lasso [30, 32, 13].

A natural extension of such methods to arbitrary graphs relies on $\ell_{1}$ analysis penalties [15] which involve the incidence matrix of the underlying graph, see for instance [25] or the Edge Lasso of [26]. Such penalties have the form

[TABLE]

where $\lambda>0$ is a tuning parameter and ${\mathbf{D}^{\top}}$ is the (edge-vertex) incidence matrix of the graph defined below, and $\beta$ represents the signal to be recovered. This approach is notably different from contributions in machine learning where $\ell_{2}$ penalties, *i.e., *Laplacian regularization, have been considered for spectral clustering [27, 20] (see also [33] for a review). Theoretical results in favor of the $\ell_{1}$ norm instead of the squared $\ell_{2}$ norm are studied in [25].

Penalties based on $\ell_{0}$ regularization with the graph incidence matrix have recently been analyzed [16], including an analysis of their approximation algorithm. They are of interest as they do not suffer from the (shrinkage) bias created by the convex $\ell_{1}$ norm. However, such methods present the difficulty that in the general case they lead to non-convex problems. Note that the 1D path is an exception since the associated optimization problem can be solve using dynamic programming [1]. Concerning the bias reduction though, simpler remedies could be used, including least-squares refitting on the model space associated, applying for instance the CLEAR method [14].

Following the introduction of the Slope regularization in the context of high dimensional regression [10], we propose Graph-Slope, its generalization to contexts where the signal is supported on a graph. In linear regression, Slope [10] is defined as follows. Given $p$ tuning parameters $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{p}\geq 0$ with at least one strict inequality, define the ordered $\ell_{1}$ norm by

[TABLE]

where for any $\theta\in\mathbb{R}^{p}$, we use the notation111following the notation considered in [7] $(|\theta|_{1}^{\downarrow},\dots,|\theta|^{\downarrow}_{p})$ for the non-increasing rearrangement of its amplitudes $(|\theta_{1}|,\dots,|\theta_{p}|)$. Then, given a design matrix $X\in\mathbb{R}^{n\times p}$ and a response vector $y\in\mathbb{R}^{n}$, the Slope estimator is defined as a solution of the minimization problem

[TABLE]

If the parameters $\lambda_{1},\dots,\lambda_{p}$ are all equal, then Slope is equal to the Lasso with tuning parameter $\lambda_{1}$.

Slope presents several advantages compared to Lasso in sparse linear regression. First, Slope provably controls the False Discovery Rate (FDR) for orthogonal design matrices [10] and experiments show that this property is also satisfied for some non-orthogonal design matrices [10]. Second, it appears that Slope has more power than Lasso in the sense that Slope will discover more nonzero coefficients of the unknown target vector [10]. An interpretation of this phenomenon is that Lasso shrinks small coefficients too heavily and may thus miss the smallest nonzero coefficients of the target vector. On the other hand, the Slope penalty induces less shrinkage on small coefficients, leading to more power. Third, while Lasso with the universal parameter is known to achieve the rate of estimation of order $(s/n)\log(p)$ (where $s$ is the sparsity of the unknown target vector, $n$ the number of measurements and $p$ the number of covariates), Slope achieves the optimal rate of estimation of order $(s/n)\log(p/s)$ [28, 4].

We propose a theoretical and experimental analysis of Graph-Slope, the counterpart estimator of Slope for signals defined on graphs. Graph-Slope is defined in the next section. Our theoretical contribution for Graph-Slope borrows some technical details recently introduced in [17] to control the Mean Squared Error (MSE) for the Generalized Lasso.

Last but not least, we provide an efficient solver to compute the Graph-Slope estimator. It relies on accelerated proximal gradient descent to solve the dual formulation [3, 12, 22]. To obtain an efficient solver, we leverage the seminal contribution made in [36] showing the link between ordered $\ell_{1}$ norm (1.2) and isotonic regression. Hence, we can use fast implementations of the PAVA algorithm (for Pool Adjacent Violators Algorithm, see for instance [6]), available for instance in scikit-learn [23] for this purpose. Numerical experiments illustrate the benefit of Graph-Slope, in particular in terms of True Discovery Rate (TDR) performance.

A high level interpretation of our simulation results is as follows. In the model considered in this paper, a sharp discontinuity of the signal corresponds to an edge of the graph with nonzero coefficient. Since Graph-Lasso uses an $\ell_{1}$-penalty, the penalty level is uniform across all edges of the graph. Edges with small coefficients are too heavily penalized with Graph-Lasso. Using Graph-Slope lets us reduce the penalty level on the edges with small coefficients. This leads to the discovery of more discontinuities of the true signal as compared to Graph-Lasso.

1.1 Model and notation

Let $\mathcal{G}=(V,E)$ be an undirected and connected graph on $n$ vertices, $V=[n]$, and $p$ edges, $E=[p]$. This graph can be represented by its edge-vertex incidence matrix ${\mathbf{D}^{\top}}=\mathbf{D}_{\mathcal{G}}^{\top}\in\mathbb{R}^{p\times n}$ (we drop the reference to $\mathcal{G}$ when no ambiguity is possible) defined as

[TABLE]

where $e=\{i,j\}$. The matrix $\mathbf{L}=\mathbf{D}{\mathbf{D}^{\top}}$ is the so-called graph Laplacian of $\mathcal{G}$. The Laplacian $\mathbf{L}$ is invariant under a change of orientation of the graph.

For any $u\in\mathbb{R}^{p}$, we denote by $\|u\|_{0}$ the pseudo $\ell_{0}$ norm of $u$ : $\|u\|_{0}=|\{j\in[p]:u_{j}\neq 0\}|$, and for any matrix $\mathbf{A}$, we denote by $\mathbf{A}^{\dagger}$ its Moore-Penrose pseudo-inverse. The canonical basis of $\mathbb{R}^{p}$ is denoted $(e_{1},\dots,e_{p})$.

For any norm $\|\cdot\|$ on $\mathbb{R}^{n}$, the associated dual norm $\|\cdot\|^{*}$ reads at $v\in\mathbb{R}^{n}$

[TABLE]

As a consequence, for every $(\beta,v)\in\mathbb{R}^{n}\times\mathbb{R}^{n}$, one has $\langle\beta,\,v\rangle\leqslant\|\beta\|\|v\|^{*}$.

In this work, we consider the following denoising problem for a signal over a graph. Assume that each vertex $i\in[n]$ of the graph carries a signal $\beta^{\star}_{i}$. For each vertex $i\in[n]$ of the graph, one observes $y_{i}$, a noisy perturbation of $\beta^{\star}_{i}$. In vector form, one observes the vector $y\in\mathbb{R}^{n}$ and aims to estimate $\beta^{\star}\in\mathbb{R}^{n}$, i.e.,

[TABLE]

where $\varepsilon\sim\mathcal{N}(0,\sigma^{2}\mathrm{Id}_{n})$ is a noise vector. We will say that an edge $e=\{i,j\}$ of the graph carries the signal $({\mathbf{D}^{\top}}\beta^{\star})_{e}$. In particular, if two vertices $i$ and $j$ are neighbours and if they carry the same value of the signal, i.e., $\beta^{\star}_{i}=\beta^{\star}_{j}$, then the corresponding edge $e=\{i,j\}$ carries the constant signal. The focus of the present paper is on signals $\beta^{\star}$ that have few discontinuities. A signal $\beta^{\star}\in\mathbb{R}^{n}$ has few discontinuities if ${\mathbf{D}^{\top}}\beta^{\star}$ has few nonzero coefficients, i.e., $\|{\mathbf{D}^{\top}}\beta^{\star}\|_{0}$ is small, or equivalently if most edges of the graph carry the constant signal. In particular, if $\|{\mathbf{D}^{\top}}\beta^{\star}\|_{0}=s$, we say that $\beta^{\star}$ is a vector of ${\mathbf{D}^{\top}}$-sparsity $s$.

1.2 The Graph-Slope estimator

We consider in this paper the so-called Graph-Slope variational scheme:

[TABLE]

where

[TABLE]

with $\lambda=(\lambda_{1},\dots,\lambda_{p})\in\mathbb{R}^{p}$ satisfying $\lambda_{1}\geqslant\lambda_{2}\geqslant\cdots\geqslant\lambda_{p}\geq 0$, and using for any vector $\theta\in\mathbb{R}^{p}$ the notation222following the notation from [7]. $(|\theta|_{1}^{\downarrow},\dots,|\theta|^{\downarrow}_{p})$ for the non-increasing rearrangement of its amplitudes $(|\theta_{1}|,\dots,|\theta_{p}|)$. According to [10], $\|\cdot\|_{[\lambda]}$ is a norm over $\mathbb{R}^{p}$ if and only if $\lambda_{1}\geqslant\lambda_{2}\geqslant\cdots\geqslant\lambda_{p}\geq 0$ with at least one strict inequality. This is a consequence of the observation that if $\lambda_{1}\geqslant\lambda_{2}\geqslant\cdots\geqslant\lambda_{p}\geq 0$ then one can rewrite the Slope-norm of $\theta$ as the maximum over all $\tau\in\mathfrak{S}_{p}$ (the set of permutations over $[p]$), of the quantity $\sum_{i=1}^{p}\lambda_{j}|\theta_{\tau(j)}|$:

[TABLE]

The Generalized Lasso (also sometimes referred to as TV denoiser) relies on $\ell_{1}$ regularization. It was recently investigated in [17, 25], and can be defined as

[TABLE]

where $\|\cdot\|_{1}$ is the standard $\ell_{1}$ norm, and $\lambda_{1}>0$ is a tuning parameter.

If $\lambda_{1}=\lambda_{2}=\dots=\lambda_{p}$ then $\|\theta\|_{[\lambda]}=\lambda_{1}\|\theta\|_{1}$ for all $\theta\in\mathbb{R}^{p}$, so that the minimization problems (1.10) and (1.7) are the same. On the other hand, if $\lambda_{j}>\lambda_{j+1}$ for some $j=1,\dots,p-1$, then the optimization problems (1.10) and (1.7) differ. For instance, if $\lambda_{1}>\lambda_{2}>0$, all coefficients of ${\mathbf{D}^{\top}}\beta$ are equally penalized in the Graph-Lasso (1.10), while coefficients of ${\mathbf{D}^{\top}}\beta$ are not uniformly penalized in the Graph-Slope optimization problem (1.7). Indeed, in the Graph-Slope optimization problem (1.7), the largest coefficient of ${\mathbf{D}^{\top}}\beta$ is penalized as in (1.10) but smaller coefficients of ${\mathbf{D}^{\top}}\beta$ receive a smaller penalization. The Graph-Slope optimization problem (1.7) is more flexible than (1.10) as it allows the smaller coefficients of ${\mathbf{D}^{\top}}\beta$ to be less penalized than its larger coefficients. We will see in the next sections that this flexibility brings advantages to both the theoretical properties of $\hat{\beta}^{\rm GS}$ as well as its performance in simulations, as compared to $\hat{\beta}^{\rm GL}$.

2 Theoretical guarantees: sharp oracle inequality

We can now state the main theoretical result of the paper, a sharp oracle inequality for the Graph-Slope. For any integer $s$ and weights $\lambda=(\lambda_{1},\dots,\lambda_{p})$, define

[TABLE]

Theorem 2.1**.**

Assume that the Graph-Slope weights $\lambda_{1}\geq\dots\geq\lambda_{p}\geq 0$ are such that the event

[TABLE]

has probability at least $1/2$. Then, for any $\delta\in(0,1)$, we have with probability at least $1-2\delta$

[TABLE]

where $\Lambda(\cdot,\cdot)$ is defined in (2.1) and the compatibility factor $\kappa(s)$ is defined as

[TABLE]

Proof.

Let $\beta$ be a minimizer of the right hand side of (2.3) and let $s=\|{\mathbf{D}^{\top}}\beta\|_{0}$. Define the function $f(\cdot)$ by

[TABLE]

Let also $w=(\hat{\beta}-\beta)/(\|\hat{\beta}-\beta\|)$. By Lemma A.3 and Lemma A.1 with $\alpha=0$, for any $z\in\mathbb{R}^{p}$, we have

[TABLE]

where for the last inequality we used the elementary inequality $2ab\leq a^{2}+b^{2}$.

By Lemma A.4 we have $\mathrm{Id}_{n}=(\mathrm{Id}_{n}-\Pi)+\Pi$ with $\mathrm{Id}_{n}-\Pi=({\mathbf{D}^{\top}})^{\dagger}{\mathbf{D}^{\top}}$ and where $\Pi$ is the orthogonal projection onto $\ker({\mathbf{D}^{\top}})$. Furthermore, $\ker({\mathbf{D}^{\top}})$ has dimension 1 so that $\|\Pi\varepsilon\|_{2}^{2}/\sigma^{2}$ is a $\chi^{2}$ random variable with 1 degree of freedom. Thus for any $v\in\mathbb{R}^{n}$ with $\|v\|=1$ we have

[TABLE]

Let us define the function $g(\cdot)$ by

[TABLE]

Then, by the definition of $f,g$ and (2.11), we have almost surely $f(\varepsilon)\leq\|\Pi\varepsilon\|_{2}+g(\varepsilon)$. By a standard bound on $\chi^{2}$ random variable with 1 degree of freedom, we have $\mathbb{P}(\|\Pi\varepsilon\|_{2}\leq\sigma+\sigma\sqrt{2\log(1/\delta)})\geq 1-\delta$. Furthermore, the function $g$ is is 1-Lipschitz and $\varepsilon\sim\mathcal{N}(0,\sigma^{2}\mathrm{Id}_{n})$. By the Gaussian concentration theorem [11, Theorem 10.17], we have

[TABLE]

where $\operatorname{Med}[g(\varepsilon)]$ is the median of the random variable $g(\varepsilon)$. Combining these two probability bounds with the union bound, we obtain $f(\varepsilon)\leq\operatorname{Med}[g(\varepsilon)]+\sigma+2\sigma\sqrt{2\log(1/\delta)}$ with probability at least $1-2\delta$.

To complete the proof, it remains to show that

[TABLE]

By definition of the median, it is enough to show that

[TABLE]

By Lemma A.4 and the fact that $\mathrm{Id}_{n}-\Pi=({\mathbf{D}^{\top}})^{\dagger}{\mathbf{D}^{\top}}$, we obtain that for all $v$,

[TABLE]

where we used the duality between $\|\cdot\|_{[\lambda]}^{*}$ and $\|\cdot\|_{[\lambda]}$ for the second term and the fact that the transpose and the Moore-Penrose pseudo-inverse commute, which implies $(\mathbf{D}^{\dagger})^{\top}={\mathbf{D}^{\top}}^{\dagger}$.

We now bound $g(\varepsilon)$ from above on the event (2.2). On the event (2.2),

[TABLE]

Consider $v\in\mathbb{R}^{n}$ such that $\|v\|=1$ and $3\Lambda(\lambda,s)\|{\mathbf{D}^{\top}}v\|_{2}>\sum_{j=s+1}^{p}\lambda_{j}|{\mathbf{D}^{\top}}v|_{j}^{\downarrow}$. Then, by the definition of $\kappa(s)$ given in defined in (2.4) we have

[TABLE]

Consider $v\in\mathbb{R}^{n}$ such that $\|v\|=1$ and $3\Lambda(\lambda,s)\|{\mathbf{D}^{\top}}v\|_{2}\leq\sum_{j=s+1}^{p}\lambda_{j}|{\mathbf{D}^{\top}}v|_{j}^{\downarrow}$, then we have trivially

[TABLE]

Thus, we have proved that on the event (2.2) that has probability at least $1/2$, (2.15) holds. This implies (2.14) by definition of the median. ∎

The constant $\kappa(s)$ is sometimes referred to as the compatibility factor of ${\mathbf{D}^{\top}}$. Bounds on the compatibility factor are obtained for a large class of random and deterministic graphs [17]. For instance, for graphs with bounded degree, the compatibility factor is bounded from below (see for instance [17, Lemma 3]). In linear regression, constants that measure the correlations of the design matrix have been proposed to study the Lasso and the Dantzig selector: [8] defined the Restricted Eigenvalue constant, [31] defined the Compatibility constant, [35] defined the Cone Invertibility factors and [13] defined the Compatibility factor, to name a few. The Weighted Restricted Eigenvalue constant was also defined in [4] to study the Slope estimator. These constants are the linear regression analogs of $\kappa(s)$ defined in (2.4).

2.1 does not provide an explicit choice for the weights $\lambda_{1}\geq\dots\geq\lambda_{p}$. These weights should be large enough so that the event (2.2) has probability at least $1/2$. These weights should also be as small as possible in order to minimize the right hand side of (2.3). Define $g_{1},\dots,g_{p}$ by

[TABLE]

and let $|g|_{1}^{\downarrow}\geq\dots\geq|g|_{p}^{\downarrow}$ be a nondecreasing rearrangement of $(|g_{1}|,\dots,|g_{p}|)$. Inequality (2.29) below reads

[TABLE]

Thus, if the event

[TABLE]

has probability greater than $1/2$, then the event (2.2) has probability greater than $1/2$ as well, and the conclusion of 2.1 holds. This observation can be used to justify the following heuristics for the choice of the tuning parameters $\lambda_{1}\geq\dots\geq\lambda_{p}$. This heuristics can be implemented provided that the Moore-Penrose pseudo-inverse $\mathbf{D}^{\dagger}$ and the probability distribution of the noise random vector $\varepsilon$ are both known. These heuristics go as follows. Assume that one has generated $N$ independent copies of the random vector $\varepsilon$, and denote by $\mathbb{\tilde{P}}_{N}$ the empirical probability distribution with respect to these independent copies of $\varepsilon$, and $\tilde{F}_{N}^{j}$ the empirical cumulative distribution function (cdf) of $|g|_{j}^{\downarrow}$. Next, define $\lambda_{j}$ as the $(1-1/(3p))$th quantile of $2|g|_{j}^{\downarrow}$, so that

[TABLE]

As $N\to+\infty$, by the Glivenko-Cantelli theorem, $\tilde{F}_{N}^{j}(t)$ converges to the cdf of $|g|_{j}^{\downarrow}$ at $t$ uniformly in $t\in\mathbb{R},j\in[p]$. Hence if $N$ is large enough, then for all $j=1,...,p$ we have $\mathbb{P}(2|g|_{j}^{\downarrow}\leq n\lambda_{j})\geq 1-1/2p\enspace.$ By the union bound over $j=1,\dots,p$,

[TABLE]

thus the event (2.2) has probability greater than $1/2$ with respect to the probability distribution $\mathbb{P}$ of $\varepsilon$. This simple scheme provides a computational heuristic to choose the weights $\lambda_{1},\dots,\lambda_{p}$.

The following corollaries propose a theoretical choice for the weights. To state these corollaries, let us write

[TABLE]

following the notation in [17].

Corollary 1**.**

Assume that the Graph-Slope weights $\lambda_{1}\geq\dots\geq\lambda_{p}\geq 0$ satisfy for any $j\in[p]$

[TABLE]

Then, for any $\delta\in(0,1)$, the oracle inequality (2.3) holds with probability at least $1-2\delta$.

Note that if $\lambda_{1}=\dots=\lambda_{p}=\lambda$, then the event (2.2) reduces to $\|({\mathbf{D}^{\top}})^{\dagger}\varepsilon\|_{\infty}\leq n\lambda/2$. The random variable $\|({\mathbf{D}^{\top}})^{\dagger}\varepsilon\|_{\infty}$ is the maximum of $p$ correlated Gaussian random variables with variance at most $\sigma^{2}\rho(\mathcal{G})^{2}$, so that (2.2) has probability at least $1/2$ provided that $\lambda$ is of order $(\rho(\mathcal{G})\sigma/n)\sqrt{\log p}$.

Proof of 1.

It is enough to show that if the weights $\lambda_{1},\dots,\lambda_{p}$ satisfy (2.26) then the event (2.2) has probability at least $1/2$.

Define $g_{1},\dots,g_{p}$ by (2.23) and let $|g|_{1}^{\downarrow}\geq\dots\geq|g|_{p}^{\downarrow}$ be a nondecreasing rearrangement of $(|g_{1}|,\dots,|g_{p}|)$. For each $j$, the random variable $g_{j}$ is a centered Gaussian random variable with variance at most $\sigma^{2}\|\mathbf{D}^{\dagger}e_{j}\|^{2}\leq\sigma^{2}\rho(\mathcal{G})^{2}$. By definition of the dual norm, and [7, Corrolary II.4.3], we have

[TABLE]

where $V=\sigma^{2}\max_{j=1,\dots,p}\|{\mathbf{D}^{\top}}^{\dagger}e_{j}\|^{2}$. Thus, by Lemma A.5 below, the event (2.2) has probability at least $1/2$. ∎

Under an explicit choice of tuning parameters, 1 yields the following result.

Corollary 2**.**

Under the same hypothesis as 2.1 but with the special choice $n\lambda_{j}=8\sigma\rho(\mathcal{G})\sqrt{\log(2p/j)}$ for any $j\in[p]$, then for any $\delta\in(0,1)$, we have with probability at least $1-2\delta$

[TABLE]

Proof.

We apply Lemma A.2 with the choice $C=8\sigma\rho(\mathcal{G})/n$ ∎

When the true signal satisfies $\|{\mathbf{D}^{\top}}\beta^{\star}\|_{0}=s^{\star}$, the previous bound reduces to

[TABLE]

2 is an improvement w.r.t. the bound provided in [17, Theorem 2] for the TV denoiser (also sometimes referred to as the Generalized Lasso) relying on $\ell_{1}$ regularization defined in Eq. (1.10).

Indeed, the contribution of the second term in 2 is reduced from $\log({ep}/{\delta})$ (in [17, Theorem 2]) to $\log({2ep}/{s})$. Thus the dependence of the right hand side of the oracle inequality in the confidence level $\delta$ is significantly reduced compared to the result of [17, Theorem 2].

A similar bound as in 2 could be obtained for $\ell_{1}$ regularization adapting the proof from [4, Theorem 4.3]. However such a better bound would be obtained for a choice of regularization parameter relying on the ${\mathbf{D}^{\top}}$-sparsity of the signal. The Graph-Slope does not rely on such a quantity, and thus Graph-Slope is adaptive to the unknown ${\mathbf{D}^{\top}}$-sparsity of the signal.

Remark 1**.**

The optimal theoretical choice of parameter requires the knowledge of the noise level $\sigma$ from the practitioner. Whenever the noise level $\sigma$ is not known, the practitioner can use the corresponding Concomitant estimator to alleviate this issue [21, 5, 29], see also [19] for efficient algorithms to compute such scale-free estimators.

3 Numerical experiments

3.1 Algorithm for Graph-Slope

In this section, we propose an algorithm to compute a solution of the highly structured optimization problem (1.7). The data fidelity term $f:\beta\mapsto\|y-\beta\|_{2}^{2}/2$ is a convex smooth function with $1$-Lipschitz gradient, and the map $\beta\mapsto\|{\mathbf{D}^{\top}}\beta\|_{[\lambda]}$ is the pre-composition by a linear operator of the norm $\|\cdot\|_{[\lambda]}$ whose proximal operator can be easily computed [36, 10]. Thus, the use of a dual or primal-dual proximal scheme can be advocated.

Problem (1.7) can be rewritten as

[TABLE]

where $f$ is a smooth, $1$-Lipschitz strictly convex function and $g=\|\cdot\|_{[\lambda]}$ is a convex, proper, lower semicontinuous function (see for instance [2, p. 275]). Its dual problem reads

[TABLE]

where $f^{\star}$ is the convex conjugate of $f$, *i.e., *for any $x\in\mathbb{R}^{n}$

[TABLE]

Classical computations leads to the following dual problem

[TABLE]

The dual formulation (3.4) can be rewritten as an unconstrained problem, using for any set $\mathcal{C}\subset\mathbb{R}^{n}$, and any $\theta\in\mathbb{R}^{n}$, the notation

[TABLE]

The quadratic term in $y$ is constant and can be dropped. Thus the optimization problem (3.4) is equivalent to

[TABLE]

The formulation in  (3.5) is now well suited to apply an accelerated version of the forward-backward algorithm such as FISTA [3]. As a stopping criterion, we use a duality gap criterion: $\Delta(\beta,\theta)\leqslant\varepsilon$, where

[TABLE]

for a feasible pair $(\beta,\theta)$ and by $\Delta(\beta,\theta)=+\infty$ for an unfeasible pair. In practice we set $\varepsilon=10^{-2}$ as a default value. Algorithm 1 summarizes the dual FISTA algorithm applied to the Graph-Slope minimization problem.

We recall that the proximity operator of a convex, proper, lower semicontinuous function $f$ is given as the unique solution of the optimization problem

[TABLE]

To compute the proximity operator of $\iota_{\{\|\cdot\|_{[\lambda]}^{*}\leqslant 1\}}$, we use the Moreau’s decomposition [22, p. 65] which links it to the proximity operator of the dual Slope-norm,

[TABLE]

where $\Pi_{\frac{1}{\tau}B_{*}}$ is the projection onto the unit ball $B_{*}$ associated to the dual norm $\|\cdot\|^{*}$ scaled by a factor $1/\tau$. The proximity operator of $\|\cdot\|_{[\lambda]}$ can be obtained obtained in several ways [36, 10]. In our numerical experiments, we use the connection between this operator and the isotonic regression following [36], which can be computed in linear time. Under the assumption that the quantity $(u_{i}-\lambda_{i})$ is positive, non-increasing (which is obtained by sorting $|u|$ and restoring the signs and ordering afterwards, see details in [36, Eq. (24)]), computing $\operatorname{Prox}_{\|\cdot\|_{[\lambda]}}(u)$ is equivalent to solving the problem

[TABLE]

We have relied on the fast implementation implementation of the PAVA algorithm (for Pool Adjacent Violators Algorithm, see for instance [6]), available in scikit-learn [23] to solve this inner problem.

The source code used for our numerical experiments is freely available at http://github.com/svaiter/gslope_oracle_inequality.

3.2 Synthetic experiments

To illustrate the behavior of Graph-Slope, we first propose two synthetic experiments in moderate dimension. The first one is concerned with the so-called “Caveman” graph and the second one with the 1D path graph.

For these two scenarios, we analyze the performance following the same protocol. For a given noise level $\sigma$, we use the bounds derived in 2.1 (we dropped the constant term 8) and in [17], i.e.,

[TABLE]

For every $n_{0}$ between 0 and $p$, we generate 1000 signals as follows. We draw $J$ uniformly at random among all the subsets of $[p]$ of size $n_{0}$. Then, we let $\Pi_{J}$ be the projection onto $\operatorname{Ker}\mathbf{D}_{J}^{\top}$ and generate a vector $g\sim\mathcal{N}(0,\mathrm{Id}_{n})$. We then construct $\beta^{\star}=c(\mathrm{Id}-\Pi_{J})g$ where $c$ is a given constant (here $c=8$). This constrains the signal $\beta^{\star}$ to be of ${\mathbf{D}^{\top}}$-sparsity at most $p-n_{0}$.

We corrupt the signals by adding a zero mean Gaussian noise with variance $\sigma^{2}$, and run both the Graph-Lasso estimator and the Graph-Slope estimator. We then compute the mean of the mean-squared error (MSE), the false detection rate (FDR) and the true detection rate (TDR). To clarify our vocabulary, given an estimator $\hat{\beta}$ and a ground truth $\beta^{\star}$, the MSE reads $(1/n)\|\beta^{\star}-\hat{\beta}\|^{2}$, while the FDR and TDR read, respectively,

[TABLE]

and

[TABLE]

where for any $z\in\mathbb{R}^{p}$, $\operatorname{supp}(z)=\left\{j\in[p]\;:\;z_{j}\neq 0\right\}$.

Example on Caveman

The caveman model was introduced in [34] to model small-world phenomenon in sociology. Here we consider its relaxed version, which is a graph formed by $l$ cliques of size $k$ (hence $n=lk$), such that with probability $q\in[0,1]$, an edge of a clique is linked to a different clique. In our experiment, we set $l=4$, $k=10$ ($n=40$) and $q=0.1$. We provide a visualisation of such a graph in Figure 1(a). For this realization, we have $p=180$. The rewired edges are indicated in blue in Figure 1(a) whereas the edges similar to the complete graph on 10 nodes are in black. The signals are generated as random vectors of given ${\mathbf{D}^{\top}}$-sparsity with a noise level of $\sigma=0.2$. Figure 1(b) shows the weights decay.

Figures 1(c)–1(e) represent the evolution of the MSE and TDR in function of the level of ${\mathbf{D}^{\top}}$-sparsity. We observe that while the MSE is close between the Graph-Lasso and the Graph-Slope estimator at low level of sparsity, the TDR is vastly improved in the case of Graph-Slope, with a small price concerning the FDR (a bit more for the Monte Carlo choice of the weights). Hence empirically, Graph-Slope will make more discoveries than Graph-Lasso without impacting the overall FDR/MSE, and even improving it.

Example on a path: 1D–Total Variation

The classical 1D–Total Variation corresponds to the Graph-Lasso estimator $\hat{\beta}^{\rm GL}$ when $\mathcal{G}$ is the path graph over $n$ vertices, hence with $p=n-1$ edges. In our experiments, we take $n=100$, $\sigma=0.6$ and a very sparse gradient ($s=4$). According to these values, and taking a random amplitude for each step, we generate a piecewise-constant signal. We display a typical realization of such a signal in Figure 2(a). Figure 2(b) shows the weights decay. Note that in this case, the Monte–Carlo weights shape differs from the one in the previous experiment. Indeed, they are adapated to the underlying graph, contrary to the theoretical weights $\lambda_{\rm GS}$ which depend only on the size of the graph. Figures 2(c)–2(e) represent the evolution of the MSE and TDR in function of the level of ${\mathbf{D}^{\top}}$-sparsity. Here, Graph-Slope does not improve the MSE significantly. However, as for the caveman experiments, Graph-Slope is more likely to make more discoveries than Graph-Lasso for a small price concerning the FDR.

3.3 Example on real data: Paris roads network

To conclude our numerical experiments, we present our result on a real-life graph, the road network of Paris, France. Thanks to the Python module osmnx [9], which downloads and simplifies OpenStreetMap data, we run our experiments on $p=20108$ streets (edges) and $n=10205$ intersections (vertices).

The ground truth signal is constructed as in [16] as follows. Starting from 30 infection sources, each infected intersections has probability $0.75$ to infect each of its neighbors. We let the infection process run for 8 iterations. The resulting graph signal $\beta^{\star}$ is represented in Figure 3(a) with ${\mathbf{D}^{\top}}$-sparsity 1586. We then corrupt this signal by a zero mean Gaussian noise with standard-deviation $\sigma=0.8$, leading to the observations $y$ represented in Figure 3(b).

Instead of using the parameters given in (3.11), we have computed the oracle parameters for the Graph-Lasso and Graph-Slope estimators by evaluating for 100 parameters of the form

[TABLE]

where $\alpha$ lives on a geometric grid inside $[10^{-5},10^{1.5}]$. The best one in terms of MSE (*i.e., *in term of $(1/n)\|\hat{\beta}-\beta^{\star}\|^{2}$) is refered to as the oracle parameter. The results are illustrated in Figure 3(c) for Graph-Lasso and in Figure 3(d) for Graph-Slope. We can see the benefit of Graph-Slope, for instance in the center of Paris where the sources of infections are better identified as shown in the close-up, see Figures 3(e)–3(g).

Appendix A Preliminary lemmas

Lemma A.1**.**

Let $s\in[p]$. For any two $\beta,\hat{\beta}\in\mathbb{R}^{n}$ such that $\|{\mathbf{D}^{\top}}\beta\|_{0}\leqslant s$ we have

[TABLE]

where $u={\mathbf{D}^{\top}}(\hat{\beta}-\beta)$ and \Lambda(\lambda,s)=\Big{(}\sum_{j=1}^{s}\lambda_{j}^{2}\Big{)}^{1/2}.

Proof.

This is a consequence of [4, Lemma A.1]. We provide a proof here for completeness. The second inequality is simple consequence of the Cauchy-Schwarz inequality. Indeed, $(\sum_{j=1}^{s}\lambda_{j}|u|_{j}^{\downarrow})^{2}\leq(\sum_{j=1}^{s}\lambda_{j}^{2})(\sum_{j=1}^{s}(|u|_{j}^{\downarrow})^{2})\leq\|u\|_{2}^{2}(\sum_{j=1}^{s}\lambda_{j}^{2})$.

Let $v={\mathbf{D}^{\top}}\beta$ and $w={\mathbf{D}^{\top}}\hat{\beta}$ and note that $u=w-v$. Let $\varphi$ be a permutation of $[p]$ such that $\|v\|_{[\lambda]}=\sum_{j=1}^{s}\lambda_{j}|v|_{\varphi(j)}$. The first inequality can be rewritten as

[TABLE]

where the supremum is taken over all permutations $\tau$ of $[p]$. We now prove (A.2). Let $\tau$ be a permutation of $[p]$ such that $\varphi(j)=\tau(j)$ for all $j=1,\dots,s$. Then by the triangle inequality, we have

[TABLE]

for each $j=1,\dots,s$ since $u=w-v$. Furthermore, for each $j>s$, it holds that $v_{\tau(j)}=0$ so that $w_{\tau(j)}=u_{\tau(j)}$. Thus,

[TABLE]

It is clear that $\sum_{j=1}^{s}\lambda_{j}|u|_{\tau(j)}\leq\sum_{j=1}^{s}\lambda_{j}|u|_{j}^{\downarrow}$. Finally, notice that it is always possible to find a permutation $\tau$ such that $(|u|_{\tau(j)})_{j>s}$ is non-decreasing. For such choice of $\tau$ we have $-\sum_{j=s+1}^{p}\lambda_{j}|u|_{\tau(j)}\leq-\sum_{j=s+1}^{p}\lambda_{j}|u|_{j}^{\downarrow}$. ∎

Lemma A.2**.**

For the choice of weights: $\forall j\in[p],\lambda_{j}=C\sqrt{\log(2p/j)}$, the following inequalities hold true

[TABLE]

Proof.

Reminding Stirling’s formula $s\log(s/e)\leqslant\log(s!)\leqslant s\log(s)$, one can check that

[TABLE]

∎

Lemma A.3**.**

Let $z,\varepsilon\in\mathbb{R}^{n}$, $y=z+\varepsilon$, and $\hat{\beta}$ a solution of Problem (1.7). Then, for all $\beta\in\mathbb{R}^{n}$,

[TABLE]

Proof.

The objective function of the minimization problem (1.7) is the sum of two convex functions. The first term, *i.e., *the function $\beta\to\frac{1}{2n}\|y-\beta\|^{2}$, is ($1/n$)-strongly convex with respect to the Euclidean norm $\|\cdot\|$. The sum of a 1-strongly convex function and a convex function is 1-strongly convex, and thus by multiplying by $n$ we have

[TABLE]

for all $\beta\in\mathbb{R}^{n}$ and for any $\mathbf{d}$ in the subdifferential of the objective function (1.7) at $\hat{\beta}$. Since $\hat{\beta}$ is a minimizer of $\eqref{eq:gslope}$, we can choose $\mathbf{d}=0$ in the above display. For $\mathbf{d}=0$, the previous display is equivalent to the claim of the Lemma. ∎

Lemma A.4**.**

Let us suppose that the graph $\mathcal{G}$ has $K$ connected components $C_{1},\dots,C_{K}$. Then,

[TABLE]

where for any $k\in[K]$, the vectors $\mathbf{1}_{C_{k}}\in\mathbb{R}^{n}$ are defined by

[TABLE]

Moreover, the orthogonal projection over $\ker({\mathbf{D}^{\top}})$ is denoted by $\Pi$ and is the component-wise averaging given by

[TABLE]

Furthermore, if $\mathcal{G}$ is a connected graph then $\ker(\Pi)=\operatorname{span}(\mathbf{1}_{n})$.

Proof.

The proof can be done for the simple case of a connected graph (i.e., $K=1$), and the result can be generalized by tensorization of graph for $K>1$ components. Hence, we assume that $K=1$. For any $\beta\in\ker({{\mathbf{D}^{\top}}})$, the definition of the incidence matrix yields that for all $(i,j)\in E,\beta_{i}=\beta_{j}$. Since all vertices are connected, by recursion all the $\beta_{j}$’s are identical, and $\beta\in\operatorname{span}(\mathbf{1}_{n})=\operatorname{span}(\mathbf{1}_{C_{1}})$. The converse is proved in the same way. ∎

Lemma A.5** (Proposition E.2 in [4]).**

Let $g_{1},\dots,g_{p}$ be centered Gaussian random variables (not necessarily independent) with variance at most $V>0$. Then,

[TABLE]