1 Introduction
The recovery of signals that lie on a manifold/surface has received extensive attention in the recent years. For example, patchbased image processing methods such as BM3D model the patches in an image as points on a manifold [1, 2], while we [3, 4] have recently used the manifold structure of images in a dynamic time series. Another area that witnessed extensive research in the recent years is the processing of signals on a graph [5]; these methods rely on the the graph Laplacian operator to denoise and postprocess signals on a graph.
The main focus of this paper is to introduce a continuous domain perspective on the recovery of points drawn from a smooth surface in very high dimensions. This work reveals fundamental links between recent advances in superresolution theory [6, 7, 8]
and kernel based machine learning methods
[9] as well as graph signal processing [10]. We assume that the high dimensional points live on an smooth surface, which is the zero level set of a bandlimited function. This is termed the annihilation relation and it is shown that this relation can be expressed as a weighted linear combination of the exponential features of the point; the dimension of the feature maps is equal to the bandwidth of the potential function. These properties enable us to determine the sampling conditions, which will guarantee the recovery of the surface from finite number of points. Our analysis also shows that when the bandwidth is overestimated, there are multiple such annihilation relations, suggesting that the exponential feature maps of the points on the surface live in a finite dimensional space. Note that similar nonlinear maps are widely used in kernel methods; our results show that these maps can be approximated by a few basis functions, when the points are restricted to a bandlimited surface.The finite dimensional nature of the maps translate to a lowrank kernel matrix, computed from the points using a shift invariant kernel such as the Dirichlet function. We minimize the nuclear norm of the feature maps of the points to recover them from noisy data. Since the direct estimation of the surface in higher dimensions suffers from the curse of dimensionality, we use the ”kernel trick” to keep the computational complexity manageable. We rely on an iterative reweighted algorithm to recover the denoised points. The resulting algorithm has similarities to iterative nonlocal methods
[11, 2, 12, 13, 5] that are widely used in image processing and graph signal processing. Specifically, it alternates between the estimation of a graph Laplacian, which specifies the connectivity of the points, and the smoothing of points guided by the graph Laplacian. Our experiments show that the Laplacian matrix obtained by solving the proposed optimization algorithm is more representative of the graph structure than classical methods, when it is estimated from noisy data.This work is built upon our prior work [14, 7, 15, 16, 17, 18] and the recent work by Ongie et al., which considered polynomial kernels [19]. Our main focus is to generalize [19] to shift invariant kernels, which are more widely used. We also introduce sampling conditions and algorithms to determine the surface, when the dimension is low. In addition, the iterative algorithm using the kernel trick shows the connections with graph Laplacian based methods used in graph signal processing.
2 Bandlimited surfaces & annihilation
We assume the point cloud to be supported on a surface in , which is the zero levelset of a bandlimited potential function:
(1) 
Here, is the smallest set of coefficients (minimal set) that satisfies the above relation. is a set of contiguous locations that indicates the support of the Fourier series coefficients of . Consider an arbitrary point on the above surface (1). By definition (1), we have the annihilation relation . We reexpress the annihilation relation as using a nonlinear mapping :
(2) 
This annihilation relation is illustrated in Fig 1.
2.1 Curve recovery: sampling conditions
The annihilation relation introduced in the previous subsection can be used to estimate the surface, or equivalently from a few number of points. The least square estimation of the coefficients from the data points can be posed as the minimization of the criterion:
(3) 
where . The coefficients can be estimated as:
(4) 
The solution is the minimum eigen vector of
.In the remainder of the section, we will restrict our attention to 2D for simplicity, even though the results in this section can be generalized to arbitrary dimensions. We will now determine the sampling conditions for the perfect recovery of the curve using (4). Specifically, we will determine the minimum number of samples for the successful recovery of the curve, when is a rectangular neighborhood in of size . In addition, we assume that is the function with the smallest Fourier support (minimal polynomial), whose zeros define the curve. We first focus on the case where is known.
Proposition 1.
Let be points on the zerolevel set of a bandlimited function , where the bandwidth of the surface is specified by and has irreducible factors. If points are sampled on the irreducible factor, then the curve can be uniquely recovered by (4), when:
(5) 
for .
Thus, the total number of points required are . We compare this setting with the sampling conditions for the recovery of a piecewise constant image, whose gradients vanish on a bandlimited curve [7]. The minimum number of Fourier measurements required to recover the function there is ; when , then complex Fourier samples are required. In contrast, we need real samples. When the true support is not known, it is a common practice to overestimate it as . In this case, will have multiple null space vectors, as shown below.
Proposition 2.
We consider the polynomial described in Proposition 1. Let with and for :
(6) 
points be sampled on the irreducible factor of . Then all nullspace vectors of the matrix will be of the form:
(7) 
where is an arbitrary function such that .
Thus, the total number of points required are . Since is the common factor of all the annihilating functions, all of them will satisfy , for any point on the original curve. Depending on the specific , they will have additional zeros. Hence, the above result provides us a means to compute the original curve, even when the original bandwidth/support of the function is unknown.
We now consider a collection of points on the curve, stacked into a matrix . Let the feature matrix of size be denoted by:
(8) 
We state a result about the rank of the above feature matrix.
Proposition 3.
We consider the polynomial described in Proposition 1 and . Then:
(9) 
with equality if the sampling conditions of Proposition 2 are satisfied.
Here, denotes the number of valid shifts of the set within as shown in Fig 2 (a). Note that as gets smaller, the number of shifts of it within increases, and hence the rank decreases. The rank of the matrix can be used as a surrogate for the bandwidth of , or equivalently the complexity of the curve. Note that may be an irregular shape in . For example, if the points lie on a line in , then could be concentrated along a line in , resulting in a small , even when the number of features in may be considerably high. The lowrank structure of the feature maps can be used to denoise the original points, while the sum of squares function obtained from the nullspace filters can be used to estimate the surface in lowdimensions when (6) is satisfied, as illustrated in Fig 2.
2.2 Recovery of noisy point clouds in high dimensions
The explicit approach of estimating the surface is feasible, when the dimension of the points is small. However, this approach suffers from the curse of dimensionality. Since the shape of the data, or equivalently the shape of the support is not known, one needs to use a large to ensure that . Note that the dimension of the feature space specified by grows exponentially with , making this approach impractical in applications involving point clouds of images or patches.
We hence rely on the right nullspace relations to recover the points from their noisy and undersampled measurements. Specifically, we are interested in the null space relations
(10) 
where the entries of the Gram matrix are
(11) 
The function in (11) is shift invariant and is dependent on the shape of . For example, when is a centered cube in , is a Dirichlet function. The kernel matrix satisfies , where is given by (9).
2.3 Dirichlet and Gaussian surface representation
The bandlimited function in (1) can equivalently be expressed as:
(12) 
where is the Dirichlet function dependent on and is the set of sampled locations on the curve. Using reciprocity, the nonlinear maps in this case can be shown to be:
(13) 
Since the implicit curve is the zero level set of a linear combination of Dirichlet functions, it may be highly oscillatory. An alternative would be to use a level set expansion in terms of weighted exponentials , which could give smoother surfaces. In this case approaches a periodized Gaussian function, as , and the Gaussian kernel matrix is theoretically full rank. However, we observe that the Fourier series coefficients of a Gaussian function can be safely approximated to be zero outside , which translates to ; i.e., the rank will be small for high values of . We choose Gaussian kernels since they are more isotropic and less oscillatory than the Dirichlet kernel.
2.4 Denoising using nuclear norm minimization
We rely on the low rank structure of the kernel matrix to recover the noisy points. Specifically, with the addition of noise, the points deviate from the zero set of . A high bandwidth potential function is needed to represent the noisy surface. We propose to use the nuclear norm of the feature matrix as a regularizer in the recovery of the points from noisy measurements:
(14) 
We use the IRLS algorithm, where is updated as:
(15) 
and . Note that the solution for (15) involves a system of nonlinear equations. Instead, we use gradient linearization to simplify our computations, where is a Gaussian kernel matrix:
(16) 
with , , and
(17) 
We note the equivalence of the above optimization strategy with widely used nonlocal means and graph optimization schemes. These schemes estimate a Laplacian matrix , followed by the minimization of the cost function (16). These approaches can thus be seen as fitting a smooth bandlimited surface to the point cloud of patches or signals that are assumed to be on the graph.
3 Results
We demonstrate the utility of (14) in a simple 2D denoising example in Fig 3. Specifically, we consider the recovery of points on the TigerHawk logo from its noisy samples. See the caption for details. The top row shows that the proposed algorithm is able to provide good denoising of the data. The bottom three rows show that the Laplacian estimated using (17) at the iteration is more representative of the shape.
The utility of the proposed method in denoising free breathing and ungated MRI data is shown in Fig 4. Since MRI is a slow imaging modality, several rapid imaging techniques were introduced to accelerate the acquisition. All of these methods trade SNR for speed, resulting in noisy images. The proposed scheme is seen to exploit the manifold structure of the data to reduce noise. See caption for details.
4 Conclusion
We introduce a continuous domain framework for the recovery of points on a bandlimited surface. We show that the exponential maps of the points lie in a lower dimensional subspace, which translates to a kernel matrix that is lowrank. We introduce a nuclear norm minimization algorithm to recover the points. The proposed framework connects the continuous domain surface recovery problem with kernel methods and approaches in graph signal processing. The application of the algorithms to noisy data reveals its great utility in practical problems.
References
 [1] A. Danielyan, V. Katkovnik, and K. Egiazarian, “Bm3d frames and variational image deblurring,” IEEE Tran. Image Processing, vol. 21, no. 4, pp. 1715–1728, April 2012.
 [2] Yasir Q. Mohsin, Sajan Goud Lingala, Edward DiBella, and Mathews Jacob, “Accelerated dynamic mri using patch regularization for implicit motion compensation,” Magnetic Resonance in Medicine, vol. 77, no. 3, pp. 1238–1248, 2017.
 [3] S. Poddar and M. Jacob, “Dynamic mri using smoothness regularization on manifolds (storm),” IEEE Tran. Medical Imaging, vol. 35, no. 4, pp. 1106–1115, April 2016.
 [4] Sunrita Poddar, Sajan Goud Lingala, and Mathews Jacob, “Joint recovery of under sampled signals on a manifold: Application to free breathing cardiac mri,” in Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on. IEEE, 2014, pp. 6904–6908.

[5]
Vassilis Kalofolias,
“How to learn a graph from smooth signals,”
in
Proc. Int. Conf. Artificial Intelligence and Statistics
, 2016, vol. 51, pp. 920–929. 
[6]
Emmanuel J Candès and Carlos FernandezGranda,
“Towards a mathematical theory of superresolution,”
Communications on Pure and Applied Mathematics, vol. 67, no. 6, pp. 906–956, 2014.  [7] Greg Ongie and Mathews Jacob, “Offthegrid recovery of piecewise constant images from few fourier samples,” SIAM Journal on Imaging Sciences, vol. 9, no. 3, pp. 1004–1041, 2016.
 [8] Geoffrey Schiebinger, Elina Robeva, and Benjamin Recht, “Superresolution without separation,” in Computational Advances in MultiSensor Adaptive Processing (CAMSAP), 2015 IEEE 6th International Workshop on. IEEE, 2015, pp. 45–48.

[9]
Bernhard Schölkopf and Alexander J Smola,
Learning with kernels: support vector machines, regularization, optimization, and beyond
, MIT press, 2002. 
[10]
David I Shuman, Sunil K Narang, Pascal Frossard, Antonio Ortega, and Pierre
Vandergheynst,
“The emerging field of signal processing on graphs: Extending highdimensional data analysis to networks and other irregular domains,”
IEEE Signal Processing Magazine, vol. 30, no. 3, pp. 83–98, 2013.  [11] Z. Yang and M. Jacob, “Nonlocal regularization of inverse problems: A unified variational framework,” IEEE Transactions on Image Processing, vol. 22, no. 8, pp. 3192–3203, Aug 2013.
 [12] Yasir Q Mohsin, Gregory Ongie, and Mathews Jacob, “Iterative shrinkage algorithm for patchsmoothness regularized medical image recovery,” IEEE transactions on medical imaging, vol. 34, no. 12, pp. 2417–2428, 2015.
 [13] Guy Gilboa and Stanley Osher, “Nonlocal operators with applications to image processing,” Multiscale Modeling & Simulation, vol. 7, no. 3, pp. 1005–1028, 2008.
 [14] Greg Ongie and Mathews Jacob, “Recovery of piecewise smooth images from few fourier samples,” in Sampling Theory and Applications (SampTA), 2015 International Conference on. IEEE, 2015, pp. 543–547.
 [15] Greg Ongie, Sampurna Biswas, and Mathews Jacob, “Convex recovery of continuous domain piecewise constant images from nonuniform fourier samples,” IEEE Transactions on Signal Processing, vol. 66, no. 1, pp. 236–250, 2017.
 [16] Arvind Balachandrasekaran, Vincent Magnotta, and Mathews Jacob, “Recovery of damped exponentials using structured low rank matrix completion,” IEEE transactions on medical imaging, vol. 36, no. 10, pp. 2087–2098, 2017.
 [17] Arvind Balachandrasekaran and Mathews Jacob, “Novel structured lowrank algorithm to recover spatially smooth exponential image time series,” in Biomedical Imaging (ISBI 2017), 2017 IEEE 14th International Symposium on. IEEE, 2017, pp. 1–4.
 [18] Gregory Ongie and Mathews Jacob, “A fast algorithm for convolutional structured lowrank matrix recovery,” IEEE Transactions on Computational Imaging, 2017.
 [19] Greg Ongie, Rebecca Willett, Robert D. Nowak, and Laura Balzano, “Algebraic variety models for highrank matrix completion,” in Proceedings of the 34th International Conference on Machine Learning, 2017, pp. 2691–2700.
Comments
There are no comments yet.