On the Global-Local Dichotomy in Sparsity Modeling

TL;DR

This paper addresses the gap between local sparsity models and global signal modeling by proposing a method to construct a global model from local assumptions, ensuring unique and stable recovery.

Contribution

It introduces a framework to derive global signal representations from local sparsity assumptions using constrained linear systems and modern optimization techniques.

Findings

Global models can be constructed from local sparsity assumptions.

Conditions for unique and stable recovery are established.

Numerical experiments support the theoretical results.

Abstract

The traditional sparse modeling approach, when applied to inverse problems with large data such as images, essentially assumes a sparse model for small overlapping data patches. While producing state-of-the-art results, this methodology is suboptimal, as it does not attempt to model the entire global signal in any meaningful way - a nontrivial task by itself. In this paper we propose a way to bridge this theoretical gap by constructing a global model from the bottom up. Given local sparsity assumptions in a dictionary, we show that the global signal representation must satisfy a constrained underdetermined system of linear equations, which can be solved efficiently by modern optimization methods such as Alternating Direction Method of Multipliers (ADMM). We investigate conditions for unique and stable recovery, and provide numerical evidence corroborating the theory.

Equations

$$\sigma \left(D\right):=\min \left\{j:\exists s\subset \left[1,\dots ,m\right],\left|s\right|=j,\mathrm{rank}{D}_s<j\right\}.$$

$$\Vert \alpha {\Vert }_0:=\#\left\{j:{\alpha }_j\ne 0\right\}.$$

$${\mu }_1\left(s\right):=\underset{i\in \left[1,\dots ,m\right]}{\max }\underset{S\subset \left[1,\dots m\right]\setminus \left\{i\right\},|S|=s}{\max }\sum _{j\in S}\left|⟨d_i,d_j⟩\right|.$$

$$\left(1-{\delta }_k\right)\Vert \alpha {\Vert }_2^2⩽\Vert D\alpha {\Vert }_2^2⩽\left(1+{\delta }_k\right)\Vert \alpha {\Vert }_2^2$$

$$R_i:=\left[\begin{array}{ccccccc}0& \dots & 0& I{d}_{n\times n}& 0& \dots & 0\end{array}\right]\in {\mathbb{R}}^{n\times N},$$

$$\mathcal{M}=\mathcal{M}\left(D,s,P,N\right):=\left\{x\in {\mathbb{R}}^N,R_{i}x=D{\alpha }_i,\Vert {\alpha }_i{\Vert }_0⩽s\forall i=1,\dots ,P\right\}.$$

$$R_{i}x=D{\alpha }_i,i=1,\dots ,P,$$

$$x$$

$$D_G:={\left[\left(\frac{1}{n}{R}_i^{T}D\right)\right]}_{i=1,\dots P}\in {\mathbb{R}}^{N\times mP}.$$

$$\Gamma :=\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ {\alpha }_P\end{array}\right].$$

$$x=D_G\Gamma .$$

$$D{\alpha }_i=R_{i}x=R_i{D}_G\Gamma .$$

$$\underset{:=Q_i}{\underbrace{\left[\begin{array}{ccccccc}0& \dots & 0& D& 0& \dots & 0\end{array}\right]}}\Gamma ={\Omega }_i\Gamma ,$$

$$\underset{:=M\in {\mathbb{R}}^{nP\times mP}}{\underbrace{\left[\begin{array}{c}{Q}_1-{\Omega }_1\\ ⋮\\ Q_P-{\Omega }_P\end{array}\right]}}\Gamma =0.$$

$$\dim \ker M=N\left(m-n+1\right).$$

$$\Vert \Gamma {\Vert }_{0,\infty }:=\underset{i=1,\dots ,P}{\max }\Vert {\alpha }_i{\Vert }_0.$$

$$\mathcal{M}=\left\{x\in {\mathbb{R}}^N:x=D_G\Gamma ,M\Gamma =0,\Vert \Gamma {\Vert }_{0,\infty }⩽s\right\}.$$

$$\rho \left(x\right):=\left\{\Gamma \in {\mathbb{R}}^{mP}:\Vert \Gamma {\Vert }_{0,\infty }⩽s,x=D_G\Gamma ,M\Gamma =0,D_{\mathrm{supp}\widetilde{R_i}\Gamma }\text{ is full rank}.\right\}$$

$$A_{\mathcal{S}}:=\left[\begin{array}{c}\left(I_n-P_{s_1}\right)R_1\\ \left(I_n-P_{s_2}\right)R_2\\ ⋮\\ \left(I_n-P_{s_P}\right)R_P\end{array}\right]\in {\mathbb{R}}^{nP\times N}.$$

$${\Sigma }_{\mathcal{M}}:=\left\{\left(s_1,\dots ,s_P\right):\exists x\in \mathcal{M},\Gamma \in \rho \left(x\right)s.t.\forall i=1,\dots ,P:s_i=\mathrm{supp}\widetilde{R_i}\Gamma \right\}.$$

$$\mathcal{M}=\bigcup _{\mathcal{S}\in {\Sigma }_{\mathcal{M}}}\ker A_{\mathcal{S}}.$$

$$S_{T(op)}$$

$$M_{\ast }:=\left[\begin{array}{cccc}{S}_{B}D& -S_{T}D& & \\ & S_{B}D& -S_{T}D& \\ & & \ddots & \ddots \\ & & & \end{array}\right].$$

$$\mathcal{M}=\left\{x\in {\mathbb{R}}^N:x=D_G\Gamma ,M_{\ast }\Gamma =0,\Vert \Gamma {\Vert }_{0,\infty }⩽s\right\}.$$

$$\mathrm{rank}\left[S_B{D}_{s_i}-S_T{D}_{s_{i+1}}\right]<\left|s_i\right|+\left|s_{i+1}\right|⩽2s,$$

$$\left[S_{B}D-S_{T}D\right]\left[\begin{array}{c}{\alpha }_i\\ {\alpha }_{i+1}\end{array}\right]=0.$$

$$\left[S_{B}D-S_{T}D\right]\left[\begin{array}{c}{\alpha }_i\\ {\alpha }_{i+1}\end{array}\right]=\left[S_B{D}_{s_i}-S_T{D}_{s_{i+1}}\right]\left[\begin{array}{c}{\alpha }_i|s_i\\ {\alpha }_{i+1}|s_{i+1}\end{array}\right]=0,$$

$$V=\left\{\left(i_1,\dots ,i_k\right)\subset \left[1,\dots ,m\right]:\mathrm{rank}{D}_{i_1,\dots ,i_k}=k<n\right\},$$

$$E=\biggl{\{}\left(s_{1},s_{2}\right)\in V\times V:\quad\operatorname{rank}\left[S_{B}D_{s_{1}}\quad-S_{T}D_{s_{2}}\right]<\min\left\{n-1,\left|s_{1}\right|+\left|s_{2}\right|\right\}\biggr{\}}.$$

$${\Sigma }_{\mathcal{M}}\subseteq {\mathcal{R}}_{\mathcal{G}}\left(P\right).$$

Full text

1 Introduction

1.1 The need for a new local-global sparsity theory

The sparse representation model elad2010sparseand provides a powerful approach to various inverse problems in image and signal processing such as denoising elad2006image ; mairal2009non , deblurring yu_solving_2012 ; dong2013nonlocally and super-resolution yang2010image ; romano2014single , to name a few mairal2014sparse . This model assumes that a signal can be represented as a sparse linear combination of a few columns (called atoms) taken from a matrix termed dictionary. Consecutively, given a signal, the sparse recovery of its representation over a dictionary is called sparse-coding or pursuit. Due to computational and theoretical aspects, when treating high dimensional data most of the existing sparsity-inspired methods utilize local-patched-based representations rather than the global ones, i.e. they divide a signal into small overlapping blocks (patches), reconstruct these patches using standard sparse recovery techniques, and subsequently average the overlapping regions bruckstein_sparse_2009 ; elad2010sparseand . While this approach leads to highly efficient algorithms producing state-of-the-art results, it is fundamentally limited because the basic sparse model applies to patches only, and does not take into account the dependencies between them.

As an attempt to tackle this flaw, methods based on the notion of structured sparsity eldar_block_2009 ; huang2010benefit ; huang2011learning ; kyrillidis2015structured ; tropp2006algorithms started to appear; for example, in mairal2009non ; dong2013nonlocally ; romano2014single the observation that a patch may have similar neighbors in its surroundings (often termed the self-similarity property) is injected to the pursuit, leading to improved local estimations. Another possibility to consider the dependencies between patches is to exploit the multi-scale nature of the signals mairal2008learning ; sulam2014image ; papyan2016multiscale . A different direction is suggested by the EPLL zoran_learning_2011 ; sulam2015expected ; papyan2016multiscale , which encourages the patches of the final estimate (i.e., after the application of the averaging step) to comply with the local prior. Also, a related work romano2015patch ; romano2015boosting suggests promoting the local estimations to agree on their shared content (the overlap) as a way to achieve a coherent reconstruction of the signal.

Recently, an alternative to the traditional patch-based prior was suggested in the form of the convolutional, or shift-invariant, sparse coding (CSC) model grosse2012shift ; bristow2013fast ; heide2015fast ; gu2015convolutional ; thiagarajan2008shiftinvariant ; rusu2014explicit . Rather than dividing the image into local patches and process each of these independently, this approach imposes a specific structure on the global dictionary – a concatenation of banded circulant matrices – and applies a global pursuit. A thorough theoretical analysis of this model was proposed very recently in papyan2016working_1 ; papyan2016working_2 ; vardan2016convolutional , providing a clear understanding of its success.

The empirical success of the above algorithms indicates the great potential of reducing the inherent gap that exists between the independent local processing of patches and the global nature of the signal at hand. However, a key and highly desirable part is still missing – a theory which would suggest how to modify the basic sparse model to take into account mutual dependence between the patches, what approximation methods to use, and how to efficiently design and learn the corresponding structured dictionary.

1.2 Content and organization of the paper

In this paper we propose a systematic investigation of the signals which are implicitly defined by local sparsity assumptions. A major theme in what follows is that the presence of patch overlaps reduces the number of degrees of freedom, which, in turn, has theoretical and practical implications. In particular, this allows more accurate estimates for uniqueness and stability of local sparse representations, as well as better bounds on performance of existing sparse approximation algorithms. Moreover, the global point of view allows for development of new pursuit algorithms, which consist of local operation on one hand, while also taking into account the patch overlaps on the other hand. Some aspects of the offered theory are still incomplete, and several exciting research direction emerge as well.

The paper is organized as follows. In Section 2 we develop the basic framework for signals which are patch-sparse, building the global model from the “bottom up”, and discuss some theoretical properties of the resulting model. In Section 3 we consider the questions of reconstructing the representation vector, and of denoising a signal in this new framework. We describe “globalized” greedy pursuit algorithms pati1993orthogonal for these tasks, where the patch disagreements play a major role. We show that the frequently used Local Patch Averaging (LPA) approach is in fact suboptimal in this case. In Section 4 we describe several instances/classes of the local-global model in some detail, exemplifying the preceding definitions and results. The examples include piecewise-constant signals, signature-type (periodic) signals, and more general bottom-up models. In Section 5 we present results of extensive numerical experiments, where in particular we show that one of the new globalized pursuits, based on the ADMM algorithm, turns out to have superior performance in all the cases considered. We conclude the paper in Section 6 by discussing possible research directions.

2 Local-global sparsity

We start with the local sparsity assumptions for every patch, and subsequently provide two complimentary characterizations of the resulting global signal space. On one hand, we show that the signals of interest admit a global “sparse-like” representation with a dictionary of convolutional type, and with additional linear constraints on the representation vector. On the other hand, the signal space is in fact a union of linear subspaces, where each subspace is a kernel of a certain linear map. Finally we connect the two points of view by showing that the original local dictionary must carry a combinatorial structure. Concluding this section, we provide some theoretical analysis of the properties of the resulting model, in particular uniqueness and stability of representation. For this task, we define certain measures of the dictionary, similar to the classical spark, coherence function, and the Restricted Isometry Property, which take the additional dictionary structure into account.

2.1 Preliminaries

Definition 1 (Spark of a matrix)

Given a dictionary $D\in\mathbb{R}^{n\times m}$, the spark of $D$ is defined as the minimal number of columns which are linearly dependent:

[TABLE]

Clearly $\sigma\left(D\right)\leqslant n+1$.

Definition 2

Given a vector $\alpha\in\mathbb{R}^{m}$, the $\ell_{0}$ pseudo-norm is the number of nonzero elements in $\alpha$:

[TABLE]

Definition 3

Let $D\in\mathbb{R}^{n\times m}$ be a dictionary with normalized atoms. The $\mu_{1}$ coherence function (Tropp’s Babel’s function) is defined as

[TABLE]

Definition 4

Given a dictionary $D$ as above, the Restricted Isometry constant of order $k$ is the smallest number $\delta_{k}$ such that

[TABLE]

for every $\alpha\in\mathbb{R}^{m}$ with $\|\alpha\|_{0}\leqslant k$.

For any matrix $M$, we denote by $\mathcal{R}\left({M}\right)$ the column space (range) of $M$.

2.2 Globalized local model

In what follows we treat one-dimensional signals $x\in\mathbb{R}^{N}$ of length $N$, divided into $P=N$ overlapping patches of equal size $n$ (so that the original signal is thought to be periodically extended). The other natural choice is $P=N-n+1$, but for simplicity of derivations we consider only the periodic case.

So we define for each $i=1,\dots P$

[TABLE]

the operator extracting $i$-th patch from the signal.

Definition 5

Given local dictionary $D\in\mathbb{R}^{n\times m}$, sparsity level $s<n$, signal length $N$, and the number of overlapping patches $P$, the *globalized local-sparse *model is the set

[TABLE]

This model suggests that each patch, $R_{i}x$ is assumed to have an $s$-sparse representation $\alpha_{i}$, and this way we have charactatized the global $x$ by describing the local nature of its patches.

Next we derive a “global” characterization of ${\cal M}$. Starting with the equations

[TABLE]

and using the equality $Id=\frac{1}{n}\sum_{i=1}^{P}R_{i}^{T}R_{i}$, we have a representation

[TABLE]

Let the global “convolutional” dictionary $D_{G}$ be defined as the horizontal concatenation of the (vertically) shifted versions of $\frac{1}{n}D$, i.e.

[TABLE]

Let $\Gamma\in\mathbb{R}^{mP}$ denote the concatenation of the local sparse codes, i.e.

[TABLE]

Given a vector $\Gamma$ as above, we will denote by $\tilde{R}_{i}$ the operator of extracting its $i$-th portion111Notice that while $R_{i}$ extracts the $i$-th patch from the signal $x$, the operator $\tilde{R_{i}}$ extracts the representation $\alpha_{i}$ of $R_{i}x$ from $\Gamma$. , i.e. $\tilde{R_{i}}\Gamma\equiv\alpha_{i}$.

Summarizing the above developments, we have the global convolutional representation for our signal as follows:

[TABLE]

Next, applying $R_{i}$ to both sides of (5) and using (3), we obtain

[TABLE]

Let $\Omega_{i}:=R_{i}D_{G}$ denote the $i$-th stripe from the global convolutional dictionary $D_{G}$. Thus (6) can be rewritten as

[TABLE]

or $\left(Q_{i}-\Omega_{i}\right)\Gamma=0$. Since this is true for all $i=1,\dots,P$, we have shown that the vector $\Gamma$ satisfies

[TABLE]

Thus, the condition that the patches $R_{i}x$ agree on the overlaps is equivalent to the global representation vector $\Gamma$ residing in the null-space of the matrix $M$.

An easy computation provides the dimension of this null-space (see proof in Appendix A: Proof of Lemma 1), or in other words the overall number of degrees of freedom of admissible $\Gamma$.

Lemma 1

For any frame $D\in\mathbb{R}^{n\times m}$ (i.e. a full rank dictionary), we have

[TABLE]

Note that in particular for $m=n$ we have $\dim\ker M=N$, i.e. every signal admits a unique representation $x=D_{G}\Gamma$ where $\Gamma=\left(D^{-1}R_{1}x,\dots,D^{-1}R_{P}x\right)^{T}$.

Definition 6

Given $\Gamma=\left[\alpha_{1},\dots,\alpha_{P}\right]^{T}\in\mathbb{R}^{mP}$, the $\|\cdot\|_{0,\infty}$ pseudo-norm is defined by

[TABLE]

Thus, every signal complying with the patch-sparse model, with sparsity $s$ for each patch, admits the following representation.

Theorem 2.1

Given $D,s,P$, and $N$, the globalized local-sparse model (3) is equivalent to

[TABLE]

Proof

If $x\in{\cal M}$ (according to (3)), then by the above construction $x$ belongs to the set defined by the RHS of (8) (let’s call it ${\cal M}^{*}$ for the purposes of this proof only). In the other direction, assume that $x\in{\cal M}^{*}$. Now $R_{i}x=R_{i}D_{G}\Gamma=\Omega_{i}\Gamma$, and since $M\Gamma=0$, we have $R_{i}x=Q_{i}\Gamma=D\tilde{R_{i}}\Gamma$. Denote $\alpha_{i}:=\tilde{R_{i}}\Gamma$, and so we have that $R_{i}x=D\alpha_{i}$ with $\|\alpha_{i}\|_{0}\leqslant s$, i.e. $x\in{\cal M}$ by definition. ∎∎

What are the values of $s$ we are interested in? In addition to the natural requirement that $s<n$ (as in Definition 5), we would like to have uniqueness of sparse representations. We say that $\alpha_{i}$ is a minimal representation of $x_{i}$ if $x_{i}=D\alpha_{i}$ such that the matrix $D_{\operatorname{supp}\alpha_{i}}$ has full rank – and therefore the atoms participating in the representation are linearly independent. While we treat uniqueness in more detail in Subsection 2.4, at this point we would like to restrict the discussion to minimal patch representations. Notice that $\alpha_{i}$ might be a minimal representation but not a unique one with minimal sparsity.

Definition 7

Given a signal $x\in{\cal M}$, let us denote by $\rho\left(x\right)$ the set of all locally sparse *and minimal *representations of $x$:

[TABLE]

Let us now go back to the definition (3). Consider a signal $x\in{\cal M}$, and let $\Gamma\in\rho\left(x\right)$. Denote $s_{i}:=\operatorname{supp}\tilde{R_{i}}\Gamma$. Then we have $R_{i}x\in\mathcal{R}\left({D_{s_{i}}}\right)$ and therefore we can write $R_{i}x=P_{s_{i}}R_{i}x$, where $P_{s_{i}}$ is the orthogonal projection operator onto $\mathcal{R}\left({D_{s_{i}}}\right)$. In fact, since $D_{s_{i}}$ is full rank, we have $P_{s_{i}}=D_{s_{i}}D_{s_{i}}^{\dagger}$ where $D_{s_{i}}^{\dagger}=\left(D_{s_{i}}^{T}D_{s_{i}}\right)^{-1}D_{s_{i}}^{T}$ is the Moore-Penrose pseudoinverse of $D_{s_{i}}$.

Definition 8

Given a support sequence ${\cal S}=\left(s_{1},\dots,s_{P}\right)$, define the matrix $A_{{\cal S}}$ as follows:

[TABLE]

The map $A_{{\cal S}}$ measures the local patch discrepancies, i.e. how “far” is each local patch from the range of a particular subset of the columns of $D$.

Definition 9

Given a model ${\cal M},$ denote by $\Sigma_{{\cal M}}$ the set of all valid supports, i.e.

[TABLE]

With this notation in place, it is immediate to see that the global signal model is a union of subspaces.

Theorem 2.2

The global model is equivalent to the union of subspaces

[TABLE]

Remark 1

Contrary to the well-known Union of Subspaces model blumensath_sampling_2009 ; lu_theory_2007 , the subspaces $\left\{\ker A_{{\cal S}}\right\}$ do not have in general a sparse joint basis, and therefore our model is distinctly different from the well-known block-sparsity model eldar_block_2009 ; eldar_robust_2009 .

An important question of interest is to estimate $\dim\ker A_{{\cal S}}$ for a given ${\cal S}\in\Sigma_{{\cal M}}$. One possible solution is to investigate the “global” structure of the corresponding signals (as is done in Subsection 4.1 and Subsection 4.2), while another option is to utilize information about “local connections” (Subsection 2.3 and Subsection 4.4).

