Error analysis for denoising smooth modulo signals on a graph
Introduction
Many modern applications involve the acquisition of noisy modulo samples of a function , i.e., we obtain for some , where denotes noise. Here, lies in the interval and is such that 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 , 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 from now. Given (1.1), one can ask whether we can recover the original samples of f, i.e., ? 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 , it is not difficult to see that we can exactly recover the samples of f when there is no noise (i.e., in (1.1)) provided the sampling is fine enough. This is achieved by sequentially traversing the samples and reconstructing the quotient corresponding to by setting it to where equals either (a) 0 (if is “close enough” to ), or (b) ±1 (if is sufficiently smaller/larger than ). If f is smooth, then a fine enough sampling density will ensure that nearby samples of f have quotients which differ by either or 0. This argument can be extended to the general multivariate setting as well (see [12]).
While exact recovery of is of course no longer possible in the presence of noise, one can show [12] that if is not too large (uniformly for all i), and if n is sufficiently large, then the estimates returned by the above procedure are such that (up to a global integer shift) for each i, is bounded by a term proportional to the noise level. This suggests a natural two stage procedure for estimating – first obtain denoised estimates of , and in the second stage, “unwrap” these denoised modulo samples to obtain the estimates .
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., ), with the exception of for . The imaginary unit is denoted by . For , we write when there exists an absolute constant such that . If and , then we write . For , we write where denote its real and imaginary parts respectively. The Hermitian conjugate of u is denoted by , and denotes the usual norm in for .
Cucuringu and Tyagi [11], [10] considered the problem of denoising modulo samples by formulating it as an optimization problem on a manifold. Assume , and suppose that 's form a uniform grid in . The main idea is to represent each sample via and observing that if is close to , then holds provided f is smooth. In other words, the samples lie on the manifold which is the product manifold of unit radius circles, i.e., . Hence, by constructing a proximity graph on the sampling points – where there is an edge between i and j if is within a specified distance of – they proposed solving a quadratically constrained quadratic program (QCQP) where L is the Laplacian of G, and is a regularization parameter. While the objective function is convex, this is a non-convex problem due to the constraint set . It is not clear whether the global solution of (QCQP) can be obtained in polynomial time. Hence they proposed relaxing to a sphere leading to the trust region subproblem (TRS)
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 stability result was shown for the solution of (TRS). To describe this result, denote to be the ground truth samples where , and define the entry-wise projection operator [21] Consider the univariate case with , and assume that the (Gaussian) noise i.i.d. for . Assuming f is M-Lipschitz, and , denote Δ to be the maximum degree of G with 's forming a uniform grid on . Then provided , the bounds in [11, Theorem 14] imply1 that with high probability, the solution of (TRS) satisfies On the other hand, one can show that with high probability if , hence we cannot conclude 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 is smooth w.r.t. G. This is formally described below, followed by a summary of our main results.
Consider an undirected, connected graph where denotes its set of edges. Denote the maximum degree of G by △, and the (combinatorial) Laplacian matrix associated with G by . We will denote the eigenvalues of L by where denotes the Fiedler value of G, and is well known to be a measure of connectivity of G. The corresponding (unit norm) eigenvectors of L will be denoted by ; where we have that .
Let be an unknown ground truth signal which is assumed to be smooth w.r.t. G in the sense that Here, 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 . Specifically, denoting ; to be i.i.d. Gaussian random variables, we are given Given access to the noisy samples , and the graph G, we aim to answer the following two questions.
- 1.
Consider the unconstrained quadratic program (UCQP) obtained by relaxing the constraints in (QCQP) to 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 . Under what conditions can we ensure that satisfies (1.4)?
- 2.
Under what conditions can we ensure that the solution of (TRS) satisfies (1.4)?
Remark 1
It is worth mentioning that the quadratic penalty term in (UCQP) and (TRS) aims to promote solutions which are smooth w.r.t G. This is natural given that the ground truth signal 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 , which is equivalent to removing the regularization term, then we obtain . The main point of the ensuing analysis is to find a suitable “intermediate” choice of γ which ensures that 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 to be a M-Lipschitz function where for each , let for be a uniform grid. For each i, we are given where are i.i.d. Then denoting and , we will consider to be the path graph (), where
Before stating our results, let us define for any , the set 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 and satisfying , suppose that . Then for the choice , the solution of (UCQP) satisfies .
We also derive conditions under which denoising occurs with high probability3 (w.h.p.). Theorem 2 below is a simplified version of Theorem 6. Theorem 2 For any given and , suppose that If , then w.h.p., the solution of (UCQP) satisfies .
In Section 3, we interpret our results for special graphs4 such as and the star graph5 (); see Corollary 2, Corollary 3. In particular, the statement for therein can be applied to the example from Section 1.2, readily yielding conditions that ensure denoising; see Corollary 4. It states that if and , then for the solution of (UCQP) satisfies (w.h.p.) Hence for any , in the noise regime , if , then (UCQP) denoises z w.h.p.
Results for (TRS). We now describe conditions under which the estimator (TRS) provably denoises w.h.p. Theorem 3 below is a simplified version of Theorem 8. Theorem 3 For any , given s.t. and with the choice , suppose that the following conditions are satisfied. , and . and
Then the (unique) solution of (TRS) satisfies w.h.p.
Corollary 1
For any and with the choice , suppose that the following conditions are satisfied.
- (i)
and .
- (ii)
.
In the setting of the example from Section 1.2, Corollary 8 roughly states that for any , in the noise regime 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 that both (TRS) and (UCQP) provably denoise z w.h.p. in the noise regime .
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.
Section snippets
Preliminaries
This section summarizes some useful technical tools that will be employed at multiple points in our analysis. We begin by deriving some simple consequences of the smoothness assumption in (1.5).
Smooth phase signals. Similar to (1.8), for , let us define the set consisting of indices corresponding to the “high frequency” part of the spectrum of L. Then using (1.5) and the fact , it is not difficult to establish that
Error bounds for (UCQP)
We now analyse the quality of the solution of (UCQP). To begin with, we have the following Lemma that bounds the error between and h for any . Lemma 1 For any and , it holds that where Proof From the definition in (1.2), observe that where we recall . We can write
Error bounds for (TRS)
We now proceed to derive a error bound for the solution of (TRS). Following the notation in [11], we can represent any via , where . Moreover, consider the matrix where ⊗ denotes the Kronecker product. Then it is easy to check that (TRS) is equivalent to where is a solution of (TRS) iff is a solution of ().
The following Lemma from [11], which in turn is a direct consequence
Simulations
We now provide simulation results on some synthetic examples. For concreteness, we consider the following functions.
- 1.
,
- 2.
.
Discussion
We now discuss in detail some related work and conclude with directions for future work.
Acknowledgements
I would like to thank Stéphane Chrétien and Michaël Fanuel for carefully reading a preliminary version of the draft, and for providing useful feedback; Alain Celisse for the very helpful technical discussions during the early stages of this work.
References (36)
- et al.
Automatic interferogram analysis techniques applied to quasi-heterodyne holography and ESPI
Opt. Lasers Eng.
(1991) - et al.
Solving the trust-region subproblem by a generalized eigenvalue problem
SIAM J. Optim.
(2017) - et al.
Regularization and semi-supervised learning on large graphs
Concentration of Quadratic Forms Under a Bernstein Moment Assumption
(2019)- et al.
On unlimited sampling
- et al.
Unlimited sampling of sparse signals
- et al.
Unlimited sampling of sparse sinusoidal mixtures
- et al.
Concentration Inequalities: A Non Asymptotic Theory of Independence
(2013) - et al.
Spectra of Graphs
(2012) - et al.
Understanding phase maps in MRI: a new cutline phase unwrapping method
IEEE Trans. Med. Imaging
(2002)
On denoising modulo 1 samples of a function
Provably robust estimation of modulo 1 samples of a smooth function with applications to phase unwrapping
J. Mach. Learn. Res.
Denoising modulo samples: k-NN regression and tightness of SDP relaxation
Inf. Inference: J. IMA
Synthetic interferometer radar for topographic mapping
Proc. IEEE
Minimizing a quadratic over a sphere
SIAM J. Optim.
A new two-dimensional phase unwrapping algorithm for MRI images
Magn. Reson. Med.
Path-independent phase unwrapping using phase gradient and total-variation (TV) denoising
Opt. Express
Mt-025 tutorial ADC architectures VI: Folding ADCs. Analog Devices
Cited by (3)
Unlimited Sampling in Phase Space
2023, ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings