Error analysis for denoising smooth modulo signals on a graph

Abstract

In many applications, we are given access to noisy modulo samples of a smooth function with the goal being to robustly unwrap the samples, i.e., to estimate the original samples of the function. In a recent work, Cucuringu and Tyagi [11] proposed denoising the modulo samples by first representing them on the unit complex circle and then solving a smoothness regularized least squares problem – the smoothness measured w.r.t. the Laplacian of a suitable proximity graph G – on the product manifold of unit circles. This problem is a quadratically constrained quadratic program (QCQP) which is nonconvex, hence they proposed solving its sphere-relaxation leading to a trust region subproblem (TRS). In terms of theoretical guarantees, ${\ell }_2$ error bounds were derived for (TRS). These bounds are however weak in general and do not really demonstrate the denoising performed by (TRS).

In this work, we analyse the (TRS) as well as an unconstrained relaxation of (QCQP). For both these estimators we provide a refined analysis in the setting of Gaussian noise and derive noise regimes where they provably denoise the modulo observations w.r.t. the ${\ell }_2$ norm. The analysis is performed in a general setting where G is any connected graph.

Introduction

Many modern applications involve the acquisition of noisy modulo samples of a function $f:{\mathbb{R}}^d\to \mathbb{R}$, i.e., we obtain

$$y_i=(f(x_i)+{\eta }_i)\mathrm{mod}\zeta ;i=1,\dots ,n$$

for some $\zeta \in {\mathbb{R}}^{+}$, where ${\eta }_i$ denotes noise. Here, $a\mathrm{mod}\zeta$ lies in the interval $[0,\zeta )$ and is such that $a=q\zeta +(a\mathrm{mod}\zeta )$ for an integer q. This situation arises, for instance, in self-reset analog to digital converters (ADCs) which handle voltage surges by simply storing the modulo value of the voltage signal [17], [26], [35]. In other words, if the voltage signal exceeds the range $[0,\zeta ]$, then its value is simply reset via the modulo operation. Another important application is phase unwrapping where one typically obtains noisy modulo 2π samples. There, one usually seeks to infer the structure of an object by transmitting waveforms, and measuring the difference in phase (in radians) between the transmitted and scattered waveforms. This is common in synthetic radar aperture interferometry (InSAR) where one aims to learn the elevation map of a terrain (see for e.g. [13], [36]), and also arises in many other domains such as MRI [15], [20] and diffraction tomography [25], to name a few.

Let us assume $\zeta =1$ from now. Given (1.1), one can ask whether we can recover the original samples of f, i.e., $f(x_i)$? Clearly, this is only possible up to a global integer shift. Furthermore, answering this question requires making additional assumptions about f, for instance, we may assume f to be smooth (for e.g., Lipschitz, continuously differentiable etc.). In this setting, when $d=1$, it is not difficult to see that we can exactly recover the samples of f when there is no noise (i.e., ${\eta }_i=0$ in (1.1)) provided the sampling is fine enough. This is achieved by sequentially traversing the samples and reconstructing the quotient $q_i\in \mathbb{Z}$ corresponding to $f(x_i)$ by setting it to $q_{i-1}+a_{i-1,i}$ where $a_{i-1,i}$ equals either (a) 0 (if $y_i$ is “close enough” to $y_{i-1}$), or (b) ±1 (if $y_i$ is sufficiently smaller/larger than $y_{i-1}$). If f is smooth, then a fine enough sampling density will ensure that nearby samples of f have quotients which differ by either $1,-1$ or 0. This argument can be extended to the general multivariate setting as well (see [12]).

While exact recovery of $f(x_i)$ is of course no longer possible in the presence of noise, one can show [12] that if ${\eta }_i\mathrm{mod}1$ is not too large (uniformly for all i), and if n is sufficiently large, then the estimates $\tilde{f}(x_i)$ returned by the above procedure are such that (up to a global integer shift) for each i, $|\tilde{f}(x_i)-f(x_i)|$ is bounded by a term proportional to the noise level. This suggests a natural two stage procedure for estimating $f(x_i)$ – first obtain denoised estimates of $f(x_i)\mathrm{mod}1$, and in the second stage, “unwrap” these denoised modulo samples to obtain the estimates $\tilde{f}(x_i)$.

This was the motivation behind the recent work of Cucuringu and Tyagi [10], [11] that focused primarily on the first (modulo denoising) stage, which is an interesting question by itself. Before discussing their result, it will be convenient to fix the notation used throughout the paper.

Notation. Vectors and matrices are denoted by lower and upper case symbols respectively, while sets are denoted by calligraphic symbols (e.g., $\mathcal{N}$), with the exception of $[n]=\{1,\dots ,n\}$ for $n\in \mathbb{N}$. The imaginary unit is denoted by $\iota =\sqrt{-1}$. For $a,b\ge 0$, we write $a\lesssim b$ when there exists an absolute constant $C>0$ such that $a\le Cb$. If $a\lesssim b$ and $a\gtrsim b$, then we write $a\asymp b$. For $u\in {\mathbb{C}}^n$, we write $u=\text{Re}(u)+\iota \text{Im}(u)$ where $\text{Re}(u),\text{Im}(u)\in {\mathbb{R}}^n$ denote its real and imaginary parts respectively. The Hermitian conjugate of u is denoted by $u^{⁎}$, and ${\Vert u\Vert }_p$ denotes the usual ${\ell }_p$ norm in ${\mathbb{C}}^n$ for $1\le p\le \infty$.

Cucuringu and Tyagi [11], [10] considered the problem of denoising modulo samples by formulating it as an optimization problem on a manifold. Assume $f:{[0,1]}^d\to \mathbb{R}$, and suppose that $x_i$'s form a uniform grid in ${[0,1]}^d$. The main idea is to represent each sample $y_i$ via $z_i=\exp (\iota 2\pi y_i)$ and observing that if $x_i$ is close to $x_j$, then $\exp (\iota 2\pi f(x_i))\approx \exp (\iota 2\pi f(x_j))$ holds provided f is smooth. In other words, the samples $z=(z_1,\dots ,z_n)$ lie on the manifold

$${\mathcal{C}}_n:=\{u\in {\mathbb{C}}^n:\left|u_i\right|=1;i=1,\dots ,n\}$$

which is the product manifold of unit radius circles, i.e., ${\mathcal{C}}_n={\mathcal{C}}_1\times \cdots \times {\mathcal{C}}_1$. Hence, by constructing a proximity graph $G=([n],E)$ on the sampling points ${(x_i)}_{i=1}^n$ – where there is an edge between i and j if $x_i$ is within a specified distance of $x_j$ – they proposed solving a quadratically constrained quadratic program (QCQP)

$$\underset{g\in {\mathcal{C}}_n}{\min }{\Vert g-z\Vert }_2^2+\gamma g^{⁎}Lg$$

where L is the Laplacian of G, and $\gamma >0$ is a regularization parameter. While the objective function is convex, this is a non-convex problem due to the constraint set ${\mathcal{C}}_n$. It is not clear whether the global solution of (QCQP) can be obtained in polynomial time. Hence they proposed relaxing ${\mathcal{C}}_n$ to a sphere leading to the trust region subproblem (TRS)

$$\underset{g\in {\mathbb{C}}^n:{\Vert g\Vert }_2^2=n}{\min }{\Vert g-z\Vert }_2^2+\gamma g^{⁎}Lg\equiv \underset{g\in {\mathbb{C}}^n:{\Vert g\Vert }_2^2=n}{\min }-2\text{Re}(g^{⁎}z)+\gamma g^{⁎}Lg.$$

Trust region subproblems are well known to be solvable in polynomial time, with many efficient solvers (e.g., [14], [1]). In [11], [10], (TRS) was shown to be quite robust to noise and scalable to large problem sizes through extensive experiments. On the theoretical front, a ${\ell }_2$ stability result was shown for the solution $\hat{g}$ of (TRS). To describe this result, denote $h\in {\mathcal{C}}_n$ to be the ground truth samples where $h_i=\exp (\iota 2\pi f(x_i))$, and define the entry-wise projection operator ${\mathbb{C}}^n\mapsto {\mathcal{C}}_n$ [21]

$${\left(\frac{g}{\left|g\right|}\right)}_i:=\{\begin{array}{cc}\hfill \frac{g_i}{\left|g_i\right|};& g_i\ne 0,\hfill \\ \hfill 1;& \text{ otherwise, }\hfill \end{array}i=1,\dots ,n.$$

Consider the univariate case with $f:[0,1]\to \mathbb{R}$, and assume that the (Gaussian) noise ${\eta }_i\sim \mathcal{N}(0,{\sigma }^2)$ i.i.d. for $i=1,\dots ,n$. Assuming f is M-Lipschitz, and $\sigma \lesssim 1$, denote Δ to be the maximum degree of G with $x_i$'s forming a uniform grid on $[0,1]$. Then provided $\gamma \mathrm{\Delta }\lesssim 1$, the bounds in [11, Theorem 14] imply¹ that with high probability, the solution $\hat{g}$ of (TRS) satisfies

$${\Vert \frac{\hat{g}}{\left|\hat{g}\right|}-h\Vert }_2^2\lesssim \sigma n+\frac{\gamma M^2{\mathrm{\Delta }}^3}{n}.$$

On the other hand, one can show that ${\Vert z-h\Vert }_2^2\asymp {\sigma }^{2}n$ with high probability if $\sigma \lesssim 1$, hence we cannot conclude

$${\Vert \frac{\hat{g}}{\left|\hat{g}\right|}-h\Vert }_2\ll {\Vert z-h\Vert }_2.$$

This motivates the present paper which seeks to identify conditions under which (1.4) holds. In fact, we will consider a more abstract problem setting where G is any connected graph, and $h\in {\mathcal{C}}_n$ is smooth w.r.t. G. This is formally described below, followed by a summary of our main results.

Consider an undirected, connected graph $G=([n],E)$ where $E\subseteq \{\{i,j\}:i\ne j\in [n]\}$ denotes its set of edges. Denote the maximum degree of G by △, and the (combinatorial) Laplacian matrix associated with G by $L\in {\mathbb{R}}^{n\times n}$. We will denote the eigenvalues of L by

$${\lambda }_n=0<{\lambda }_{n-1}={\lambda }_{\min }\le {\lambda }_{n-2}\le \cdots \le {\lambda }_1$$

where ${\lambda }_{\min }$ denotes the Fiedler value of G, and is well known to be a measure of connectivity of G. The corresponding (unit ${\ell }_2$ norm) eigenvectors of L will be denoted by $q_j\in {\mathbb{R}}^n$; $j=1,\dots ,n$ where we have that $q_n=\frac{1}{\sqrt{n}}{(1,\dots ,1)}^T$.

Let $h={[h_1,\dots ,h_n]}^T\in {\mathcal{C}}_n$ be an unknown ground truth signal which is assumed to be smooth w.r.t. G in the sense that

$$h^{⁎}Lh=\sum _{\{i,j\}\in E}{|h_i-h_j|}^2\le B_n.$$

Here, $B_n$ will typically be “small” with the subscript n depicting possible dependence on n. We will assume that information about h is available in the form of noisy $z\in {\mathcal{C}}_n$. Specifically, denoting ${\eta }_i\sim \mathcal{N}(0,{\sigma }^2)$; $i=1,\dots ,n$ to be i.i.d. Gaussian random variables, we are given

$$z_i=h_i\exp (\iota 2\pi {\eta }_i);i=1,\dots ,n.$$

Given access to the noisy samples $z\in {\mathcal{C}}_n$, and the graph G, we aim to answer the following two questions.

• Consider the unconstrained quadratic program (UCQP) obtained by relaxing the constraints in (QCQP) to ${\mathbb{C}}^n$

  $$\underset{g\in {\mathbb{C}}^n}{\min }{\Vert g-z\Vert }_2^2+\gamma g^{⁎}Lg.$$

  This is a convex program, also referred to as Tikhonov regularization in the literature (e.g., [31]), and its solution is given in closed form by $\hat{g}={(I+\gamma L)}^{-1}z$. Under what conditions can we ensure that $\hat{g}$ satisfies (1.4)?

• Under what conditions can we ensure that the solution $\hat{g}$ of (TRS) satisfies (1.4)?

As we will see in Section 4, the solution of (TRS) is closely related to that of (UCQP) but requires a more careful analysis.

Remark 1

It is worth mentioning that the quadratic penalty term $\gamma g^{⁎}Lg$ in (UCQP) and (TRS) aims to promote solutions $\hat{g}$ which are smooth w.r.t G. This is natural given that the ground truth signal $h\in {\mathcal{C}}_n$ is assumed to be smooth w.r.t. G, as stated in (1.5). Moreover, the choice of the regularization parameter γ is important since larger values of γ increase the smoothness of the estimate (larger bias), while smaller values increase the variance of the estimate. For e.g., if we set $\gamma =0$, which is equivalent to removing the regularization term, then we obtain $\hat{g}=z$. The main point of the ensuing analysis is to find a suitable “intermediate” choice of γ which ensures that $\hat{g}$ satisfies (1.4).

Example: denoising modulo samples of a function. We will also seek to instantiate the above results for the setup in [11] described earlier. More precisely, denoting $f:[0,1]\to \mathbb{R}$ to be a M-Lipschitz function where $|f(x)-f(y)|\le M|x-y|$ for each $x,y\in [0,1]$, let $x_i=\frac{i-1}{n-1}$ for $i=1,\dots ,n$ be a uniform grid. For each i, we are given

$$y_i=(f(x_i)+{\eta }_i)\mathrm{mod}1;{\eta }_i\sim \mathcal{N}(0,{\sigma }^2),$$

where ${\eta }_i$ are i.i.d. Then denoting $z_i=\exp (\iota 2\pi y_i)$ and $h_i=\exp (\iota 2\pi f(x_i))$, we will consider $G=([n],E)$ to be the path graph ($P_n$), where

$$E=\{\{i,i+1\}:i=1,\dots ,n\}.$$

Before stating our results, let us define for any $\lambda \in [{\lambda }_{\min },{\lambda }_1]$, the set

$${\mathcal{L}}_{\lambda }:=\{j\in [n-1]:{\lambda }_j<\lambda \}$$

consisting of indices corresponding to the “low frequency” part of the spectrum of L. Moreover, while all our results are non-asymptotic, we suppress constants in this section for clarity. The reader is referred to Appendix A for a tabular summary of the main notation used in the paper.

Results for (UCQP). We begin by outlining our main results for the estimator (UCQP). Theorem 1 below identifies conditions on σ under which (UCQP) provably denoises the samples in expectation. The statement is a simplified version of Theorem 4.

Theorem 1

For any given $\epsilon \in (0,1)$ and $\bar{\lambda }\in [{\lambda }_{\min },{\lambda }_1]$ satisfying $1+\left|{\mathcal{L}}_{\bar{\lambda }}\right|\lesssim \epsilon n$, suppose that $\frac{1}{\epsilon \bar{\lambda }}\sqrt{\frac{△B_n}{n}}\lesssim \sigma \lesssim 1$. Then for the choice $\gamma \asymp {(\frac{{\sigma }^{2}n}{△B_n{\bar{\lambda }}^2})}^{1/4}$, the solution $\hat{g}$ of (UCQP) satisfies $\mathbb{E}{\Vert \frac{\hat{g}}{\left|\hat{g}\right|}-h\Vert }_2^2\le \epsilon \mathbb{E}{\Vert z-h\Vert }_2^2$.

The parameter $\bar{\lambda }$ depends on the spectrum of the Laplacian and can be thought of as the “cut-off frequency”. As can be seen from Theorem 1, we would ideally like $\bar{\lambda }$ to be “large” and $\left|{\mathcal{L}}_{\bar{\lambda }}\right|$ to be “small”; in particular, $\left|{\mathcal{L}}_{\bar{\lambda }}\right|=o(n)$. For e.g., when G is the complete graph² ($K_n$), then ${\lambda }_n=0$ and ${\lambda }_i=n$ for $i=1,\dots ,n-1$. Hence, the choice $\bar{\lambda }=n$ is ideal as it leads to $\left|{\mathcal{L}}_{\bar{\lambda }}\right|=0$. Furthermore, the noise regime in the Theorem involves a lower bound on σ which might seem unnatural. We believe this is likely an artefact of the analysis involving the estimation error. Specifically, the error bound we obtain is of the form (see Lemma 3.1)

$$\mathbb{E}{\Vert \frac{\hat{g}}{\left|\hat{g}\right|}-h\Vert }_2^2\lesssim \frac{\sigma }{\bar{\lambda }}\sqrt{△B_{n}n}+{\sigma }^2(1+\left|{\mathcal{L}}_{\bar{\lambda }}\right|).$$

The problematic term is the first term which depends linearly on σ. Indeed, since ${\Vert z-h\Vert }_2^2\asymp {\sigma }^{2}n$ w.h.p., we end up with the lower bound requirement when $\sigma \ll 1$ for ensuring the denoising guarantee. Nevertheless, if $\frac{1}{\epsilon \bar{\lambda }}\sqrt{\frac{△B_n}{n}}=o(1)$ as n increases, the requirement on σ is of the form $o(1)\le \sigma \lesssim 1$ which becomes progressively mild as n increases.

We also derive conditions under which denoising occurs with high probability³ (w.h.p.). Theorem 2 below is a simplified version of Theorem 6.

Theorem 2

For any given $\epsilon \in (0,1)$ and $\bar{\lambda }\in [{\lambda }_{\min },{\lambda }_1]$, suppose that

$$\max \{\frac{1}{\epsilon \bar{\lambda }}\sqrt{\frac{△B_n}{n}},\frac{\log n}{\sqrt{\epsilon n}}\}\lesssim \sigma \lesssim 1\mathit{\text{and}}1+\left|{\mathcal{L}}_{\bar{\lambda }}\right|+\sqrt{(1+\left|{\mathcal{L}}_{\bar{\lambda }}\right|)\log n}\lesssim \epsilon n.$$

If $\gamma \asymp {(\frac{{\sigma }^{2}n}{△B_n{\bar{\lambda }}^2})}^{1/4}$, then w.h.p., the solution $\hat{g}$ of (UCQP) satisfies ${\Vert \frac{\hat{g}}{\left|\hat{g}\right|}-h\Vert }_2^2\le \epsilon {\Vert z-h\Vert }_2^2$.

The conditions in Theorem 2 are slightly more stringent compared to those in Theorem 1 due to the appearance of extra log factors. These arise due to the concentration inequalities used in the analysis for bounding the estimation error, see Theorem 5.

In Section 3, we interpret our results for special graphs⁴ such as $K_n,P_n$ and the star graph⁵ ($S_n$); see Corollary 2, Corollary 3. In particular, the statement for $P_n$ therein can be applied to the example from Section 1.2, readily yielding conditions that ensure denoising; see Corollary 4. It states that if $\frac{\log n}{\sqrt{n}}\lesssim \sigma \lesssim 1$ and $n\gtrsim 1$, then for $\gamma \asymp {\left(\frac{{\sigma }^2{n}^{10/3}}{M^2}\right)}^{1/4}$ the solution $\hat{g}$ of (UCQP) satisfies (w.h.p.)

$${\Vert \frac{\hat{g}}{\left|\hat{g}\right|}-h\Vert }_2^2\lesssim (\sigma M+{\sigma }^2)n^{2/3}+\log n.$$

Hence for any $\epsilon \in (0,1)$, in the noise regime $\max \{\frac{M}{\epsilon n^{1/3}},\frac{\log n}{\sqrt{\epsilon n}}\}\lesssim \sigma \lesssim 1$, if $n\gtrsim {(1/\epsilon )}^3$, then (UCQP) denoises z w.h.p.

Results for (TRS). We now describe conditions under which the estimator (TRS) provably denoises $z\in {\mathcal{C}}_n$ w.h.p. Theorem 3 below is a simplified version of Theorem 8.

Theorem 3

For any $\epsilon \in (0,1)$, given $k\in [n-1]$ s.t. ${\lambda }_{n-k+1}<{\lambda }_{n-k}$ and $\bar{\lambda }\in [{\lambda }_{\min },{\lambda }_1]$ with the choice $\gamma \asymp {(\frac{{\sigma }^{2}n}{△B_n{\bar{\lambda }}^2})}^{1/4}$, suppose that the following conditions are satisfied.

• $B_n\lesssim \min \{n{\lambda }_{n-k},n\bar{\lambda }\}$, and $1+\left|{\mathcal{L}}_{\bar{\lambda }}\right|+\sqrt{(1+\left|{\mathcal{L}}_{\bar{\lambda }}\right|)\log n}\lesssim \epsilon n$.

• $\sigma \lesssim \min \{\sqrt{\epsilon },\frac{\bar{\lambda }}{{\lambda }_{n-k+1}^2}\sqrt{\frac{△B_n}{n}}\}$ and

  $$\sigma \gtrsim \max \{{\left(\frac{\log n}{\sqrt{n}}\right)}^{1/2},\frac{B_n}{n{\lambda }_{n-k}\sqrt{\epsilon }},\sqrt{\frac{\log n}{n\epsilon }},\frac{1}{\epsilon \bar{\lambda }}(\sqrt{\frac{△B_n}{n}}+{\lambda }_{n-k+1}^2\sqrt{\frac{n}{△B_n}})\}.$$

Then the (unique) solution $\hat{g}\in {\mathbb{C}}^n$ of (TRS) satisfies ${\Vert \frac{\hat{g}}{\left|\hat{g}\right|}-h\Vert }_2^2\le \epsilon {\Vert z-h\Vert }_2^2$ w.h.p.

The above statement is admittedly more convoluted than that for (UCQP) which is primarily due to the more intricate nature of the estimation error analysis. An interesting aspect of the above result is that we can consider the “best” choice of k (satisfying ${\lambda }_{n-k+1}<{\lambda }_{n-k}$) which leads to the mildest constraints on σ and $B_n$. Note that since G is connected (by assumption), $k=1$ always satisfies this condition ($0={\lambda }_n<{\lambda }_{n-1}$). This leads to the following useful Corollary which is a simplified version of Corollary 6.

Corollary 1

For any $\epsilon \in (0,1)$ and $\bar{\lambda }\in [{\lambda }_{\min },{\lambda }_1]$ with the choice $\gamma \asymp {(\frac{{\sigma }^{2}n}{△B_n{\bar{\lambda }}^2})}^{1/4}$, suppose that the following conditions are satisfied.

• $B_n\lesssim n{\lambda }_{\min }$ and $1+\left|{\mathcal{L}}_{\bar{\lambda }}\right|+\sqrt{(1+\left|{\mathcal{L}}_{\bar{\lambda }}\right|)\log n}\lesssim \epsilon n$.

• $\max \{{\left(\frac{\log n}{\sqrt{n}}\right)}^{1/2},\frac{B_n}{n{\lambda }_{\min }\sqrt{\epsilon }},\sqrt{\frac{\log n}{n}},\frac{1}{\epsilon \bar{\lambda }}\sqrt{\frac{△B_n}{n}}\}\lesssim \sigma \lesssim \sqrt{\epsilon }$.

Then the (unique) solution $\hat{g}\in {\mathbb{C}}^n$ of (TRS) satisfies ${\Vert \frac{\hat{g}}{\left|\hat{g}\right|}-h\Vert }_2^2\le \epsilon {\Vert z-h\Vert }_2^2$ w.h.p.

It is natural to ask whether Corollary 1 is – for all practical purposes – sufficient, compared to Theorem 3. For the complete graph $K_n$, observe that the only valid choice is $k=1$, hence we need Corollary 1. However, as we will see in Section 4, it turns out that for the path graph $P_n$, Corollary 1 leads to unreasonably strict conditions on σ and $B_n$ (see Remark 5). In fact for the path graph, ${\lambda }_{n-k+1}<{\lambda }_{n-k}$ holds for each $k\in [n-1]$, hence a careful choice of k in Theorem 3 is key to obtaining satisfactory conditions, see Corollary 8.

In the setting of the example from Section 1.2, Corollary 8 roughly states that for any $\epsilon \in (0,1)$, in the noise regime ${(\frac{\log n}{\sqrt{n}})}^{1/2}\lesssim \sigma \lesssim \sqrt{\epsilon }$ with n large enough and a suitably chosen γ, (TRS) denoises z w.h.p. This condition on σ is visibly stricter than that for (UCQP); we believe that this is likely an artefact of the analysis. Nevertheless, for this example, we see for ε fixed and $n\to \infty$ that both (TRS) and (UCQP) provably denoise z w.h.p. in the noise regime $o(1)\le \sigma \lesssim 1$.

Outline of the paper. The rest of the paper is organized as follows. Section 2 introduces some preliminaries involving certain intermediate facts and technical results that will be needed in our analysis. Section 3 contains the analysis for (UCQP) while Section 4 analyzes (TRS). Section 5 contains some simulation results on synthetic examples for (UCQP) and (TRS). We conclude with Section 6 which contains a discussion with related work from the literature, and directions for future work.