2.3 Local support dependencies

In this section we highlight the importance of the local connections (briefly mentioned above) between the neighboring patches of the signal, and therefore between the corresponding subspaces containing those patches. This in turn allows to characterize $\Sigma_{{\cal M}}$ as the set of all “realizable” paths in a certain dependency graph derived from the dictionary $D$. This point of view allows to describe the model ${\cal M}$ using only the intrinsic properties of the dictionary, in contrast to Theorem 2.2.

First we show the equivalence of the condition $M\Gamma=0$ to equality of pairwise overlaps.

Definition 10

Define the “extract from top/bottom” operators $S_{T}\in\mathbb{R}^{\left(n-1\right)\times n}$ and $S_{B}\in\mathbb{R}^{\left(n-1\right)\times n}$:

[TABLE]

The following result is proved in Appendix B: Proof of Lemma 2.

Lemma 2

Let $\Gamma=\left[\alpha_{1},\dots,\alpha_{P}\right]^{T}$. Under the above definitions, the following are equivalent:

$M\Gamma=0;$** 2. 2.

For each $i=1,\dots,P,$ we have $S_{B}D\alpha_{i}=S_{T}D\alpha_{i+1}$.

Definition 11

Let the matrix $M_{*}\in\mathbb{R}^{\left(n-1\right)P\times mP}$ be defined as

[TABLE]

Corollary 1

The global model is equivalent to

[TABLE]

Proposition 1

Let $0\neq x\in{\cal M}$ and $\Gamma\in\rho\left(x\right)$ with $\operatorname{supp}\Gamma=\left(s_{1},\dots,s_{P}\right)$. Then for $i=1,\dots,P$

[TABLE]

where by definition $\operatorname{rank}\emptyset=-\infty$.

Proof

$x\in{\cal M}$ implies by Lemma 2 that for every $i=1,\dots P$

[TABLE]

But

[TABLE]

and therefore the matrix $\left[S_{B}D_{s_{i}}\quad-S_{T}D_{s_{i+1}}\right]$ must be rank-deficient. Note in particular that the conclusion still holds if one (or both) of the $\left\{s_{i},s_{i+1}\right\}$ is empty. ∎∎

The preceding result suggests a way to describe all the supports in $\Sigma_{{\cal M}}$ .

Definition 12

Given a dictionary $D$, we define an abstract directed graph ${\cal G}_{D,s}=\left(V,E\right)$, with the vertex set

[TABLE]

and the edge set

[TABLE]

In particular, $\emptyset\in V$ and $\left(\emptyset,\emptyset\right)\in E$ with $\operatorname{rank}\left[\emptyset\right]:=-\infty$.

Remark 2

It might be impossible to compute ${\cal G}_{D,s}$ in practice. However we set this issue aside for now, and only explore the theoretical ramifications of its properties.

Definition 13

The set of all directed paths of length $P$ in ${\cal G}_{D,s}$, not including the self-loop $\underbrace{\left(\emptyset,\emptyset,\dots\emptyset\right)}_{\times P}$, is denoted by ${\cal C_{G}}\left(P\right)$.

Definition 14

A path ${\cal S}\in{\cal C_{G}}\left(P\right)$ is called realizable if $\dim\ker A_{{\cal S}}>0$. The set of all realizable paths in ${\cal C_{G}}\left(P\right)$ is denoted by ${\cal R_{G}}\left(P\right)$.

Thus we have the following result.

Theorem 2.3

Suppose $0\neq x\in{\cal M}$. Then

Every representation $\Gamma=\left(\alpha_{i}\right)_{i=1}^{P}\in\rho\left(x\right)$ satisfies $\operatorname{supp}\Gamma\in{\cal C_{G}}\left(P\right)$, and therefore

[TABLE] 2. 2.

The model ${\cal M}$ can be characterized “intrinsically” by the dictionary as follows:

[TABLE]

Proof

Let $\operatorname{supp}\Gamma=\left(s_{1},\dots,s_{P}\right)$ with $s_{i}=\operatorname{supp}\alpha_{i}$ if $\alpha_{i}\neq\bm{0}$, and $s_{i}=\emptyset$ if $\alpha_{i}=\bm{0}$. Then by Proposition 1 we must have that

[TABLE]

Furthermore, since $\Gamma\in\rho\left(x\right)$ we must have that $D_{s_{i}}$ is full rank for each $i=1,\dots,P$. Thus $\left(s_{i},s_{i+1}\right)\in{\cal G}_{D,s}$, and so $\operatorname{supp}\Gamma\in{\cal R_{G}}\left(P\right)$. Since by assumption $\operatorname{supp}\Gamma\in\Sigma_{{\cal M}}$, this proves (10).

To show (11), notice that if $\operatorname{supp}\Gamma\in{\cal R_{G}}\left(P\right)$, then for every $x\in\ker A_{\operatorname{supp}\Gamma}$ we have $R_{i}x=P_{s_{i}}R_{i}x$, i.e. $R_{i}x=D\alpha_{i}$ for some $\alpha_{i}$ with $\operatorname{supp}\alpha_{i}\subseteq s_{i}$. Clearly in this case $\left|\operatorname{supp}\alpha_{i}\right|\leqslant s$ and therefore $x\in{\cal M}$. The other direction of (11) follows immediately from the definitions. ∎∎

Definition 15

The dictionary $D$ is called “$\left(s,P\right)$-good” if

[TABLE]

Theorem 2.4

The set of “$\left(s,P\right)$-good” dictionaries has measure zero in the space of all $n\times m$ matrices.

Proof

Every low rank condition defines a finite number of algebraic equations on the entries of $D$ (given by the vanishing of all the $2s\times 2s$ minors of $\begin{bmatrix}S_{B}D_{s_{i}},&S_{T}D_{s_{j}}\end{bmatrix}$ ). Since the number of possible graphs is finite (given fixed $n,m$ and $s$), the resulting solution set is a finite union of semi-algebraic sets of low dimension, and hence has measure zero. ∎∎

The above considerations suggest that the good dictionaries are hard to come by; we provide explicit constructions in Section 4.

Now suppose the graph ${\cal G}$ is known (or can be easily constructed). Then this gives a simple procedure to generate signals from ${\cal M}$, presented in Algorithm 0.1.

An interesting question arises: given ${\cal S}\in{\cal C_{G}}\left(P\right)$, can we say something about $\dim\ker A_{{\cal S}}$? In particular, when is it strictly positive (i.e. when ${\cal S}\in{\cal R_{G}}\left(P\right)$?) While in general the question seems to be difficult, in some special cases this number can be estimated using only the properties of the local connections $\left(s_{i},s_{i+1}\right)$, by essentially counting the additional “degrees of freedom” when moving from patch $i$ to patch $i+1$. We discuss this in more details in Subsection 4.4 (in particular see Proposition 11), while here we show the following easy result.

Proposition 2

For every ${\cal S}\in{\cal R_{G}}\left(P\right)$ we have

[TABLE]

Proof

Notice that

[TABLE]

and therefore $\dim\ker A_{{\cal S}}\leqslant\dim\ker M_{*}^{\left({\cal S}\right)}$. Furthermore, the map D_{G}^{\left({\cal S}\right)}\bigl{|}_{\ker M_{*}^{\left({\cal S}\right)}} is injective, because if $D_{G}^{\left({\cal S}\right)}\Gamma_{{\cal S}}=0$ and $M_{*}^{\left({\cal S}\right)}\Gamma_{{\cal S}}=0$, we must have that $D_{s_{i}}\alpha_{i}|_{s_{i}}=0$ and, since $D_{s_{i}}$ has full rank, also $\alpha_{i}=0$. The conclusion follows. ∎∎

2.4 Uniqueness and stability

Given a signal $x\in{\cal M}$, it has a globalized representation $\Gamma\in\rho\left(x\right)$ according to Theorem 2.1. When is such a representation unique, and under what conditions can it be recovered when the signal is corrupted with noise?

In other words, we study the problem

[TABLE]

and its noisy version

[TABLE]

For this task, we define certain measures of the dictionary, similar to the classical spark, coherence function, and the Restricted Isometry Property, which take the additional dictionary structure into account. In general, the additional structure implies possibly better uniqueness as well as stability to perturbations, however it is an open question to show they are *provably *better in certain cases.

The key observation is that the global model ${\cal M}$ imposes a constraint on the allowed local supports.

Definition 16

Denote the set of allowed local supports by

[TABLE]

Recall the definition of the spark (1). Clearly $\sigma\left(D\right)$ can be equivalently rewritten as

[TABLE]

Definition 17

The *globalized spark *$\sigma^{*}\left(D\right)$ is

[TABLE]

The following proposition is immediate by comparing (12) with (13).

Proposition 3

$\sigma^{*}\left(D\right)\geqslant\sigma\left(D\right).$**

The globalized spark provides a uniqueness result in the spirit of donoho2003optimally .

Theorem 2.5 (Uniqueness)

Let $x\in{\cal M}\left(D,s,N,P\right)$. If $\exists\Gamma\in\rho\left(x\right)$ for which $\|\Gamma\|_{0,\infty}<\frac{1}{2}\sigma^{*}\left(D\right)$ (i.e. it is a sufficiently sparse solution of $P_{0,\infty}$), then it is the unique solution (and so $\rho\left(x\right)=\left\{\Gamma\right\}$).

Proof

Suppose that $\exists\Gamma_{0}\in\rho\left(x\right)$ which is different from $\Gamma$. Put $\Gamma_{1}:=\Gamma-\Gamma_{0}$, then $\|\Gamma_{1}\|_{0,\infty}<\sigma^{*}\left(D\right)$ , while $D_{G}\Gamma_{1}=0$ and $M\Gamma_{1}=0$. Denote $\beta_{j}:=\tilde{R_{j}}\Gamma_{1}$. By assumption, there exists an index $i$ for which $\beta_{i}\neq 0$, but we must have $D\beta_{j}=0$ for every $j$, and therefore $D_{\operatorname{supp}\beta_{i}}$ must be rank deficient – contradicting the fact that $\|\beta_{i}\|<\sigma^{*}\left(D\right)$. ∎∎

In classical sparsity, we have the bound

[TABLE]

where $\mu_{1}$ is given by Definition 3. In a similar fashion, the globalized spark $\sigma^{*}$ can be bounded by an appropriate analog of “coherence” – however, computing this new coherence appears to be in general intractable.

Definition 18

Given the model ${\cal M}$, we define the following globalized coherence function

[TABLE]

where ${\cal T}\cup{\cal T}:=\left\{s_{1}\cup s_{2}:\;s_{1},s_{2}\in{\cal T}\right\}.$

Theorem 2.6

The globalized spark $\sigma^{*}$ can be bounded by the globalized coherence as follows222In general $\min\left\{s:\;\mu_{1}^{*}\left(s-1\right)\geqslant 1\right\}\neq\max\left\{s:\;\mu_{1}^{*}\left(s\right)<1\right\}$ because the function $\mu_{1}^{*}$ need not be monotonic.:

[TABLE]

Proof

Following closely the corresponding proof in donoho2003optimally , assume by contradiction that

[TABLE]

Let $s^{*}\in{\cal T}\cup{\cal T}$ with $\left|s^{*}\right|=\sigma^{*}\left(D\right)$ for which $D_{s^{*}}$ is rank-deficient. Then the restricted Gram matrix $G:=D_{s^{*}}^{T}D_{s^{*}}$ must be singular. On the other hand, $\mu_{1}^{*}\left(\left|s^{*}\right|\right)<1$, and so in particular

[TABLE]

But that means that $G$ is diagonally dominant and therefore $\det G\neq 0$, a contradiction. ∎∎

We see that $\mu_{1}^{*}\left(s+1\right)\leqslant\mu_{1}\left(s\right)$ since the outer maximization is done on a smaller set. Therefore, in general the bound of Theorem 2.6 appears to be sharper than (14).

A notion of globalized RIP can also be defined as follows.

Definition 19

The globalized RIP constant of order $k$ associated to the model ${\cal M}$ is the smallest number $\delta_{k,{\cal M}}$ such that

[TABLE]

for every $\alpha\in\mathbb{R}^{m}$ with $\operatorname{supp}\alpha\in{\cal T}$.

Immediately one can see the following (recall Definition 4).

Proposition 4

The globalized RIP constant is upper bounded by the standard RIP constant:

[TABLE]

Definition 20

The generalized RIP constant of order $k$ associated to signals of length $N$ is the smallest number $\delta_{k}^{\left(N\right)}$ such that

[TABLE]

for every $\Gamma\in\mathbb{R}^{mN}$ satisfying $M\Gamma=0,\;\|\Gamma\|_{0,\infty}\leqslant k$.

Proposition 5

We have

[TABLE]

Proof

Obviously it is enough to show only the leftmost inequality. If $\Gamma=\left(\alpha_{i}\right)_{i=1}^{N}$ and $\|\Gamma\|_{0,\infty}\leqslant k$, this gives $\|\alpha_{i}\|_{0}\leqslant k$ for all $i=1,\dots,P$. Further, setting $x:=D_{G}\Gamma$ we clearly have $\Gamma\in\rho\left(x\right)$ and so $\operatorname{supp}\Gamma\in\Sigma_{{\cal M}}$. Thus $\operatorname{supp}\alpha_{i}\in{\cal T}$, and therefore

[TABLE]

By Corollary 4 we know that for every $\Gamma$ satisfying $M\Gamma=0$, we have

[TABLE]

Now for the lower bound,

[TABLE]

For the upper bound,

[TABLE]

∎∎

Theorem 2.7 (Uniqueness and stability of $P_{0,\infty}$ via RIP)

Suppose that $\delta_{2s}^{\left(N\right)}<1$, and suppose further that $x=D_{G}\Gamma_{0}$ with $\|\Gamma_{0}\|_{0,\infty}=s$ and $\|D_{G}\Gamma_{0}-x\|_{2}\leqslant\varepsilon$. Then every solution $\hat{\Gamma}$ of the noise-constrained $P_{0,\infty}^{\varepsilon}$ problem

[TABLE]

satisfies

[TABLE]

In particular, $\Gamma_{0}$ is the unique solution of the noiseless $P_{0,\infty}$ problem.

Proof

Immediate using the definition of the globalized RIP:

[TABLE]

∎∎

3 Pursuit algorithms

In this section we consider the problem of efficient projection onto the model ${\cal M}$. First we treat the “oracle” setting, i.e. when the supports of the local patches (and therefore of the global vector $\Gamma$) is known. We show that the patch averaging (LPA) method is not a good projector, however repeated application of it does achieve the desired result.

For the non-oracle setting, we consider “local” and “globalized” pursuits. The former type does not use any dependencies between the patches, and tries to reconstruct the supports $\alpha_{i}$ completely locally, using standard methods such as OMP – and as we demonstrate, it can be guaranteed to succeed in more cases than the standard analysis would imply. However a possibly better alternative exists, namely a “globalized” approach with the patch disagreements as a major driving force.

3.1 Global (oracle) projection, local patch averaging (LPA) and the local-global

gap

Here we briefly consider the question of efficient projection onto the subspace $\ker A_{{\cal S}}$, given ${\cal S}$.

As customary in the literature candes2006modernstatistical , the projector onto $\ker A_{{\cal S}}$ can be called *an oracle. *In effect, we would like to compute

[TABLE]

given $y\in\mathbb{R}^{N}$.

To make things concrete, let us assume the standard Gaussian noise model:

[TABLE]

The following is well-known.

Proposition 6

In the Gaussian noise model (16), the performance of the oracle estimator (15) is

[TABLE]

Let us now turn to the LPA method. The (linear part of) LPA is the solution to the minimization problem:

[TABLE]

where $y$ is the noisy signal. This has a closed-form solution

[TABLE]

Again, the following fact is well-established.

Proposition 7

In the Gaussian noise model (16), the performance of the averaging estimator (17) is

[TABLE]

where $\left\{\lambda_{1},\dots,\lambda_{N}\right\}$ are the eigenvalues of $M_{A}M_{A}^{T}$.

Thus, there exists a local-global gap in the oracle setting, illustrated in Figure 1 on page 1. In Subsection 4.1 we estimate this gap for a specific case of piecewise-constant signals.

The following result is proved in Appendix C: Proof of Theorem 3.1.

Theorem 3.1

Iterations of (17) converge to $x_{G}$.

To conclude, the iterated LPA algorithm provides an efficient method for computing the global oracle projection $x_{G}$.

3.2 Local pursuit guarantees

Now we turn to the question of projection onto the model ${\cal M}$ when the support of $\Gamma$ is not known.

Here we show that the local OMP chen1989orthogonal ; pati1993orthogonal in fact succeeds in more cases than can be predicted by the classical unconstrained sparse model for each patch. We use the modified coherence function (which is unfortunately intractable to compute)

[TABLE]

The proof of the following theorem is very similar to proving the guarantee for the standard OMP via the Babel function (Definition 3), see e.g. (foucart2013amathematical, , Theorem 5.14) – and therefore we do not reproduce it here.

Theorem 3.2

If $\eta_{1}^{*}\left(s\right)<1$ then the local OMP will recover the supports of any $x\in{\cal M}$.

Since the modified coherence function takes the allowed local supports into consideration, one can readily conclude that

[TABLE]

and therefore Theorem 3.2 gives in general a better guarantee than the one based on $\mu_{1}$.

3.3 Globalized pursuits

We now turn to consider several pursuit algorithms, aiming at solving the $P_{0,\infty}/P_{0,\infty}^{\varepsilon}$ problems, in the globalized model. The main question is how to project the patch supports onto the nonconvex set $\Sigma_{{\cal M}}$.

The core idea is to relax the constraint $M_{*}\Gamma=0,\;\|\Gamma\|_{0,\infty}\leqslant s$, and allow for some patch disagreements, so that the term $\|M_{*}\Gamma_{k}\|$ is not exactly zero. Intuitive explanation is as follows: the disagreement term “drives” the pursuit, and the probability of success is higher because we only need to “jump-start” it with the first patch, and then by strengthening the weight of the penalty related to this constraint the supports will “align” themselves correctly. Justifying this intuition, at least in some cases, is a future research goal.

Q-OMP

Given $\beta>0$, we define

[TABLE]

The main idea of the Q-OMP algorithm is to substitute the matrix $Q_{\beta}$ as a proxy for the constraint $M_{*}\Gamma=0$, by plugging it as a dictionary to the OMP algorithm. Then, given the obtained support ${\cal S}$, as a way to ensure that this constraint is met, one can construct the matrix $A_{{\cal S}}$ and project the signal onto the subspace $\ker A_{{\cal S}}$ (in Subsection 3.1 we show how such a projection can be done efficiently). The Q-OMP algorithm is detailed in Algorithm 0.2. Let us re-emphasize the point that various values of $\beta$ correspond to different weightings of the model constraint $M_{*}\Gamma=0$, and this might possibly become useful when considering relaxed models (see Section 6).

ADMM-based approach

In what follows we extend the above idea and develop an ADMM-type pursuit boyd2011distributed .

We start with the following global objective:

[TABLE]

Clearly, it is equivalent to $\hat{x}=D_{G}\hat{\Gamma}$, where

[TABLE]

Applying Corollary 4, we have the following result.

Proposition 8

The following problem is equivalent to (18):

[TABLE]

We propose to approximate solution of the nonconvex problem (19) as follows. Define new variables $z_{i}$ (which we would like to be equal to $\alpha_{i}$ eventually), and rewrite the problem in ADMM form (here $Z$ is the concatenation of all the $z_{i}$’s):

[TABLE]

The constraints can be written in concise form

[TABLE]

and so globally we would have the following structure (for $N=3)$

[TABLE]

Our ADMM-based method is defined in Algorithm 0.3.

4 Examples

We now turn to present several classes of signals that belong to the proposed globalized model, where each of these is obtained by imposing a special structure on the local dictionary. Then, we demonstrate how one can sample from ${\cal M}$ and generate such signals.

4.1 Piecewise constant signals

The (unnormalized) Heaviside $n\times n$ dictionary $H_{n}$ is the upper triangular matrix with 1’s in the upper part (see Figure 2 on page 2). Formally, each local atom $d_{i}$ of length $n$ is expressed as a step function, given by $d_{i}^{T}=[\mathbf{1_{i}}\ ,\ \mathbf{0_{n-i}}]^{T}$, $1\leq i\leq n$, where $\mathbf{1_{i}}$ is a vector of ones of length $i$. Similarly, $\mathbb{\mathbf{0}}_{n-i}$ is a zero-vector of length $n-i$. The following property is verified by noticing that $H_{n}^{-1}$ is the discrete difference operator.

Proposition 9

If a patch $x_{i}\in\mathbb{R}^{n}$ has $L-1$ steps, then its (unique) representation in the Heaviside dictionary $H_{n}$ has at most $L$ nonzeros.

Corollary 2

Let $x\in\mathbb{R}^{N}$ be a piecewise-constant signal with at most $L-1$ steps per each segment of length $n$(in the periodic sense). Then

[TABLE]

Remark 3

The model ${\cal M}\left(H_{n},L,N,P=N\right)$ contains also some signals having exactly $L$ steps in a particular patch, but those patches must have their last segment with zero height.

As an example, one might synthesize signals with sparsity $\|\Gamma\|_{0,\infty}\leq 2$ according to the following scheme:

Draw at random the support of $\Gamma$ with the requirement that the distance between the jumps within the signal will be at least the length of a patch (this allows at most 2 non-zeros per patch, one for the step and the second for the bias/DC). 2. 2.

Multiply each step by a random number.

The global subspace $A_{{\cal S}}$ and the corresponding global oracle denoiser $x_{G}$ (15) in the PWC model can be explicitly described.

Proposition 10

Let $x\in\mathbb{R}^{N}$ consist of $s$ constant segments with lengths $\ell_{r}$, $r=1,\dots,s$, and let $\Gamma$ be the (unique) global representation of $x$ in ${\cal M}$ (i.e. $\rho\left(x\right)=\left\{\Gamma\right\}$). Then

[TABLE]

where $B_{r}=\frac{1}{\ell_{r}}\mathbf{1}_{\ell_{r}\times\ell_{r}}$. Therefore, $\dim\ker A_{\operatorname{supp}\Gamma}=s$ and $MSE\left(\hat{x}_{G}\right)=s\sigma^{2}$ under the Gaussian noise model (16).

Proof

Every signal $y\in\ker A_{\operatorname{supp}\Gamma}$ has the same “local jump pattern” as $x$, and therefore it also has the same global jump pattern. That is, every such $y$ consists of $s$ constant segments with lengths $\ell_{r}$. It is an easy observation that such signals satisfy

[TABLE]

This proves (20). It is easy to see that $\dim\ker\left(I_{\ell_{r}}-B_{r}\right)=1$, and therefore

[TABLE]

The proof is finished by invoking Proposition 6. ∎∎

In other words, the global oracle is the averaging operator within the constant segments of the signal, which is quite intuitive.

It turns out that the LPA performance (and the local-global gap) can be accurately described by the following result. We provide an outline of proof in Appendix D: Proof of Theorem 4.1.

Theorem 4.1

Let $x\in\mathbb{R}^{N}$ consist of $s$ constant segments with lengths $\ell_{r}$, $r=1,\dots,s$, and assume the Gaussian noise model (16). Then

There exists a function $R\left(n,\alpha\right):\mathbb{N}\times\mathbb{N}\to\mathbb{R}^{+}$, with $R\left(n,\alpha\right)>1$, such that

[TABLE] 2. 2.

The function $R\left(n,\alpha\right)$ satisfies:

1. (a)

$R\left(n,\alpha\right)=1+\frac{\alpha\left(2\alpha H_{\alpha}^{\left(2\right)}-3\alpha+2\right)-1}{n^{2}}$* if $n\geqslant\alpha$, where $H_{\alpha}^{\left(2\right)}=\sum_{k=1}^{\alpha}\frac{1}{k^{2}}$;* 2. (b)

$R\left(n,\alpha\right)=\frac{11}{18}+\frac{2\alpha}{3n}+\frac{6\alpha-11}{18n^{2}}$* if $n\leqslant\frac{\alpha}{2}$.*

Corollary 3

The function $R\left(n,\alpha\right)$ is monotonically increasing in $\alpha$ (with $n$ fixed) and monotonically decreasing in $n$ (with $\alpha$ fixed). Furthermore,

$\lim_{n\to\infty}R\left(n,n\right)=\frac{\pi^{2}}{3}-2\approx 1.28968$; 2. 2.

$\lim_{n\to\infty}R\left(n,2n\right)=\frac{35}{18}\approx 1.9444.$**

Thus, for reasonable choices of the patch size, the local-global gap is roughly a constant multiple of the number of segments, reflecting the global complexity of the signal.

For numerical examples of reconstructing the PWC signals using our local-global framework, see Subsection 5.2.

4.2 Signature-type dictionaries

Another type of signals that comply with our model are those represented via a signature dictionary, which has been shown to be effective for image restoration aharon2008sparse . This dictionary is constructed from a small signal, $x\in\mathbb{R}^{m}$, such that its every patch (in varying location, extracted in a cyclic fashion), $R_{i}x\in\mathbb{R}^{n}$, is a possible atom in the representation, namely $d_{i}=R_{i}x$. As such, every consecutive pair of atoms $(i,i+1)$ is essentially a pair of overlapping patches that satisfy $S_{B}d_{i}=S_{T}d_{i+1}$ (before normalization). The complete algorithm is presented for convenience in Algorithm 0.4.

Given $D$ as above, one can generate signals $y\in\mathbb{R}^{N}$, where $N$ is an integer multiple of $m$, with $s$ non-zeros per patch, by the easy procedure outlined below.

Init: Construct a base signal $b\in\mathbb{R}^{N}$ by replicating $x\in\mathbb{R}^{m}$ $N/m$ times (note that $b$ is therefore periodic). Set $y=0$. 2. 2.

Repeat for $j=1,\dots,s$:

1. (a)

Shift: Circularly shift the base signal by $t_{j}$ positions, denoted by $\mbox{shift}(b,t_{j})$, for some $t_{j}=0,1,\dots,m-1$ (drawn at random). 2. (b)

Aggregate: y=y+\omega_{j}\mbox{\cdotshift}(b,t_{j}), where $\omega$ is an arbitrary random scalar.

Notice that a signal constructed in this way must be periodic, as it is easily seen that

[TABLE]

while the support sequence ${\cal S}$ is

[TABLE]

Assuming that there are no additional relations between the single atoms of $D$ except those from the above construction, the dependency graph ${\cal G}_{D,1}$ of the resulting dictionary is easily seen to be *cyclic, *and all ${\cal S}\in\Sigma_{{\cal M}}$ are of the above form.

In Figure 3 on page 3 we give an example of a signature-type dictionary $D$ for $\left(n,m\right)=\left(6,10\right)$, its dependency graph ${\cal G}_{D}$, and a signal $x$ with $N=P=30$ together with its corresponding sparse representation $\Gamma$.

Remark 4

It might seem that every $n\times m$ Hankel matrix such as the one shown in Figure 3 on page 3 produces a signature-type dictionary with a nonempty signal space ${\cal M}$. However this is not the case, because such a dictionary will usually fail to generate signals of length larger than $n+m-1$, as its dependency graph will not be cyclic (but rather consist of a single chain of nodes).

Multi-signature dictionaries

One can generalize the construction of Subsection 4.2 and consider $s$-tuples of initial base signals $x_{i},\dots,x_{s}$, instead of a single $x$. The desired dictionary $D$ will consist of corresponding $s$-tuples of atoms, which are constructed from those base signals. In order to avoid ending up with the same structure as the case $s=1$, we also require a “mixing” of the atoms. The complete procedure is outlined in Algorithm 0.5.

In order to generate a signal of length $N$ from ${\cal M}$, one can follow these steps (again we assume that $m$ divides $N$ ):

Create a base signal matrix $X^{G}\in\mathbb{R}^{N\times s}$ by stacking $s\frac{N}{m}$ copies of the original basis matrix $X$. Set $y=0$. 2. 2.

Repeat for $j=1,\dots,k$:

1. (a)

Select a base signal $b_{j}\in\mathcal{R}\left({X^{G}}\right)$ and shift it (in a circular fashion) by some $t_{j}=0,1,\dots,R-1$. 2. (b)

Aggregate: $y=y+\mbox{shift}(b_{j},t_{j})$ (note that here we do not need to multiply by a random scalar).

This procedure will produce a signal $y$ of local sparsity $k\cdot s$. The corresponding support sequence can be written as

[TABLE]

where $s_{i}=s_{1}+i\;\left(\mod m\right)$ and

[TABLE]

Here $\left(t_{j},i\right)$ denotes the atom $d_{t_{j},i}$ in the notation of Algorithm 0.5. The corresponding signal space is

[TABLE]

and it is of dimension $k\cdot s$.

An example of a multi-signature dictionary and corresponding signals may be seen in Figure 4 on page 4.

4.3 Convolutional dictionaries

An important class of signals is the sparse convolution model, where each signal $x\in\mathbb{R}^{N}$ can be written as a linear combination of shifted “waveforms” $\bm{d_{i}}\in\mathbb{R}^{n}$, each $\bm{d_{i}}$ being a column in the local dictionary $D^{\prime}\in\mathbb{R}^{n\times m}$. More conveniently, any such $x$ can be represented as a circular convolution of $\bm{d_{i}}$ with a (sparse) “feature map” $\bm{\psi_{i}}\in\mathbb{R}^{N}$:

[TABLE]

Such signals arise in various applications, such as audio classification blumensath_sparse_2006 ; grosse_shift-invariance_2012 ; smith_efficient_2005 , neural coding ekanadham_unified_2014 ; quiroga_spike_2007 , mid-level image representation and denoising kavukcuoglu_learning_2010 ; zeiler_deconvolutional_2010 ; zeiler_adaptive_2011 .

Formally, the convolutional class can be re-cast into the patch-sparse model of this paper as follows. First, we can rewrite (21) as

[TABLE]

where each $\bm{C}_{i}\in\mathbb{R}^{N\times N}$ is a banded circulant matrix with its first column being equal to $\bm{d_{i}}$, and $\bm{\Psi}\in\mathbb{R}^{Nm}$ is the concatenation of the $\bm{\psi_{i}}$’s. It is easy to see that by permuting the columns of $\bm{E}$ one obtains precisely the global convolutional dictionary $nD_{G}$ based on the local dictionary $D^{\prime}$ (recall (4)). Therefore we obtain

[TABLE]

While it is tempting to conclude from comparing (22) and (5) that the convolutional model is equivalent to the patch-sparse model, an essential ingredient is missing, namely the requirement of equality on overlaps, $M\bm{\Gamma^{\prime}}=0.$ Indeed, nothing in the definition of the convolutional model restricts the representation $\bm{\Psi}$ (and therefore $\bm{\Gamma^{\prime}}),$ therefore in principle the number of degrees of freedom remains $Nm$, as compared to $N\left(m-n+1\right)$ from Proposition 2.

To fix this, following papyan2016working_1 ; papyan2016working_2 we apply $R_{i}$ to (22) and obtain $R_{i}x=R_{i}D_{G}^{\prime}\bm{\Gamma^{\prime}}$. The “stripe” $\Omega_{i}^{\prime}=R_{i}D_{G}^{\prime}$ has only $\left(2n-1\right)m$ nonzero consecutive columns, and in fact the nonzero portion of $\Omega_{i}^{\prime}$ is equal for all $i$. This implies that every $x_{i}$ has a representation $x_{i}=\Theta\bm{\gamma_{i}}$ in the “pseudo-local” dictionary

[TABLE]

where the operators $Z_{B}^{\left(k\right)}$ and $Z_{T}^{\left(k\right)}$ are given by Definition 10 in Appendix B: Proof of Lemma 2. If we now assume that our convolutional signals satisfy

[TABLE]

then we have shown that they belong to ${\cal M}\left(\Theta\left(D^{\prime}\right),s,P,N\right)$, and thus can be formally treated by the framework we have developed.

It turns out that this direct approach is quite naive, as the dictionary $\Theta\left(D^{\prime}\right)$ is extremely bad-equipped for sparse reconstruction (for example it has repeated atoms, and therefore $\mu\left(\Theta\left(D^{\prime}\right)\right)=1$). To tackle this problem, a convolutional sparse coding framework was recently developed in papyan2016working_1 ; papyan2016working_2 , where the explicit dependencies between the sparse representation vectors $\bm{\gamma_{i}}$ (and therefore the special structure of the corresponding constraint $M\left(D^{\prime}\right)\bm{\Gamma^{\prime}}=0$) were exploited quite extensively, resulting in efficient recovery algorithms and nontrivial theoretical guarantees. We refer the reader to papyan2016working_1 ; papyan2016working_2 for further details and examples.

4.4 Arbitrary dependency graphs

The examples considered in the previous sections are somewhat special. For the most general case, one can define an abstract graph ${\cal G}$ with some desirable properties, and subsequently look for a nontrivial realization $D$ of the graph, so that in addition ${\cal R}_{{\cal G}}\neq\emptyset$. Let us therefore discuss each one of those steps, along with a specific example.

Defining an abstract ${\cal G}$ with desirable properties

In this context, we would want ${\cal G}$ to contain sufficiently many different long cycles, which would correspond to long signals and a rich resulting model ${\cal M}$. In contrast with the models from Subsection 4.2, one therefore should allow for some branching mechanism. An example of a possible ${\cal G}$ is given in Figure 5 on page 5. It differs only slightly from the example in Figure 3 on page 3. Notice that due to the structure of ${\cal G}$ there are many possible paths in ${\cal C_{G}}\left(P\right)$. In fact, a direct search algorithm yields $\left|{\cal C_{G}}\left(70\right)\right|=37614$.

Finding a dictionary $D$ which has ${\cal G}$ as its dependency

graph

Every edge in ${\cal G}$ corresponds to a conditions of the form (9) imposed on the entries of $D$. As discussed in Theorem 2.4, this in turn translates to a set of algebraic equations. So the natural idea would be to write out the large system of such equations and look for a solution over the field $\mathbb{R}$ by well-known algorithms in numerical algebraic geometry basu2006algorithms . However, this approach is highly impractical because these algorithms have (single or double) exponential running time. We consequently propose a simplified, more direct approach to the problem.

In detail, we replace the low-rank conditions (9) with more explicit and restrictive ones below.

Assumptions(*)

For each $\left(s_{i},s_{j}\right)\in{\cal G}$ we have $\left|s_{i}\right|=\left|s_{j}\right|=k$. We require that $\operatorname{span}S_{B}D_{s_{i}}=\operatorname{span}S_{T}D_{s_{j}}=\Lambda_{i,j}$ with $\dim\Lambda_{i,j}=k$. Thus there exists a non-singular transfer matrix $C_{i,j}\in\mathbb{R}^{k\times k}$ such that

[TABLE]

In other words, every column in $S_{B}D_{s_{i}}$ must be a specific linear combination of the columns in $S_{T}D_{s_{j}}$. This is much more restrictive than the low-rank condition, but on the other hand, given the matrix $C_{i,j}$, it defines a set of linear constraints on $D$. To summarize, the final algorithm is presented in Algorithm 0.6. In general, nothing guarantees that for a particular choice of ${\cal G}$ and the transfer matrices, there is a nontrivial solution $D$, however in practice we do find such solutions. For example, taking the graph from Figure 5 on page 5 and augmenting it with the matrices $C_{i,j}$ (scalars in this case), we obtain a solution over $\mathbb{R}^{6}$ which is shown in Figure 6 on page 6. Notice that while the resulting dictionary has a Hankel-type structure similar to what we have seen previously, the additional dependencies between the atoms produce a rich signal space structure, as we shall demonstrate in the following section.

The resulting signal space

Given ${\cal G}_{D}$ and the signal length $N=P$, the signals $x$ can be generated according to Algorithm 0.1 0.1. Not all paths in ${\cal C_{G}}$ are realizable, but it turns out that in our example we have $\left|{\cal R_{G}}\left(70\right)\right|=17160$. Three different signals and their supports ${\cal S}$ are shown in Figure 7 on page 7. As can be seen from these examples, the resulting model ${\cal M}$ is indeed much richer than the signature-type construction from Subsection 4.2.

Using the restricted construction of this section, the following estimate can be easily shown.

Proposition 11

Assume that the model satisfies **Assumptions() *above. Then for every ${\cal S}\in{\cal R_{G}}\left(P\right)$

[TABLE]

Proof

The idea is to construct a spanning set for $\ker M_{*}^{\left({\cal S}\right)}$ and invoke Proposition 2. Let us relabel the nodes along ${\cal S}$ to be $1,2,\dots,P$. Starting from an arbitrary $\alpha_{1}$ with support $\left|s_{1}\right|=k$, we use (23) to obtain, for $i=1,2,\dots,P-1$, a formula for the next portion of the global representation vector $\Gamma$

[TABLE]

This gives a set $\Delta$ consisting of overall $k$ linearly independent vectors $\Gamma_{i}$ with $\operatorname{supp}\Gamma_{i}={\cal S}$. It may happen that equation (24) is not satisfied for $i=P$. However, every $\Gamma$ with $\operatorname{supp}\Gamma={\cal S}$ and $M_{*}^{\left({\cal S}\right)}\Gamma_{{\cal S}}=0$ must belong to $\operatorname{span}\Delta$, and therefore

[TABLE]

∎∎

We believe that Proposition 11 can be extended to more general graphs, not necessarily satisfying Assumptions(*).** **In particular, the following estimate appears to hold for a general model ${\cal M}$ and ${\cal S}\in{\cal R_{G}}\left(P\right)$:

[TABLE]

We leave the rigorous proof of this result to a future work.

Further remarks

While the presented model is the hardest to analyze theoretically, even in the restricted case of Assumptions(*) (when does a nontrivial realization of a given ${\cal G}$ exist? how does the answer depend on $n$? When ${\cal R_{G}}\left(P\right)\neq\emptyset$? etc?), we hope that this construction will be most useful in applications such as denoising of natural signals.

5 Numerical experiments

In this section, we test the effectiveness of the globalized model for recovering the signals from Section 4, both in the noiseless and noisy cases. These results are also compared to the classical LPA approach.

5.1 Signature-type signals

In this section we investigate the performance of the pursuit algorithms on signals complying with the signature dictionary model elaborated in Subsection 4.2, constructed from one base signal and allowing for varying values of $s$.

Constructing the dictionary

In the context of the LPA algorithm, the condition for its success in recovering the representation is a function of the mutual coherence of the local dictionary – the smaller this measure the larger the number of non-zeros that are guaranteed to be recovered. Leveraging this, we aim at constructing $D\in\mathbb{R}^{n\times m}$ of a signature type that has a small coherence. This can be cast as an optimization problem

[TABLE]

where $\tilde{D}\left(x\right)$ is computed by Algorithm 0.4 and $\mu$ is the (normalized) coherence function.

In our experiments, we choose $n=15$ and $m=20$, and minimize the above loss function via gradient descent, resulting in $\mu(D(x))=0.26$. We used the TensorFlow open source package tensorflow2015-whitepaper . As a comparison, the coherence of a random signature dictionary is about $0.5.$

Noiseless case

In this setting, we test the ability of the globalized OMP to perfectly recover the sparse representation of clean signature-type signals. Figure 8 compares the proposed algorithm (for different choices of $\beta\in\left\{0.25,0.5,1,2,5\right\}$) with the LPA one by providing their probability of success in recovering the true sparse vectors, averaged over $10^{3}$ randomly generated signals of length $N=100$.

From a theoretical perspective, since $\mu(D)=0.26$, the LPA algorithm is guaranteed to recover the representation when $\|\Gamma\|_{0,\infty}\leq 2$, while as can be seen in practice it successfully recover these for $\|\Gamma\|_{0,\infty}\leq 3$. Comparing the LPA approach to the globalized OMP, one can observe that for $\beta\geq 1$ the latter consistently outperforms the former, having a perfect recovery for $\|\Gamma\|_{0,\infty}\leq 4$. Another interesting insight of this experiment is the effect of $\beta$ on the performance; roughly speaking, a relatively large value of this parameter results in a better success-rate than the very small ones, thereby emphasizing importance of the constraint $M_{*}\Gamma=0$. On the other hand, $\beta$ should not be too large since then the importance of the signal is reduced compared to the constraint, which might lead to deterioration in the success-rate (see the curve that corresponds to$\beta=5$ in Figure 8).

Noisy case

In what follows, the stability of the proposed globalized OMP and the ADMM-pursuit are tested and compared to the traditional LPA algorithm. In addition, we provide the projected versions of these algorithms, which by definition satisfy the constraint of our model. More specifically, given the estimated support $\hat{S}$ of each pursuit, we ensure that the constraint $M_{*}\hat{\Gamma}=0$ is met by constructing the matrix $A_{\hat{S}}$ and then projecting the signal onto the subspace $\ker A_{\hat{{\cal S}}}$. In addition to the above, we provide the restoration performance of the oracle estimator, serving as an indication for the best possible denoising that can be achieved. In this case, the oracle projection matrix $A_{S}$ is constructed according to the ground-truth support $S$.

Per each local cardinality $1\leq\|\Gamma\|_{0,\infty}\leq 4$ we generate 10 random signature-type signals, where each of these is corrupted by white additive Gaussian noise with standard deviation $\sigma$, ranging from 0.05 up-to 0.5. The global number of non-zeros is injected to the globalized OMP, and the information regarding the local sparsity is utilized both by the LPA algorithm and our ADMM-pursuit (which is based on local sparse recovery operations). Following Figure 9, which plots the Mean Squared Error (MSE) of the estimation as a function of the noise level, the ADMM-pursuit achieves the best denoising performance, having similar results to the oracle estimator for all noise-levels and sparsity factors. The source of superiority of the ADMM pursuit might be its inherent ability to obtain an estimation that perfectly fits to the globalized model, i.e., a reconstruction that is identical to its projected version. The second best algorithm being the globalized OMP; for relatively small noise levels, its projected version performs as good as the oracle one, indicating that it successfully recovers the true supports. For large noise levels, however, this algorithm tends to err and results in local estimations that do not “agree” with each other on the overlaps. Yet, the denoising performance of the globalized OMP is better than the one of the LPA algorithm. Notice that the latter performs similarly to the oracle estimator only for very low noise levels and relatively small sparsity factors. This sheds light on the difficulty of finding the true supports, the non trivial solution of this problem, and the great advantage of the proposed globalized model.

Similar conclusion holds for the stable recovery of the sparse representations. Per each pursuit algorithm, Figure 10 illustrates the $\ell_{2}$ distance between the original sparse vector $\Gamma$ and its estimation $\hat{\Gamma}$, averaged over the different noise realizations. As can be seen, the ADMM-pursuit achieves the most stable recovery, the globalized OMP is slightly behind it, and both of them outperform the LPA algorithm especially in the challenging cases of high noise levels and/or large sparsity factors.

5.2 Denoising PWC Signals

In this scenario, we test the ability of the globalized ADMM-pursuit to restore corrupted PWC signals, and compare these to the outcome of the LPA algorithm. Similarly to the previous subsection, the projected versions of the two pursuit algorithms are provided along with the one of the oracle estimator. Following the description in Section 4.1, we generate a signal of length $N=200$, composed of patches of size $n=m=20$ with a local sparsity of at most $2$ non-zeros in the $\ell_{0,\infty}$-sense. These signals are then contaminated by a white additive Gaussian noise with $\sigma$ in the range of $0.1$ to $0.9$.

The restoration performance (in terms of MSE) of the above-mentioned algorithms is illustrated in Figure 11 and the stability of the estimates is demonstrated in Figure 12, where the results are averaged over 10 noise realizations. As can be seen, the globalized approach significantly outperforms the LPA algorithm for all noise levels. Furthermore, when $\sigma\leq 0.5$, the ADMM-pursuit performs similarly to the oracle estimator. One can also notice that, similarly to the previous subsection, the ADMM-pursuit and its projected version resulting in the very same estimation, i.e. this algorithm forces the signal to conform with the patch-sparse model globally. On the other hand, following the visual illustration given in Figure 13, the projected version of the LPA algorithm has zero segments, which are the consequence of a complete disagreement in the support (local inconsistency). This is also reflected in Figure 11, illustrating that even for a very small noise level ($\sigma=0.1$) the projected version of the LPA algorithm has a very large estimation error ($\text{{MSE}}\approx 0.18$) compared to the one of the ADMM-pursuit ($\text{{MSE}}\approx 0.0004$), indicating that the former fails in obtaining a consistent representation of the signal.

6 Discussion

In this work we have presented an extension of the classical theory of sparse representations to signals which are locally sparse, together with novel pursuit algorithms. We envision several promising research directions which might emerge from this work.

6.1 Relation to other models

Viewed globally, the resulting signal model can be considered a sort of “structured sparse” model, however, in contrast to other such constructions (tropp2006algorithms ; huang2010benefit ; huang2011learning ; kyrillidis2015structured and others), our model incorporates both structure in the representation coefficients and a structured dictionary.

The recently developed framework of Convolutional Sparse Coding (CSC) papyan2016working_1 ; papyan2016working_2 ; vardan2016convolutional bears some similarities to our work, in that it, too, has a convolutional representation of the signal via a dictionary identical in structure to $D_{G}$. However, the underlying local sparsity assumptions are drastically different in the two models, resulting in very different guarantees and algorithms. That said, we believe that it would be important to provide precise connections between the results, possibly leading to their deeper understanding. First steps in this direction are outlined in Subsection 4.3.

6.2 Further extensions

The decomposition of the global signal $x\in\mathbb{R}^{N}$ into its patches:

[TABLE]

is a special case of a more general decomposition, namely

[TABLE]

where ${\cal P}_{i}$ is the (orthogonal) projection onto a subspace $W_{i}$ of $\mathbb{R}^{N}$, and $w_{i}$ are some weights. This observation naturally places our theory, at least partially, into the framework of fusion frames, a topic which is generating much interest recently in the applied harmonic analysis community (finiteframes, , Chapter 13). In fusion frame theory, which is motivated by applications such as distributed sensor networks, the starting point is precisely the decomposition (26). Instead of the reconstruction formula $x=\sum_{i}\frac{1}{n}R_{i}^{T}R_{i}x$, in fusion frame theory we have

[TABLE]

where $S_{{\cal W}}$ is the associated fusion frame operator. The natural extension of our work to this setting would seek to enforce some sparsity of the projections ${\cal P}_{i}x$. Perhaps the most immediate variant of (25) in this respect would be to drop the periodicity requirement, resulting in a slightly modified $R_{i}$ operators near the endpoints of the signal. We would like to mention some recent works which investigate different notions of fusion frame sparsity boufounos2011sparserecovery ; aceska2015localsparsity ; ayaz2016uniform .

Another intriguing possible extension of our work is to relax the complete overlap requirement between patches and consider an “approximate patch sparsity” model, where the patch disagreement vector $M\Gamma$ is not zero but “small”. In some sense, one can imagine a full “spectrum” of such models, ranging from a complete agreement (this work) to an arbitrary disagreement (such as in the CSC framework mentioned above).

6.3 Learning models from data

The last point above brings us to the question of how to obtain “good” models, reflecting the structure of the signals at hand (such as speech/images etc.) We hope that one might use the ideas presented here in order to create novel learning algorithms. In this regard, the main difficulty is how to parametrize the space of allowed models in an efficient way. While we presented some initial ideas in Section 4, in the most general case (incorporating the approximate sparsity direction above) the problem remains widely open.

Acknowledgements.

The research leading to these results has received funding from the European Research Council under European Union’s Seventh Framework Programme, ERC Grant agreement no. 320649. The authors would also like to thank Jeremias Sulam, Vardan Papyan, Raja Giryes and Gitta Kutinyok for inspiring discussions.

Appendix A: Proof of Lemma 1

Proof

Denote $Z:=\ker M$ and consider the linear map $A:Z\to\mathbb{R}^{N}$ given by the restriction of the “averaging map” $D_{G}:\mathbb{R}^{mP}\to\mathbb{R}^{N}$ to $Z$.

Let us see first that $im\left(A\right)=\mathbb{R}^{N}$. Indeed, for every $x\in\mathbb{R}^{N}$ consider its patches $x_{i}=R_{i}x$. Since $D$ is full rank, there exist $\left\{\alpha_{i}\right\}$ for which $D\alpha_{i}=x_{i}$. Then setting $\Gamma:=\left(\alpha_{1},\dots,\alpha_{P}\right)$ we have both $D_{G}\Gamma=x$ and $M\Gamma=0$ (by construction, see Section 2), i.e. $\Gamma\in Z$ and the claim follows. 2. 2.

Define

[TABLE]

We claim that $J=\ker A$.

1. (a)

In one direction, let $\Gamma=\left(\alpha_{1},\dots,\alpha_{P}\right)\in\ker A$, i.e. $M\Gamma=0$ and $D_{G}\Gamma=0$. Immediately we see that $\frac{1}{n}D\alpha_{i}=0$ for all $i$, and therefore $\alpha_{i}\in\ker D$ for all $i$, thus $\Gamma\in J$. 2. (b)

In the other direction, let $\Gamma=\left(\alpha_{1},\dots,\alpha_{P}\right)\in J$, i.e. $D\alpha_{i}=0$. Then the local representations agree, i.e. $M\Gamma=0$, thus $\Gamma\in Z$. Furthermore, $D_{G}\Gamma=0$ and therefore $\Gamma\in\ker A$. 3. 3.

By the fundamental theorem of linear algebra we conclude

[TABLE]

∎∎

Appendix B: Proof of Lemma 2

We start with an easy observation.

Proposition 12

For any vector $\rho\in\mathbb{R}^{N}$, we have

[TABLE]

Proof

Since

[TABLE]

we can rearrange the sum and get

[TABLE]

∎∎

Corollary 4

Given $M\Gamma=0$, we have

[TABLE]

Proof

Using Proposition 12, we get

[TABLE]

Now since $M\Gamma=0$, then by definition of $M$ we have $\Omega_{j}\Gamma=D\alpha_{j}$ (see (7)), and this completes the proof. ∎∎

Recall Definition 10 and (2). Multiplying the corresponding matrices gives

Proposition 13

We have the following equality for all $i=1,\dots P$:

[TABLE]

To facilitate the proof, we introduce extension of Definition 10 to multiple shifts as follows.

Definition 21

Let $n$ be fixed. For $k=0,\dots,n-1$ let

$S_{T}^{\left(k\right)}:=\begin{bmatrix}I_{n-k}&\bm{0}\end{bmatrix}$ and $S_{B}^{\left(k\right)}:=\begin{bmatrix}\bm{0}&I_{n-k}\end{bmatrix}$ denote the operators extracting the top (resp. bottom) $n-k$ entries from a vector of length $n$; the matrices have dimension $\left(n-k\right)\times n$. 2. 2.

$Z_{B}^{\left(k\right)}:=\begin{bmatrix}S_{B}^{\left(k\right)}\\ \bm{0}_{k\times n}\end{bmatrix}$ and $Z_{T}^{\left(k\right)}:=\begin{bmatrix}\bm{0}_{k\times n}\\ S_{T}^{\left(k\right)}\end{bmatrix}$ . 3. 3.

$W_{B}^{\left(k\right)}:=\begin{bmatrix}\bm{0}_{k\times n}\\ S_{B}^{\left(k\right)}\end{bmatrix}$ and $W_{T}^{\left(k\right)}:=\begin{bmatrix}S_{T}^{\left(k\right)}\\ \bm{0}_{k\times n}\end{bmatrix}$ .

Note that $S_{B}=S_{B}^{\left(1\right)}$ and $S_{T}=S_{T}^{\left(1\right)}$. We have several useful consequences of the above definitions. The proofs are carried out via elementary matrix identities and are left to the reader.

Proposition 14

For any $n\in\mathbb{N}$ the following hold:

$Z_{T}^{\left(k\right)}=\left(Z_{T}^{\left(1\right)}\right)^{k}$* and $Z_{B}^{\left(k\right)}=\left(Z_{B}^{\left(1\right)}\right)^{k}$ for $k=0,\dots,n-1$;* 2. 2.

$W_{T}^{\left(k\right)}W_{T}^{\left(k\right)}=W_{T}^{\left(k\right)}$* and $W_{B}^{\left(k\right)}W_{B}^{\left(k\right)}=W_{B}^{\left(k\right)}$ for $k=0,\dots,n-1$;* 3. 3.

$W_{T}^{\left(k\right)}W_{B}^{\left(j\right)}=W_{B}^{\left(j\right)}W_{T}^{\left(k\right)}$* for $j,k=0,\dots,n-1$;* 4. 4.

$Z_{B}^{\left(k\right)}=Z_{B}^{\left(k\right)}W_{B}^{\left(k\right)}$* and $Z_{T}^{\left(k\right)}=Z_{T}^{\left(k\right)}W_{T}^{\left(k\right)}$ for $k=0,\dots,n-1$;* 5. 5.

$W_{B}^{\left(k\right)}=Z_{T}^{\left(1\right)}W_{B}^{\left(k-1\right)}Z_{B}^{\left(1\right)}$* and $W_{T}^{\left(k\right)}=Z_{B}^{\left(1\right)}W_{T}^{\left(k-1\right)}Z_{T}$ for $k=1,\dots,n-1$;* 6. 6.

$Z_{B}^{\left(k\right)}Z_{T}^{\left(k\right)}=W_{T}^{\left(k\right)}$* and $Z_{T}^{\left(k\right)}Z_{B}^{\left(k\right)}=W_{B}^{\left(k\right)}$ for $k=0,\dots,n-1$;* 7. 7.

$\left(n-1\right)I_{n\times n}=\sum_{k=1}^{n-1}\left(W_{B}^{\left(k\right)}+W_{T}^{\left(k\right)}\right).$**

Proposition 15

If the vectors $u_{1},\dots,u_{N}\in\mathbb{R}^{n}$ satisfy pairwise

[TABLE]

then they also satisfy for each $k=0,\dots,n-1$ the following:

[TABLE]

Proof

It is easy to see that the condition $S_{B}u_{i}=S_{T}u_{i+1}$ directly implies

[TABLE]

Let us first prove (28) by induction on $k$. The base case $k=1$ is precisely (30). Assuming validity for $k-1$ and $\forall i$, we have

[TABLE]

To prove (29) we proceed as follows:

[TABLE]

This finishes the proof of Proposition 15.

∎∎

Example 1

To help the reader understand the claim of Proposition 15, consider the case $k=2$, and take some three vectors $u_{i},u_{i+1},u_{i+2}.$ We have $S_{B}u_{i}=S_{T}u_{i+1}$ and also $S_{B}u_{i+1}=S_{T}u_{i+2}$. Then clearly $S_{B}^{\left(2\right)}u_{i}=S_{T}^{\left(2\right)}u_{i+2}$ (see Figure 14 on page 14) and therefore $W_{B}^{\left(2\right)}u_{i}=Z_{T}^{\left(2\right)}u_{i+2}$.

Let us now present the proof of Lemma 2.

Proof

We show equivalence in two directions.

• •

$\left(1\right)\Longrightarrow\left(2\right)$: Let $M\Gamma=0$. Define $x:=D_{G}\Gamma$, and then further denote $x_{i}:=R_{i}x$. Then on the one hand:

[TABLE]

On the other hand, because of (27) we have $S_{B}R_{i}x=S_{T}R_{i+1}x$, and by combining the two we conclude that $S_{B}D\alpha_{i}=S_{T}D\alpha_{i+1}$.

• •

$\left(2\right)\Longrightarrow\left(1\right)$: In the other direction, suppose that $S_{B}D\alpha_{i}=S_{T}D\alpha_{i+1}$. Denote $u_{i}:=D\alpha_{i}$. Now consider the product $\Omega_{i}\Gamma$ where $\Omega_{i}=R_{i}D_{G}$. One can easily be convinced that in fact

[TABLE]

Therefore

[TABLE]

Since this holds for all $i$, we have shown that $M\Gamma=0$.

∎∎

Appendix C: Proof of Theorem 3.1

Recall that $M_{A}=\frac{1}{n}\sum_{i}R_{i}^{T}P_{s_{i}}R_{i}$. We first show that $M_{A}$ is a contraction.

Proposition 16

$\left\|M_{A}\right\|_{2}\leqslant 1$.

Proof

Closely following a similar proof in romano2015boosting , divide the index set $\left\{1,\dots,N\right\}$ into $n$ groups representing *non-overlapping *patches: for $i=1,\dots,n$ let

[TABLE]

Now

[TABLE]

By construction, $R_{j}R_{k}^{T}=\bm{0}_{n\times n}$ for $j,k\in K\left(i\right)$ and $j\neq k$. Therefore for all $i=1,\dots,n$ we have

[TABLE]

Substituting in back into the preceding inequality finally gives

[TABLE]

∎∎

Now let us move on to prove Theorem 3.1.

Proof

Define

[TABLE]

It is easy to see that

[TABLE]

Let the SVD of $A_{{\cal S}}$ be

[TABLE]

Now

[TABLE]

Therefore $T=nI-V\Sigma V^{T}$, and

[TABLE]

This shows that the eigenvalues of $M_{A}$ are $\tau_{i}=1-\frac{\sigma_{i}^{2}}{n}$ where $\left\{\sigma_{i}\right\}$ are the singular values of $A_{{\cal S}}$. Thus we obtain

[TABLE]

If $\sigma_{i}=0$ then $\tau_{i}=1$, and in any case, by Proposition 16 we have $\left|\tau_{i}\right|\leqslant 1$. Let the columns of the matrix $W$ consist of the singular vectors of $A_{{\cal S}}$ corresponding to $\sigma_{i}=0$ (and so $\operatorname{span}W=\mathcal{N}\left({A_{{\cal S}}}\right)$), then

[TABLE]

Thus, as $k\to\infty$, $M_{A}^{k}$ tends to the orthogonal projector onto $\mathcal{N}\left({A_{{\cal S}}}\right)$. ∎∎

Appendix D: Proof of Theorem 4.1

Recall that the signal consists of $s$ constant segments of corresponding lengths $\ell_{1},\dots,\ell_{s}$. We would like to compute the MSE for every pixel within every such segment of length $\alpha:=\ell_{r}$. For each patch, the oracle provides the locations of the jump points within the patch.

Let us calculate the MSE for pixel with index [math] inside a constant (nonzero) segment $\left[-k,\alpha-k-1\right]$ with value $v$ (Figure 15 on page 15 might be useful). The oracle estimator has the explicit formula

[TABLE]

where $j=1,\dots,n$ corresponds to the index of the overlapping patch containing the pixel, intersecting the constant segment on $\left[a_{j},b_{j}\right]$, so that

[TABLE]

Now, the oracle error for the pixel is

[TABLE]

where the coefficients $c_{i,\alpha,n,k}$ are some positive rational numbers depending only on $i,\alpha,n$ and $k$. It is easy to check by rearranging the above expression that

[TABLE]

and furthermore, denoting $d_{i}:=c_{i,\alpha,n,k}$ for fixed $\alpha,n,k$, we also have that

[TABLE]

Example 2

$n=4,\;\alpha=3$

• •

For $k=1$:

[TABLE]

• •

For $k=2$:

[TABLE]

Now consider the optimization problem

[TABLE]

It can be easily verified that it has the optimal value $\frac{1}{\alpha}$, attained at $c^{*}=\alpha\mathbf{1}$. From this, (32) and (33) it follows that

[TABLE]

Since the $z_{i}$ are i.i.d., we have

[TABLE]

while for the entire nonzero segment of length $\alpha=\ell_{r}$

[TABLE]

Defining

[TABLE]

we obtain that $R\left(n,\alpha\right)>1$ and furthermore

[TABLE]

This proves item $\left(1\right)$ of Theorem 4.1. For showing the explicit formulas for $R\left(n,\alpha\right)$ in item $\left(2\right)$, we have used automatic symbolic simplification software MAPLE maplesoft .

By construction (31), it is not difficult to see that if $n\geqslant\alpha$ then

[TABLE]

where $H_{k}:=\sum_{i=1}^{k}\frac{1}{i}$ is the $k$-th harmonic number. This simplifies to

[TABLE]

where $H_{k}^{\left(2\right)}=\sum_{i=1}^{k}\frac{1}{i^{2}}$ is the $k$-th harmonic number of the second kind.

On the other hand, for $n\leqslant\frac{\alpha}{2}$ we have

[TABLE]

where

[TABLE]

and

[TABLE]

Automatic symbolic simplification of the above gives

[TABLE]