Solving Partial Differential Equations on Manifolds From Incomplete Inter-Point Distance

TL;DR

This paper introduces an efficient method for solving PDEs on manifolds using incomplete inter-point distance data, leveraging local reconstruction and low-rank matrix completion without full manifold reconstruction.

Contribution

It proposes a novel framework that discretizes PDEs directly from incomplete distance data using local reconstruction and low-rank matrix completion, avoiding global coordinate reconstruction.

Findings

Effective PDE solutions on incomplete distance data.

Successful manifold reconstruction from partial distances.

Numerical experiments validate the approach.

Abstract

Solutions of partial differential equations (PDEs) on manifolds have provided important applications in different fields in science and engineering. Existing methods are majorly based on discretization of manifolds as implicit functions, triangle meshes, or point clouds, where the manifold structure is approximated by either zero level set of an implicit function or a set of points. In many applications, manifolds might be only provided as an inter-point distance matrix with possible missing values. This paper discusses a framework to discretize PDEs on manifolds represented as incomplete inter-point distance information. Without conducting a time-consuming global coordinates reconstruction, we propose a more efficient strategy by discretizing differential operators only based on point-wisely local reconstruction. Our local reconstruction model is based on the recent advances of low-rank matrix completion theory, where only a very small random portion of distance information is required. This method enables us to conduct analyses of incomplete distance data using solutions of special designed PDEs such as the Laplace-Beltrami (LB) eigen-system. As an application, we demonstrate a new way of manifold reconstruction from an incomplete distance by stitching patches using the spectrum of the LB operator. Intensive numerical experiments demonstrate the effectiveness of the proposed methods.

Tables

Table 1. Table 1: Rate of The successful reconstruction ρ 𝜌 \rho and the average relative error E B subscript 𝐸 𝐵 E_{B} of the Gram matrix out of 50 50 50 tests by the proposed model ( 2 ) from distances with information availability rate γ 𝛾 \gamma .

		1%	2%	3%	5%	10%	20%
$S^2$	$E_B$	7.157E-1	1.376E-3	4.791E-4	2.474E-4	1.342E-5	4.262E-5
	$\rho$	0%	92%	100%	100%	100%	100%
Cow	$E_B$	4.9427E-5	3.980E-4	1.837E-4	5.319E-5	1.4072E-5	2.155E-5
	$\rho$	100%	100%	100%	100%	100%	100%
Swiss roll	$E_B$	2.722E-4	2.894E-4	1.633E-4	5.054E-5	1.704E-5	1.114E-5
	$\rho$	100%	100%	100%	100%	100%	100%

Table 2. Table 2: Comparisons of time consumption (minutes) of solving the LB eigenvalue problem based on local/global reconstruction methods.

number of points
1002	1962	4002	7842	16002
$\gamma =100\%,\ell =6$
0.26	0.51	1.01	2.03	4.05
$\gamma =80\%,\ell =9$
2.28	5.60	11.17	22.28	45.02
$\gamma =50\%,\ell =18$
4.03	8.09	16.14	32.44	64.71
$\gamma =30\%,\ell =30$
15.13	30.19	60.42	120.95	241.63
global reconstruction using 3% distance ($\ell =6$ for MLS)
2.09	9.86	40.13	154.40	597.06

Table 3. Table 3: Time comparisons (seconds) between global reconstruction and the stitching model.

	Armadillo	Kitten	Swiss roll
Global reconstruction	35321.84	1315.50	622.65
Stitching	760.18	138.76	199.75

Equations

$$\underset{B\in \mathbb{S}}{\min }\mathrm{Rank}(B),\text{s.t.}B⪰0\&{\mathcal{R}}_{\Omega }(B)=D_{\Omega }$$

$$\underset{B\in \mathbb{S}}{\min }\Vert B{\Vert }_{\ast },\text{s.t.}B⪰0,{\mathcal{R}}_{\Omega }(B)=D_{\Omega }\&B1=0$$

$$\mathcal{A}:{\mathbb{R}}^{l(l+1)/2}\to {\mathbb{R}}^{\Omega }\times {\mathbb{R}}^{\ell },\hat{B}\mapsto \left({\mathcal{R}}_{\Omega }\circ \iota (\hat{B}),\iota (\hat{B})1\right)$$

$$\underset{\hat{B}\in {\mathbb{R}}^{k(k+1)/2}}{\min }\text{Trace}(\iota (\hat{B})),\text{s.t.}\mathcal{A}\hat{B}={\tilde{D}}_{\Omega },\iota (\hat{B})⪰0$$

$$\underset{\hat{B},\hat{C}\in {\mathbb{R}}^{n(n+1)/2}}{\min }\text{Trace}(\iota (\hat{B})),\text{s.t.}\mathcal{A}\hat{B}={\tilde{D}}_{\Omega },\iota (\hat{C})⪰0,\hat{B}=\hat{C}$$

$${\mathcal{L}}_{{\mu }_1,{\mu }_2}(\hat{B},\hat{C};H_1,H_2)=\text{Trace}(\iota (\hat{B}))+\frac{{\mu }_1}{2}\Vert \mathcal{A}\hat{B}-{\tilde{D}}_{\Omega }+H_1{\Vert }_2^2+\frac{{\mu }_2}{2}\Vert \hat{B}-\hat{C}+H_2{\Vert }_2^2.$$

$$\{\begin{array}{l}{\hat{B}}^{k+1}=\arg \underset{\hat{B}}{\min }\text{Trace}(\iota (\hat{B}))+\frac{{\mu }_1}{2}\Vert \mathcal{A}\hat{B}-{\tilde{D}}_{\Omega }+H_1{\Vert }_2^2+\frac{{\mu }_2}{2}\Vert \hat{B}-\hat{C}+H_2{\Vert }_2^2\\ {\hat{C}}^{k+1}=\arg \underset{\hat{C}}{\min }\frac{{\mu }_2}{2}\Vert {\hat{B}}^{k+1}-\hat{C}+H_2^k{\Vert }^2,\iota (\hat{C})⪰0,\\ H_1^{k+1}=H_1^k+(\mathcal{A}\hat{B}-{\tilde{D}}_{\Omega }),\\ H_2^{k+1}=H_2^k+\hat{B}-\hat{C},\end{array}$$

$$({\mu }_1{\mathcal{A}}^{\ast }\mathcal{A}+{\mu }_2){\hat{B}}^{k+1}={\mu }_1{\mathcal{A}}^{\ast }({\tilde{D}}_{\Omega }-H_1^k)+{\mu }_2({\hat{C}}^k-H_2^k)-\zeta (I_n).$$

$$U_{ij}=\{\begin{array}{ll}\sum _{(i,k)\in \Omega }{u}_{ik}+\sum _{(k,i)\in \Omega }{u}_{ki},& \text{if}i=j\\ -2u_{ij},& \text{if}i>j\end{array}\text{and}{Q}_{ij}=\{\begin{array}{ll}{q}_i,& \text{if}i=j\\ q_i+q_j,& \text{if}i>j\end{array}$$

$${\mathcal{T}}_E(X)=V{\Lambda }_{+}{V}^{\top },{\Lambda }_{+}=\text{diag}(\max \{\Lambda (i,i),0{\}}_{i=1}^n).$$

$${\hat{C}}^{k+1}=\zeta \circ {\mathcal{T}}_E\circ \iota ({\hat{B}}^{k+1}+H_2^k).$$

$$\underset{\tilde{D}\in \mathbb{S}}{\min }\Vert \tilde{D}{\Vert }_{\ast }\text{s.t.}{\tilde{D}}_{ij}=D_{ij},(i,j)\in \Omega$$

$$P_i=\sum _{k\in N(i)}({\mathbf{x}}_k-c_i{)}^T({\mathbf{x}}_k-c_i).$$

$$\underset{{\mathbf{c}}_{\mathbf{\xi }}^1,\cdots ,{\mathbf{c}}_{\mathbf{\xi }}^{p-d}}{\min }\sum _{k\in N(i)}w(\Vert {\mathbf{x}}_k-{\mathbf{x}}_i\Vert )\Vert {\mathbf{Z}}_i({\mathbf{u}}_k)-{\mathbf{v}}_k{\Vert }^2,$$

$${\nabla }_{\mathcal{M}}f=g^{st}\frac{\partial f}{\partial u_t}{\partial }_{u_s},{\mathrm{div}}_{\mathcal{M}}\mathbf{V}=\frac{1}{\sqrt{g}}\frac{\partial }{\partial u_s}\left(\sqrt{g}{V}^s\right),{\Delta }_{\mathcal{M}}f=\frac{1}{\sqrt{g}}\frac{\partial }{\partial u_s}\left(\sqrt{g}{g}^{st}\frac{\partial f}{\partial u_t}\right).$$

$${\nabla }_{M}f({\mathbf{x}}_i)=\frac{{\sum }_{S\in \mathcal{R}(i)}|S|{\nabla }_{S}f({\mathbf{x}}_i)}{{\sum }_{S\in \mathcal{R}(i)}|S|},{\mathrm{div}}_M\overrightarrow{V}(p_i)=\frac{{\sum }_{S\in \mathcal{R}(i)}|S|{\mathrm{div}}_S\overrightarrow{V}({\mathbf{x}}_i)}{{\sum }_{S\in \mathcal{R}(i)}|S|}$$

$${\Delta }_{\mathcal{M}}{\varphi }_k=-{\lambda }_k{\varphi }_k,k=0,1,2,\cdots .$$

$$\{\begin{array}{l}|{\nabla }_{\mathcal{M}}d(x)|=1\\ d(x)=0,x\in \Gamma \subset \mathcal{M}\end{array}$$

$$\underset{P,\{R_j\in O(d),b_j\in {\mathbb{R}}^d{\}}_{j=1}^L}{\min }\sum _{j=1}^L\Vert P_{{\Omega }_j}-Q_j{R}_j-1_j{b}_j{\Vert }_2^2,$$

$$\underset{\mathbf{\alpha }}{\min }\Vert P-\sum _{i=1}^N{\varphi }_i{\alpha }_i{\Vert }_F^2=\underset{\mathbf{\alpha }}{\min }\Vert P-\Phi \mathbf{\alpha }{\Vert }_F^2,$$

$$\underset{\mathbf{\alpha }\in {\mathbb{R}}^{N\times d},\{R_j\in O(d),b_j\in {\mathbb{R}}^d{\}}_{j=1}^L}{\min }\mathcal{E}(\mathbf{\alpha },R_j,b_j)=\frac{1}{2}\sum _{j=1}^L\Vert Q_j-{\Phi }_{{\Omega }_j}\mathbf{\alpha }{R}_j-1_j{b}_j{\Vert }_F^2,$$

$$\{\begin{array}{l}{\mathbf{\alpha }}^{k+1}={\mathbf{\alpha }}^k-(\sum _{j=1}^L{\Phi }_{{\Omega }_j}^{\top }{\Phi }_{{\Omega }_j}{)}^{-1}{\nabla }_{\mathbf{\alpha }}\mathcal{E}({\mathbf{\alpha }}^k,R_j^k,b_j^k),\\ R_j^{k+1}=\mathcal{K}(R_j^k,{\nabla }_{R_j}\mathcal{E}({\mathbf{\alpha }}^{k+1},R_j^k,b_j^k)),\forall 1\le j\le L,\\ b_j^{k+1}=b_j^k-\frac{1}{\text{Card}({\Omega }_j)}{\nabla }_{b_j}\mathcal{E}({\mathbf{\alpha }}^{k+1},R_j^{k+1},b_j^k),\forall 1\le j\le L.\end{array}$$

$$\{\begin{array}{l}{\nabla }_{\mathbf{\alpha }}\mathcal{E}({\mathbf{\alpha }}^k,R_j^k,b_j^k)=-\sum _{j=1}^L{\Phi }_{{\Omega }_j}^{\top }(Q_j-{\Phi }_{{\Omega }_j}{\mathbf{\alpha }}^k{R}_j^k-1_j{b}_j^k)(R_j^k{)}^{\top },\\ {\nabla }_{R_j}\mathcal{E}({\mathbf{\alpha }}^{k+1},R_j^k,b_j^k)=-({\mathbf{\alpha }}^k{)}^{\top }{\Phi }_{{\Omega }_j}^{\top }(Q_j-{\Phi }_{{\Omega }_j}{\mathbf{\alpha }}^k{R}_j^k-1_j{b}_j^k),\\ {\nabla }_{b_j}\mathcal{E}({\mathbf{\alpha }}^{k+1},R_j^{k+1},b_j^k)=-1_j^{\top }(Q_j-{\Phi }_{{\Omega }_j}{\mathbf{\alpha }}^k{R}_j^k-1_j{b}_j^k).\end{array}$$

$$G_j^k$$

$$S_j^k$$

$$R_j^{k+1}=\mathcal{K}(R_j^k,G_j^k)=(1+\frac{\delta }{2}{S}_j^k{)}^{-1}(1-\frac{\delta }{2}{S}_j^k)R_j^k,$$

$$D_{\gamma,\ell}=\Big{\{}D(i,j)~{}|~{}j\in N_{\ell}(i)~{}\&~{}(i,j)\in\Omega_{\gamma},\quad i=1,\cdots,n\Big{\}}$$

$$S=\{(t+0.1)\cos (t),(t+0.1)\sin (t),8\pi s)|t\in [0,4\pi ],s\in [0,1]\}$$

Full text

{AMS}

65D18, 65D25, 65N25

1 Introduction

With the rapid development of advanced data acquisition technology, processing and analyzing data sampled on 3D shapes or even higher dimensional geometric objects becomes ubiquitous tasks such as those used in a 3D camera, medical imaging, protein structuring, social network analysis and many others [21, 2, 58, 50, 17, 6, 44, 22, 9]. PDE and variational PDE based methods have made great success to handle problems in signal and image processing which can be viewed as data on Euclidean domains. It is natural to consider PDE based methods to analyze and process signals on a general manifold and to understand geometric structures hidden in the data. Besides the classical implicit methods [38, 8, 7], finite difference methods [39, 52, 36, 62], finite element methods [41, 27, 20] and parameterization methods [49, 48, 57, 35] for solving differential equations on surfaces in $\mathbb{R}^{3}$, there has been increasing interests of solving PDEs on general $d$-dimensional manifold in $\mathbb{R}^{p}$ and their applications to data analysis. For instance, a diffusion geometry framework is developed to investigate the geometric structure of data based on solving Laplace-Beltrami (LB) eigenproblem using integral kernel methods [3, 4, 16]. More recently, a moving least square method and a local mesh method are considered to intrinsically solve different types of PDEs on manifolds represented as point clouds and its applications to the geometric understanding of point cloud data [29, 33, 34, 31].

All the aforementioned methods of solving PDEs on a general manifold $\mathcal{M}\subset\mathbb{R}^{p}$ are commonly considered $\mathcal{M}$ is sampled as a set of points $\{\bm{x}_{i}\in\mathbb{R}^{p}\}_{i=1}^{n}$, referred as a point cloud, and discretization of differential operators or approximation of integral equations are relied on available coordinates information of $\{\bm{x}_{i}\}_{i=1}^{n}$. However, there are many applications whose input information has no point coordinates but only an incomplete inter-point distance $(d(\bm{x}_{i},\bm{x}_{j}))$. Examples include incomplete distance information from sensor network localization [22, 9], protein structuring from NMR spectroscopy [17, 6] and global position from local distance configuration of cities [47]. A well-known distance geometry problem [17, 37] is to find the global configuration of data based on the incomplete distance information. The objective in this paper is different from the canonical distance geometry problem. We would like to develop numerical methods for solving PDEs on manifolds represented as incomplete inter-point distance. One natural way to approach this problem is first to apply a global reconstruction algorithm to obtain a point cloud representation of the input manifold, then established numerical methods for PDEs on point clouds can be directly applied. However, the global coordinate reconstruction might be very time consuming as it involves with semi-positive definite programming whose size depends on the number of points and could be very large in practice [9].

In this paper, we propose a different strategy to solve PDEs on manifolds represented as incomplete inter-point distance without conducting global coordinates reconstruction. Our idea is based on two intrinsic properties of differential operators on manifolds. Namely, the definition of a differential operator is only point-wisely depending on local information of the manifold and is invariant under different choices of local coordinates. This motivates us to only consider to conduct point-wisely local coordinate reconstruction for the associated local neighborhood, then we can point-wisely approximate differential operators based on the reconstructed local coordinates. After that, numerical solver for differential equations can be constructed.

Inspired by the classical multi-dimensional scaling [26], the full distance matrix is one-to-one corresponding to the Gram matrix which can be further used to determine coordinates by its eigen-decomposition. More importantly, the Gram matrix can be essentially viewed as an inner product matrix after certain coordinates shift. As long as the matrix size, namely the size of a given local neighborhood, is larger than the embedding dimension of the local neighborhood, then the corresponding Gram matrix has to be a low-rank matrix. This suggests us to use low rank as a prior knowledge to reconstruct the Gram matrix based on constraints of available distance information. Adapted from the recent advances of the low-rank matrix completion theory [13], we consider a nuclear norm regularized convex optimization problem to reconstruct local coordinates based on available distance information. Once local coordinates can be obtained, we apply the intrinsic methods of approximating differential operators on point clouds developed in [29, 33, 34] to discretize the desired differential equation. These intrinsic methods can be used to discretize different types of PDEs including parabolic, elliptic and hyperbolic PDEs. Our method can be viewed as natural extensions of these two methods to a new data structure where no coordinate information but only partial inter-point distance information of point clouds is provided. The outline of our strategy can be as summarized as follows:

Step 0

For the $i$-th point, chose its K-nearest neighborhood (KNN) $N(i)$ based on the given incomplete distance or from the prior information if its available.

Step 1

Applying the matrix completion method discussed in Section 2 to reconstruct local coordinates for KNN of the $i$-th point.

Step 2

Applying the MLS method or the local mesh method in [29, 33, 34] to approximate the desired differential operators at the $i$-th point. This provides the $i$-th row of the discretized matrix representation of the desired differential equation.

An immediate advantage of this approach is to save computation time by avoiding the global coordinate reconstruction, which can reduce the complexity quadratically to linearly scaling with the total number of points. More details about this method will be discussed in Section 3 and will be verified in our numerical experiments. Furthermore, this approach enables us to conduct geometric understanding of data without global coordinates reconstruction. Examples include global pattern extraction, comparison and classification as many existing methods conducted using results of differential equations [46, 30, 31]. As a byproduct of global information from PDEs, we also propose a new method for reconstructing manifolds by stitching its local patches. This new method is much more efficient than the way of direct global reconstruction using matrix completion method. Moreover, it can also overcome possible reconstruction failure using global matrix completion due to the coherent missing information.

The rest of this paper is organized as follows. In section 2, we propose a low-rank matrix completion model to reconstruct local coordinates and design an algorithm to solve the proposed convex optimization problem based on operator splitting and the alternating direction method. After that, section 3 is devoted to discuss point-wisely approximating differential operators based on the intrinsic methods proposed in [29, 33, 34]. We also discuss our two model PDEs, the Laplace-Beltrami (LB) eigenvalue problem and the Eiknoal equation for solve geodesic distance on manifolds. As an application of global information using solutions of LB eigenvalue problem, we propose a new method for reconstructing manifolds by stitching local patches in the LB frequency space in section 4. All numerical experiments about the effectiveness of our local coordinates reconstruction model, the accuracy and robustness of the PDEs solvers, and demonstration of manifold reconstruction are discussed in section 5. Finally, we conclude our work in section 6.

2 Local coordinates reconstruction via matrix completion

In this section, we first review the concept of the classical multidimensional scaling (MDS). Inspired by the classical MDS, we propose a method of reconstructing local coordinates from the given incomplete inter-point distance matrix using the low-rank matrix completion [13]. A numerical algorithm based on the augmented Lagrangian is also designed to solve the proposed convex optimization problem. As we described before, our strategy is to conduct local coordinates reconstruction for the KNN $N(i)$ of the $i$-th point. For convenience, we write $\{\bm{x}_{1},\bm{x}_{2},\ldots\bm{x}_{\ell}\}$ as $\ell$ points for $N(i)$ as the following model does not rely on a typical choice of $N(i)$.

2.1 Classical multidimensional scaling

The classical MDS is to find a configuration of a set of points $\{\bm{x}_{1},\bm{x}_{2},\ldots\bm{x}_{\ell}\}$ in $\mathbb{R}^{p}$ from the given distance matrix $D=(d_{ij}^{2})_{\ell\times\ell}$ 111For convenience, each entry of $D$ is specified by the square of distance. such that the Euclidean distance matrix $\left(\|\bm{x_{i}}-\bm{x}_{j}\|^{2}_{2}\right)_{ij}$ among these points is as close as possible to the given distance matrix $D$ [26, 10].

Definition 2.1**.**

*A distance matrix $D=(D_{ij})_{\ell\times\ell}$ is called Euclidean if there exist some points $\bm{x}_{1},\bm{x}_{2},\ldots,\bm{x}_{\ell}\in\mathbb{R}^{p}$ such that $D_{ij}=(\bm{x}_{i}-\bm{x}_{j})^{\top}(\bm{x}_{i}-\bm{x}_{j})$. Let’s write $X=(\bm{x}_{1},\bm{x}_{2},\ldots\bm{x}_{\ell})^{\top}$ as an $\ell\times p$ matrix and denote the centering matrix by $\displaystyle H=I_{l}-\frac{1}{l}\bm{1}\bm{1}^{\top}$. A Gram matrix $B=(B_{ij})$ is defined as $B=-\frac{1}{2}HDH$, equivalently, $D_{ij}=B_{ii}+B_{jj}-B_{ij}-B_{ji}$. *

Note that the Euclidean distance matrix is translation invariant, thus, the matrix $H$ essentially translates the point set centering at the origin. It is well-known that $D$ is Euclidean if and only if the Gram matrix $B$ is positive semidefinite. Furthermore, since $B=HXX^{\top}H$, it is clear that $\text{rank}(B)=\text{rank}(XX^{\top})=\text{rank}(X)=p$. In our case, we consider the number of points is much larger the dimension of points. This suggests the low-rank property of $B$ as its prior knowledge. If $D$ is Euclidean, and thus $B$ is positive semidefinite, then the eigen-decomposition of $B=U\Lambda U^{\top}$ provides the coordinate reconstruction $X=U_{p}\Lambda^{1/2}_{p}$, where $\Lambda_{p}$ is a diagonal matrix formed by the largest $p$ eigenvalues of $X$ and $U_{p}$ are the corresponding eigenfunctions.

2.2 Local coordinates reconstruction via positive-semidefinite matrix completion

Different from the MDS described as the above, one of the major challenges in our assumption is that only a small portion of $D$ is available. Generally speaking, it is impossible to recover the whole matrix only from its small portion of information. We here assume the information $D$ is randomly missing. Different models of semi-definite programming have been considered [1, 9] in this scenario. Here, we consider a different model inspired by the low-rank prior information of the Gram matrix and the recent advances from the low-rank matrix completion theory [13]. In other words, we are seeking for a symmetric positive-semidefinite low-rank matrix $B$ satisfying $B_{ii}+B_{jj}-B_{ij}-B_{ji}=D_{ij}$ for $(i,j)$ in the available portion. Mathematically, let’s write $\Omega\subset\{(i,j)~{}|~{}1\leq j<i<\ell\}$ as the available index set and write $\mathbb{S}=\{B=(b_{ij})\in\mathbb{R}^{\ell\times\ell}~{}|~{}B=B^{T}\}$. We further write the available part of $D$ as a vector $D_{\Omega}\in\mathbb{R}^{|\Omega|}$ and define a restriction operator $\mathcal{R}_{\Omega}:\mathbb{S}\rightarrow\mathbb{R}^{|\Omega|}$ by $\mathcal{R}_{\Omega}(B)_{ij}=B_{ii}+B_{jj}-B_{ij}-B_{ji}$ with $(i,j)\in\Omega$. Our problem becomes to find a positive semidefinite matrix $B\in\mathbb{S}$ satisfying $\mathcal{R}_{\Omega}(B)=D_{\Omega}$. It is clear that this problem is underdetermined and there are infinitely many solutions, most of which are far from the ground truth of $B$ and thus not applicable. Regularization method is a common way to handle this type of situation. Based on the prior knowledge of $B$ is a low rank matrix without knowing the exact rank $p$, we here consider the following optimization problem.

[TABLE]

Inspired by the matrix completion theory [13], we relax the above NP-hard problem by replacing $\mathrm{Rank}(B)$ as the nuclear norm $\|B\|_{*}$. In addition, the input distance matrix is invariant under translation of the reconstructed coordinates $\{\bm{x}_{i}\}$. This translation ambiguity can be fixed by requiring $\sum_{i}\bm{x}_{i}=\bm{0}$. This condition is equivalent to row sum of $B$ is zero, which can also be obtained from the definition of the Gram matrix $B$. Therefore, we propose the following matrix completion model:

[TABLE]

where $\bm{1}$ denotes a column vector with all elements constant $1$. Once the above Gram matrix is obtained, coordinates can be reconstructed by the eigen-decomposition of $B$ as the method used in the classical MDS described in section 2.1. Simultaneously, the eigen-decomposition of the reconstructed $B$ clearly shows the exact dimension $p$ by counting the cardinality of non-zero singular value. We would like to remark that the constraints and the restriction operator are different from the setting in [13] and it is highly nontrivial to check the restricted isometry property [40], therefore, theoretical analysis conducted in [13, 40] can not be directly applied in our problem. By considering a dual basis method, we can also theoretically show the exact completion can be achieved under certain coherence condition. More details about the theoretical analysis of the model will appear in our ongoing work [28]. Additionally, in a related work [59], the authors propose a semidefinite embedding (SDE) model as a kernel learning method for nonlinear dimension reduction, where no random missing but only KNN information of $D$ is considered, and the rank minimization is not considered but maximization of $Trace(B)$ is proposed for maximizing the variance in feature space.

2.3 Numerical Algorithm

For a convenience of designing a numerical algorithm of model (2), we first introduce some notations. Note that the set of symmetric matrix $\mathbb{S}$ can be viewed as $\mathbb{R}^{\ell(\ell+1)/2}$ due to a natural isomorphism $\zeta:\mathbb{S}\rightarrow\mathbb{R}^{\ell(\ell+1)/2},B\mapsto\hat{B}=(b_{11},\cdots,b_{\ell 1},b_{22},\cdots,b_{\ell 2},\cdots,b_{\ell\ell})^{T}$. We denote $\iota=\zeta^{-1}$ and further merge the two linear constraints by defining the following linear operator:

[TABLE]

and extend the vector $D_{\Omega}$ as $\tilde{D}_{\Omega}=(D_{\Omega},\vec{0})$. Note that $\|B\|_{*}$ is the same as $\text{Trace}(B)$ as $B$ is positive semi-definite. Therefore, the proposed model (2) can be written as

[TABLE]

The model (4) is a semi-definite programming (SDP) problem. It is well-known that the SDP problem can be solved by interior-points method [53, 56] as a canonical choice. However, when the number of constraint is with the order $\mathscr{O}(\ell^{2})$ as in (4), the computational time can be with the complexity $\mathscr{O}(\ell^{6})$ [60]. Thus, we use the alternating direction method [60] to save both the time and memory consumption. We first introduce an auxiliary variable $\hat{C}=\hat{B}$ and write (4) as:

[TABLE]

Then we introduce an augmented Lagrangian of model (5) as follows:

[TABLE]

where $H_{1}$ and $H_{2}$ are two dual variables. The variable splitting enables us to alternatively optimize $\hat{B}$ and $\hat{C}$. As (5) is convex, the saddle point of $\mathcal{L}_{\mu_{1},\mu_{2}}(\hat{B},\hat{C};H_{1},H_{2})$ is the solution which can be obtained by the following iterative procedure:

[TABLE]

Note that the model (5) is convex, therefore convergence of the above method can be guaranteed [54, 60].

The first subproblem is convex and differentiable. From its first order optimality condition, $\hat{B}^{k+1}$ solves the following linear system:

[TABLE]

where $\mathcal{A}^{*}$ is the adjoint operator of $\mathcal{A}$ defined as follows. It is clear that for any given $\tilde{U}=(U_{\Omega},\vec{q})\in\mathbb{R}^{\Omega}\times\mathbb{R}^{\ell}$ with $U_{\Omega}=\{u_{ij}\}_{(i,j)\in\Omega}$, the conjugate of $\mathcal{A}$ is given by $\mathcal{A}^{*}\tilde{U}=\zeta(U)+\zeta(Q)$. Here we define symmetric matrices $U$ and $Q$ by:

[TABLE]

and $U_{ij}=0$ for any $(i,j)\notin\Omega$. As $\mu_{1}\mathcal{A}^{*}\mathcal{A}+\mu_{2}I_{\ell}$ is a symmetric and positive definite matrix, there are numerous numerical solvers to approach the solution of the above problem. Here, we typically chose the conjugate gradient method to solve (8).

The second subproblem of (7) is to find the nearest point of $\hat{B}^{k+1}+H_{2}^{k}$ in the positive semi-definite cone. Its closed-form solution can be simply realized by truncating all the information in the negative eigenvalue part. More precisely, for a symmetric matrix $X$ with eigen-decomposition $X=V\Lambda V^{\top}$, we define an eigenvalue hard thresholding (EVHT) operator as follows:

[TABLE]

With this definition, the second step of (7) can be calculated as follows:

[TABLE]

We summarize the above iterative method in Algorithm 1. After finding $\hat{B}$ by Algorithm 1, eigen-decomposition of $\iota(\hat{B})$ can ultimately reconstruct the coordinates of points $\{\bm{x}_{1},\bm{x}_{2},\ldots\bm{x}_{\ell}\}$ from the incomplete distance $D_{\Omega}$.

Remark 2.2**.**

*In Algorithm 1, Step 2 is usually time consuming as the full eigen-decomposition has time complexity $\mathscr{O}(\ell^{3})$. Because ultimately the solution of $B$ and $C$ should be low rank, it is safe to only carry out the largest $m$ eigenfunctions by setting $m$ is relatively big enough comparing with the rank of $B$ but also small in terms of the matrix size. For instance, in our later numerical experiments for 3D point clouds, we use the MATLAB function ”eigs” to find the largest $20$ eigenvalues and truncate all the other parts. In this way, the time consumption of Step 2 can be apparently decreased, and numerical results suggest no difference between this simplification and calculating full eigen-decomposition. *

Remark 2.3**.**

In addition, we would like to remark that the distance matrix is also a low rank matrix provided by $\text{Rank}(D)\leq\text{Rank}(B)+2$ due to the relationship between $B$ and $D$ satisfying $D=-2B+\text{diag}(B)\bm{1}^{\top}+\bm{1}\text{diag}(B)^{\top}$, where $\text{diag}:\mathbb{R}^{\ell\times\ell}\rightarrow\mathbb{R}^{\ell\times 1}$ is the linear operator to project diagonal component of $B$ to a column vector. Therefore, we can also directly reconstruct $D$ based on its available information. This leads to a low rank minimization model as follows:

[TABLE]

*The above model is exactly the same as the general matrix completion theory discussed in [13], which also indicates that the minimum desired information of $D$ for successful reconstruction is proportion to the matrix rank. Note that the rank of $D$ is higher than the rank of matrix $B$. Plus, the less flexibility of $B$ due to the positive semi-definite condition can also help to reduce the requirement of valid distance information for successful reconstruction. Therefore, we would expect that model (12) requires a higher rate of valid distance than the proposed model (2). In fact, our numerical experiments in Section 5 also verify this observation. *

3 Solving PDEs on Manifolds from Incomplete Distance

In this section, after a brief review of the moving least square (MLS) method [33, 34] and the local mesh (LM) method [29] for solving PDEs on manifolds represented as point clouds, we demonstrate our proposed methods of solving an elliptic eigenvalue problem and a nonlinear hyperbolic equation on manifolds represented as incomplete inter-point distance information.

3.1 Approximating differential operators from reconstructed local coordinates

For traditional implicit or triangle mesh representation of surfaces in $\mathbb{R}^{3}$, implicit methods, parameterization methods, finite difference methods and finite element methods [38, 39, 8, 52, 7, 36, 49, 48, 57, 35, 11, 62, 41, 27, 20] have been proposed to solve PDEs on surfaces. However, it is not straightforward to use these methods in our case as the incomplete distance data structure has no global mesh structure and the ground manifold could be with dimension higher than two and in a high co-dimensional ambient space. As the definitions of differential operators on manifolds are coordinate invariant, we can approximate differential operators point-wisely based on its local coordinate reconstruction obtained from the matrix completion method proposed in Section 2. Namely, each given point, its reconstructed local neighborhood can be viewed as a point cloud, therefore, we apply moving least square (MSL) method [33, 34] or local mesh (LM) method [29] to approximate the differential operators. These methods can achieve high order accuracy and enjoy more flexibility. They can be applied to manifolds with arbitrary dimensions and codimensions. To make the paper self-sufficient, we briefly discuss these two methods in this section.

Denote the index set of the KNN of the $i$-th point as $N(i)$ and write $\mathcal{X}(i)=\{\bm{x}_{k}\in\mathbb{R}^{p}~{}|~{}k\in N(i)\}$ as the reconstructed Euclidean coordinates of the neighborhood $N(i)$. Based on $\mathcal{X}(i)$, its tangent space and normal space can be determined by standard principle component analysis (PCA) [23], which is provided by eigensystem of the covariance matrix $P_{i}$ of $N(i)$ defined on:

[TABLE]

where $c_{i}=\frac{1}{|N(i)|}\sum\limits_{k\in N(i)}\bm{x}_{k}$ is the centroid of $\mathcal{X}(i)$. If the intrinsic dimension of the manifold is $d$, then the jump of $P_{i}$’s eigenvalues guides the splitting $\mathbb{R}^{p}=\mathcal{T}_{i}\oplus\mathcal{N}_{i}$. Here $\mathcal{T}_{i}$ represents the tangent space spanned by $\{e_{i}^{1},\cdots,e^{d}_{i}\}$ corresponding to the $d$ largest eigenvectors of $P_{i}$, and $\mathcal{N}_{i}$ represents the normal space spanned by the rest of the eigenvectors of $P_{i}$. As $\{e_{i}^{1},\cdots,e_{i}^{p}\}$ forms an orthonormal basis near the point $\bm{x}_{i}$, a new virtual coordinates of $\bm{x}_{k}$ can be obtained by $(\bm{u}_{k},\bm{v}_{k})=(\{\langle\bm{x}_{k},e_{i}^{\alpha}\rangle\}_{\alpha=1}^{d},\{\langle\bm{x}_{k},e_{i}^{\beta}\rangle\}_{\beta={d+1}}^{p})$. Therefore, the manifold structure near $\bm{x}_{i}$ can be approximated by a degree 2 polynomial map $\displaystyle\bm{Z}_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p-d},\bm{Z}_{i}(\bm{u})=(\sum_{|\bm{\xi}|\leq 2}\bm{c}^{1}_{\bm{\xi}}\bm{u}^{\bm{\xi}},\cdots,\sum_{|\bm{\xi}|\leq 2}\bm{c}^{p-d}_{\bm{\xi}}\bm{u}^{\bm{\xi}})$ whose coefficients are given by the following moving least square problem:

[TABLE]

where the standard multi-index notation is used for $\bm{u}^{\bm{\xi}}=u_{1}^{\xi_{1}}\cdots u_{d}^{\xi_{d}}$ and $w(\cdot)$ is a weight function whose typical choice can be $w(d)=\exp(-\frac{d^{2}}{h^{2}})$ and $\displaystyle h=\max_{k\in N(i)}\|\bm{x}_{k}-\bm{x}_{i}\|$. Similarly, any function $f$ or vector filed $\bm{V}$ defined on $\mathcal{X}(i)$ can also be interpolated as a degree 2 polynomial $f(\bm{u})$ or $\bm{V}(\bm{u})$. Overall, we obtain a local polynomial interpolation of $\mathcal{X}(i)$ as $\bm{Z}_{i}(\bm{u})$ and a local polynomial interpolation of function defined on $\mathcal{X}(i)$ as $f(\bm{u})$. Therefore, the metric tensor $G=(g_{st})$ near $\bm{x}_{i}$ can be obtained by $g_{st}=\delta_{st}+\langle\frac{\partial{\bm{Z}_{i}}}{\partial u_{s}},\frac{\partial{\bm{Z}_{i}}}{\partial u_{t}}\rangle$. This enables us to approximate intrinsic differential operators at $\bm{x}_{i}$ such as:

[TABLE]

where $g=det(G)$, $(g^{st})=G^{-1}$ and the Einstein summation is used. Other differential operators can also be approximated using the similar method. After that, a finite difference type of method can be applied to solve differential equations on the given data. We refer [33, 34] for more details about this approach.

The local mesh method proposed in [29] is another way to approximate the above differential operators from the virtual local coordinates. Namely, at $i$-th point, a local connectivity can be constructed for the projection image of $\mathcal{X}(i)$ on $\mathcal{T}_{i}$ through the standard Delaunay triangulation. Thus, this connectivity can be directly inherited on $\mathcal{X}(i)$. If we write the simplex of the first ring of $\bm{x}_{i}$ by $\mathcal{R}(i)$, the differential operators at each point can be approximated by weighted average as follows:

[TABLE]

Moreover, a finite element type of method can be applied to estimate the mass matrix and stiffness matrix. Therefore, differential equations like Laplace-Beltrami eigenvalues problems can be solved. We refer [29] for more detailed discussion about construction of mass matrix and stiffness matrix.

We remark that the above constructions are conducted point-wisely. As long as a differential operator is well-defined on a manifold, namely, its definition does not depend on the choice of local coordinate, then the above procedure can consistently approximate the desired differential operators. A different strategy that can also be considered is to reconstruct coordinates first using the proposed Gram matrix completion algorithm, then approximate differential operators based on the global reconstructed coordinates. We would like to point out that the proposed methods of approximating differential operators based on local reconstruction enjoys advantages of computation efficiency and memory consumption. In fact, the most time-consuming part of the coordinate reconstruction Algorithm 1 is the eigenvalue hard thresholding. Consider a incomplete distance data with $n$ nodes, the complexity of each eigenvalue thresholding step is $\mathscr{O}(n^{2}m)$ for global reconstruction if only the largest $m$ eigenvalues are computed. In our local reconstruction strategy, we only reconstruct coordinates of $\ell$ nearest points near $\bm{x}_{i}$ with complexity $\mathscr{O}(n\ell^{2}m)$. When $n$ is very large (such as $16002$ used in our experiments) and $\ell$ is small ( $\leq$ 30 in our experiments ), the local reconstruction strategy is apparently much more efficient than global reconstruction strategy. Our numerical experiments discussed in Section 5 also support this point (see Table 2). Moreover, the local reconstruction method supports parallel computation which can further reduce time consumption. Additionally, the memory consumption of global reconstruction is $\mathscr{O}(n^{2})$ and the matrix is not sparse, which often crashes the program due to insufficient memory for large size data. However, our local reconstruction strategy only requires the memory $\mathscr{O}(\ell^{2})$ which avoids the problem of exceeding memory limit.

3.2 Solving differential equations based on incomplete distance

The first PDE we consider is an elliptic eigenvalue problem of the Laplace-Beltrami (LB) operator $\Delta_{\mathcal{M}}$. The LB operator is self-adjoint and elliptic, so its spectrum is discrete. We denote the eigenvalues of $-\Delta_{\mathcal{M}}$ as $0=\lambda_{0}<\lambda_{1}<\lambda_{2}<\cdots$ and the corresponding eigenfunctions as $\phi_{0},\phi_{1},\phi_{2},\cdots$ satisfying the following equations [15]:

[TABLE]

The set of LB eigenfunctions $\{\phi_{k}\}_{k=0}^{\infty}$ forms an orthonormal basis of the space of $L^{2}$ functions on $\mathcal{M}$. The set $\{\lambda_{k},\phi_{k}\}_{k=0}^{\infty}$ is called LB eigensystem of $\mathcal{M}$. Due to the intrinsic definition of the LB operator $\Delta_{\mathcal{M}}$, the induced LB eigensystem $\{\lambda_{k},\phi_{k}\}_{k=0}^{\infty}$ is also completely intrinsic to the manifold geometry and provides an intrinsic and systematic characterization of the underlying manifold geometry [5]. Recently, there have been increasing interests in using the LB eigen-geometry for 3D shape analysis as well as point clouds analysis [41, 32, 55, 46, 24, 51, 30, 12, 33, 31]. Therefore, it would be also important to consider solve the LB eigensystem for manifolds represented as incomplete inter-point distance, then all existing methods of conducting data analysis using the LB eigensystem can also be adapted to the incomplete distance data structure. Our numerical solvers are based on the numerical solvers for on manifolds represented as point clouds. First, we apply the matrix completion model for point-wise local coordinate reconstruction. Either MLS or LM can be then applied to obtain the discretized matrix form of the equation, where the only step we need to conduct is to use the reconstructed local coordinates for the LB operator discretization as we discussed in the section 3.1. As the LB operator is invariant under different choices of the local coordinates, the rigid motion ambiguity from the local coordinate reconstruction will not introduce inconsistency across different coordinate patches. This guarantees that our method can provide satisfactory numerical results as we illustrate in the numerical section. We would like to remark that it is relatively straightforward to have the local consistency of the differential operator as long as the available distance information satisfying certain incoherence condition as we discussed in our theoretical validation work [28]. However, the stability of the discretization is still open and will be explored in our future work.

The second equation we consider is the Eikonal equation, a special type of nonlinear hyperbolic PDE on manifolds. This equation is used to characterize the geodesic distance, an intrinsic measurement between two points, on a manifold. The Eikonal equation for the distance map $d$ to a given set $\Gamma$ on $\mathcal{M}$ can be stated as follows:

[TABLE]

when $\Gamma$ only includes a single point $p$, the distance map can be denoted as $d_{p}$. In practice, provided points coordinates and mesh structure, model (16) can be solved by fast marching and fast sweeping methods [45, 25, 63]. For manifolds represented as point clouds, a local mesh method based on fast marching is proposed in [29]. As long as the local coordinate reconstruction is obtained using the proposed matrix completion model, then a local Delaunay mesh structure can be constructed. After that, we repeat the local mesh method to conduct the fast marching as this approach only depends on the distance information of the first ring structure of the local mesh, which is again rigid motion invariant and will not be affected due to the rigid motion ambiguity from the local coordinate reconstruction.

4 Applications on Manifold Reconstruction and Dimension Reduction

As we mentioned before, solutions of differential equations on manifolds can provide global information for understanding data structure. Thus, PDEs can be viewed as “bridges” linking between local information and global information. This shares the same sprit as “think globally, fit locally” discussed in [43] although no PDEs are involved there. Therefore, without global reconstruction, some global analysis of point clouds such as pattern recognition, comparisons and classification can be further considered for data represented as incomplete distance. In this section, we illustrate results on applications of manifold reconstruction and dimension reduction based on solutions of the Laplace-Beltrami eigenvalue problem on incomplete distance data.

4.1 Manifold reconstruction from distance via patch stitching using LB eigenfunctions

The patch stitching problem is to reconstruct the coordinate of a point cloud $P\in\mathbb{R}^{n\times d}$ from coordinates of $L$ subsets (referred as patches) $\{\Omega_{j}\}_{j=1}^{L}$, where $\Omega_{j}$ denotes the index set for patch $j$. In practice, for each patch $\Omega_{j}$, the local coordinates $Q_{j}\in\mathbb{R}^{n_{j}\times d}$, which are provided as input information or are reconstructed from distance information, have the same pair-wise distance and geometry as the restricted global coordinate $P_{\Omega_{j}}$ with possible rigid motion difference. In other words, $P_{\Omega_{j}}=Q_{j}R_{j}+\bm{1}_{j}b_{j}$, where $R_{j}\in O(d)=\{R\in\mathbb{R}^{d\times d}~{}|~{}R^{\top}R=I_{d}\}$ is an orthogonal matrix for rotation and reflection, $\bm{1}_{j}\in\mathbb{R}^{n_{j}\times 1}$ is a vector with all element equal to 1, and $b_{j}\in\mathbb{R}^{1\times d}$ is the translation vector. A straightforward approach to solve $\{R_{j}\}$ and $\{b_{j}\}$ is to minimize the quadratic loss:

[TABLE]

Such model is a non-convex problem since $\{R_{j}\}$ are required to optimize over the non-convex domain of orthogonal transforms. In [18], the authors proposed a three-stage patch stitching method from the local coordinate by synchronizing the reflection, rotation and translation successively. However, the quality of synchronization is strongly relies on the proportion of the overlapping index. More recently, a spectral relaxation method is proposed to solve the above problem in [14]. This method essentially relax the set of orthogonality constraints and perform well if neighborhood patches has reasonable enough overlapping points.

Here, we propose an alternative method to reconstruct the global coordinate from local coordinate, or merely incomplete distance information of each patch. Based on solution of LB eigen-problem incomplete distance information discussed in Section 3.2, we consider the LB eigen-system as a ”bridge” to connect the information of different patches. This idea of this new method is essentially align each local patch in the spectral domain instead of the original Euclidean coordinate domain. As the global information has been captured from the LB eigen-system, the overlapping of different patches is not directly required here. This approach is inspired by approximating coordinates using LB eigenfunctions. For a given point cloud $P\in\mathbb{R}^{n\times d}$, it can be approximated by the linear combination of the LB eigenfunctions as $P\approx\Phi\bm{\alpha}$. Here, $\Phi=(\phi_{1},\cdots,\phi_{N})\in\mathbb{R}^{n\times N}$ represents the first $N$ LB eigenfunctions $\{\phi_{i}\}_{i=1}^{N}$ and $\bm{\alpha}=(\alpha_{1},\cdots,\alpha_{d})\in\mathbb{R}^{N\times d}$ represents coefficients of coordinates as column vectors. Since the mass matrix is generally unknown without the prior knowledge of the structure information, the coefficients matrix $\bm{\alpha}$ can be obtained in the following least square sense:

[TABLE]

Consider the relationship between restricted global coordinate $\{P_{\Omega_{j}}\}$ and the input local coordinate $\{Q_{j}\}$, we already have that for any $1\leq j\leq L$, $P_{\Omega_{j}}=Q_{j}R_{j}+\bm{1}b_{j}$. Therefore the norm $\|P-\Phi\bm{\alpha}\|_{2}^{2}$ in local patch $Q_{j}$ can be equivalently written as $\|P_{\Omega_{j}}-\Phi_{\Omega_{j}}\bm{\alpha}\|_{2}^{2}=\|Q_{j}R_{j}+\bm{1}_{j}b_{j}-\Phi_{\Omega_{j}}\bm{\alpha}\|_{2}^{2}=\|Q_{j}+\bm{1}_{j}b_{j}R_{j}^{\top}-\Phi_{\Omega_{j}}\bm{\alpha}R_{j}^{\top}\|_{2}^{2}$. Since $b_{j}R_{j}^{\top}$ is still a translation vector in $\mathbb{R}^{d}$, the norm can be simplified as $\|Q_{j}-\Phi_{\Omega_{j}}\bm{\alpha}R_{j}-\bm{1}_{j}b_{j}\|_{2}^{2}$ without loss of generality. As a result, with the given $\Phi$ and all the local coordinate $\{Q_{j}\}$, we propose the following model for stitching point clouds $\{Q_{j}\}_{j=1}^{L}$:

[TABLE]

Once $\{R_{j}\}$ and $\{b_{j}\}$ are obtained, we can find all the adjusted coordinates $P_{\Omega_{j}}=Q_{j}R_{j}+\bm{1}_{j}b_{j}$, which leads to the global coordinate $P$ ultimately.

Different from the synchronization method and the spectral relaxation method considered in [18, 14], we design the following method to solve the nonconvex problem (19) by updating $\bm{\alpha}$, $\{R_{j}\}$ and $\{b_{j}\}$ iteratively:

[TABLE]

where we update $\bm{\alpha}$ and $\{b_{j}\}$ using Newton’s method and it is straightforward to check

[TABLE]

We next describe a method of updating the orthogonality constrained variable $\{R_{j}\}$. The operator $\mathcal{K}$ proposed in [61] is designed for searching the gradient descent direction on the Stiefel manifold, the geometric description of the set of orthogonal matrices. To realize the operator $\mathcal{K}(\nabla R_{j}^{k})$, we first define the skew-symmetric operator as

[TABLE]

then the new trial point $R_{j}^{k+1}$ satisfing $R_{j}^{k+1}(R_{j}^{k+1})^{\top}=I$ can be generated by

[TABLE]

where the step size $\delta$ can be obtained by setting fixed value or line search methods. More detailed discussion about the orthogonality preserving property and convergence of this approach can be found in [61]. With the above explanation of (20), we summarize an algorithm of solving (19) as Algorithm 2.

Based on the above method, we therefore propose a global coordinate reconstruction algorithm from merely random missing distance in each local patches. Given the local patches $\{\Omega_{j}\},1\leq j\leq L$ with patch size not larger than $K$, we assume the local distance information include $r\%$ of the local distance of each patches. The full reconstrution procedure can be summarized as follows:

Step 1

Using the method discussed in Section 3.2 with local mesh reconstruction method to compute the first $N$ Laplace-Beltrami eigenfunctions $\Phi$.

Step 2

Using Algorithm 1 to generate the virtual local coordinate $\{Q_{j}\},1\leq j\leq L$ for each patch.

Step 3

Using Algorithm 2 to find the global coordinate $P$ from $\Phi$ and $\{Q_{j}\}$.

5 Numerical Results

In this section, numerical tests are presented to illustrate the proposed methods for solving PDEs on manifold represented as incomplete inter-point distance information and the application to manifold stitching. First, we demonstrate the effectiveness of the proposed matrix completion method for coordinate reconstruction. We also illustrate the phase transition curve of successful reconstruction. Second, we test our method of computing the LB eigen-problem and illustrate that the proposed method of solving PDEs through local coordinate reconstruction is much more efficient than solving PDEs from global coordinate reconstruction. In addition, we also test the problem of solving Eiknoal equations and illustrate the accuracy improvement of this approach to the canonical Dijkstra method. As applications, we also demonstrate preliminary results on point clouds reconstruction using global information from LB eigenfunctions and its potential application to dimension reduction problem. All numerical experiments are implemented by MATLAB in a PC with a 32G RAM and 16 dual-core 2.7 GHz CPUs.

5.1 Coordinates reconstruction from incomplete distances

The local coordinate reconstruction method based on matrix completion is a crucial part of solving differential equations on a manifold represented as incomplete inter-point distance data. To demonstrate the effectiveness of the method proposed in Section 2, we first test this method for global coordinate reconstruction, although our numerical method for solving differential equations only conduct coordinate reconstruction for a certain neighborhood of each point. Starting with a full distance matrix $D\in\mathbb{R}^{n\times n}$, the actual number of freedom in $D$ is $(n^{2}-n)/2$ as the diagonal elements of $D$ are [math] and $D$ is symmetric. We randomly chose $\gamma(n^{2}-n)/2,(\gamma\in[0,~{}1])$ number of entries of $D$ as the input incomplete distance matrix and apply Algorithm 1 to find a Gram matrix $B$. After that, global coordinates can be consequently generated from $B$ by its eigen-decomposition as we described in Section 2.1.

Our first numerical test is conducted for checking successful reconstruction rate of reconstructing the input incomplete distance matrices. We test our model for different size matrices with fixed rank. By uniformly choosing $n$ points on the unit sphere, we obtained the corresponding pairwise Euclidean distance matrix $D_{T}$ and the associated Gram matrix $B_{T}$ as ground truth. After that, we randomly choose $\gamma n(n-1)/2$ entries of $D_{T}$ as available information and apply our reconstruction algorithm. We mean that the reconstruction is successful if the relative error between $B_{T}$ and the reconstructed Gram matrix under the Frobenius norm is less than $10^{-3}$. Given $n$ and $\gamma$, we run 50 tests and record the rate of successful reconstruction as $\displaystyle\rho=\frac{\#\mbox{ successful reconstruction}}{50}$. Based on the same input, Figure 1 plots the phase transition of $\rho$ using Gram matrix completion model (2) (the left image in Figure 1) and distance completion model (12) (the right image in Figure 1) in terms of the sampling rate $\gamma$ and the logarithm of total number of points $\log(n)$. Both images clearly show that for a fixed rank problem, successful reconstruction needs less portion of distance information as number of points increases. As we expect, the Gram matrix completion model (2) has a larger domain of successful reconstruction than the model (12) using distance matrix completion. This is because the rank of Gram matrix is less than the rank of distance matrix and an additional semi-positive definite constraint is also impose for Gram matrix which further decrease the number of freedom. Therefore, we use model (2) as our major local coordinate reconstruction tool to approximate differential operators in the rest of numerical experiments.

We further test distance matrices from different data including point clouds sampled from the unit sphere(1002 points), a cow surface (2602 points), and a Swiss roll surface (2048 points) using our Gram matrix completion model. Table 1 reports the relative error of the Gram matrix and rates of successful reconstruction out of 50 tests for each $\gamma$. Figure 2 shows that with 1% distance information, the reconstructed coordinates cannot formulate an outline of the unit sphere while the same rate of distance information can successfully reconstruct the cow shape. This is compatible with our observation showed in Figure 1, which indicates less percentage information is need for successfully reconstructing matrices with larger size. In addition, Table 1 and Figure 2 also demonstrate both quantitatively and qualitatively that the proposed reconstruction model can provide very good coordinate reconstruction once the available distance information is sufficient.

5.2 Solve LB eigenvalue problem from distance data with missing values

Our second numerical experiment illustrates the proposed method of solving the LB eigenvalue problem on manifolds represented as incomplete inter-point distance information. Considering a point cloud $\{p_{1},\cdots,p_{n}\}$ uniformly sampled on a given manifold $\mathcal{M}\subset\mathbb{R}^{p}$, we compute the associated distance matrix $D$ provided by pairwise Euclidean distance for this point cloud. Our numerical experiments illustrate the proposed method for computing the LB eigenvalue problem only based on the incomplete distance. More precisely, we first uniformly randomly choose a subset $\Omega_{\gamma}\subset\{(i,j)~{}|~{}1\leq j<i\leq n\}$, such that $|\Omega_{\gamma}|=\left\lceil\gamma n(n-1)/2\right\rceil$. Then incompletion distance information is provided as

[TABLE]

where $N_{\ell}(i)$ denotes the index set of the nearest $\ell$ points to $p_{i}$. Thus, $D_{\gamma,\ell}$ roughly contains $\gamma\ell/n$ portion of $D$, which contains a very small portion of $D$ if $\ell\ll n$ as we showed in our experiments. In addition, we measure computation accuracy using the normalized error $\displaystyle E_{\max,k}=\max_{i}\frac{|\tilde{\lambda}_{k,i}-\lambda_{k}|}{|\lambda_{k}|}$, where $\tilde{\lambda}_{k,i}$ are numerical approximation from our method, and $i$ represents the multiplicity.

We first demonstrate our method on different data sampled from the unit sphere. It is known that the $k$-th LB eigenvalue of the unit sphere is given by $\lambda_{k}=k(k+1)$ with multiplicity $2k+1$. This allowss us to illustrate numerical accuracy of our methods on the unit sphere case. In the first test, we assume the $K$-nearest neighborhood (KNN) is provided and sample distance matrices for point clouds on the unit sphere by setting different values of $\gamma=100\%,80\%,50\%,30\%$ and with different KNN size $\ell=6,9,18,30$, respectively. Given a KNN size $\ell$, the settings of local sampling rate should above the border line as shown in Figure 1 to guarantee the exact local coordinate reconstruction. Conversely, given a fixed sampling rate of local distance, we need to set $\ell$ large enough to guarantee the successful reconstruction. With these input data, we apply the proposed strategy for approximating the LB operator only from local coordinates reconstruction and compute the first 100 LB eigenvalues. We typically chose approximation results of LB eigenvalues $\lambda=20$ and $\lambda=72$ to conduct accuracy analysis by illustrating the numerical error $E_{\max,n}$ in the right table of Figure 3. In addition, we list results from global coordinate reconstruction based on $3\%$ of distance information. As a reference, we also list computation results from classical finite element method from all the exact coordinates of points and the mesh information, which can be experimentally regarded as the best computation result based on the whole distance information. It is clear that our method provides satisfactory approximation of the LB eigenvalue problem on the unit sphere. Convergence of our method is illustrated by the error curves showed in the left two plots in Figure 3. Moreover, we also observe the second order convergence numerically as $E_{\max,k}$ is approximately reduced by half if the number of total point doubled. This is compatible with convergence behavior of the moving least square method for solving the LB eigenvalue problem on point clouds discussed in [33, 34].

We also report total time consumption in Table 2 for computing the first 100 LB eigenvalues with different setting of $\gamma$ and $\ell$ in this test. Note that the major time-consuming part is local coordinates reconstruction based on incomplete distance information, whose complexity is dependent on the local size parameter $\ell$ and the total number of points $n$. Theoretically, our local coordinate reconstruction method generally has the time complexity $\mathcal{O}(\min(\ell^{2}nm,\ell^{3}n))$, while the global coordinate reconstruction has time complexity $\mathcal{O}(n^{2}m)$. Here, $m$ is denoted as the maximum number of eigenvalues used in the eigenvalue hard thresholding step of algorithm 1, which is typically chosen as 20 in our experiments. For large data sets satisfying $\ell\ll n$, local reconstruction method can save a huge amount of time. This advantage is also illustrated in Table 2, where local reconstruction method for solving LB eigenvalue problem is about 150 times faster than global reconstruction method if we choose $n=16002$ and $\ell=6$.

In practice, the input condition might be even weaker than our previous experiment as the exact KNN information may not be directly obtained from incomplete local distances. Therefore, we also test our method for reconstructing local coordinates only relying on KNN information from incomplete distance information. In this scenario, we use the same setting as the previous test but require a larger rate of available distance information to have successfully local coordinate reconstruction. Our numerical results reported in Figure 4 also demonstrate nearly second order convergence of LB eigenvalue problem on the unit sphere. We further test the robustness of our method to distance with Gaussian perturbation. In this setting, we assume the input data is a incomplete distance $D_{80\%,30}$ corrupted by different levels of Gaussian noise with standard deviation $\sigma=2\%d_{\max},5\%d_{\max},10\%d_{\max},15\%d_{\max}$ respectively, where $d_{\max}=\max_{i,j}\{D(i,j)\}$. The left table in Figure 5 reports numerical accuracy for different levels of Gaussian noise. It is clear that our method still provides a reasonable good approximation of LB eigenvalues. Moreover, we also plot the corresponding LB eigenfunctions in Figure 6, where LB eigenfunctions are color-coded on the unit sphere by setting red for positive values and blue for negative values. This figure illustrates consistent distribution patterns of LB eigenfunctions for noisy distance case. In addition, we also test our method of computing LB eigenvalue problems for incomplete distance information from point clouds sampled from more complicated manifolds such as armadillo (16519 points) and kitten (2884 points) surfaces. Figure 7 reports serval LB eigenfunctions on armadillo and kitten surfaces based on incomplete distance matrix $D_{50\%,30}$.

As we discussed before, our methods can also be applied on manifolds with dimension more than 2 and co-dimension more than 1. Our third experiment tests the proposed method for computing LB eigenvalue problem on a 2 dimensional torus $T^{2}=\{(\cos\theta_{1},\sin\theta_{1},\cos\theta_{2},\sin\theta_{2})\in\mathbb{R}^{4}~{}|~{}\theta_{1},\theta_{2}\in[0,2\pi)\}$ in $\mathbb{R}^{4}$ and a 3 dimension torus $T^{3}=\{(\cos\theta_{1},\sin\theta_{1},\cos\theta_{2},\sin\theta_{2},\cos\theta_{3},\sin\theta_{3})\in\mathbb{R}^{6}~{}|~{}\theta_{1},\theta_{2},\theta_{3}\in[0,2\pi)\}$ in $\mathbb{R}^{6}$. Solutions of the corresponding LB eigenvalue problem for both objects have closed forms, which enable us to conduct computation accuracy comparisons. For $T^{2}$, we choose $D_{80\%,30}$ from the corresponding distance matrix obtained from different sample size on $T^{2}$. Similarly setting is considered for $T^{3}$ but using $D_{85\%,30}$ with slightly large portion of information as the rank for each local Gram matrix for $T^{3}$ is higher than the one for $T^{2}$. This requires more information of the distance matrix from the matrix completion theory. Nevertheless, we only use less than $1\%$ of the distance information for our experiments as $\ell=30,n\geq 2500$. We remark that our computation does not assume the exact KNN information is available. This is the reason that we need higher portion of local information for accurate local coordinate reconstruction. In addition, we do not assume the dimension information is available in our computation. As we demonstrated in the first experiment, the value of $\gamma$ can be even lower if the exact KNN information is provided. We compute the first 100 LB eigenvalues for $T^{2}$ and $T^{3}$ using incomplete distance matrices with different sample size . Figure 8 reports the convergence and relative errors of our method. Both cases demonstrate approximately second order convergence.

5.3 Solve Eikonal equation from incomplete distance

In this subsection, we test our method for solving a special hyperbolic equation, the Eikonal equation, on manifolds represented as incomplete distance information. In our experiments, assuming the KNN information is given as prior knowledge, we considered the incomplete local distance of uniform/non-uniform sampled point cloud of a unit sphere, in which the number of points vary from $1002$ to $16002$. In particular, we assume the point cloud always contains two points with coordinate $(0,0,1)$ and $(0,0,-1)$, namely, the north pole and the south pole (See Figure 9). To evaluate the quality of the reconstruction of distance map, we compare the relative error of the geodesic distance from the north pole to south pole (referred as $E_{se}$), and from the north pole to all other points (referred as $E_{av}$). The table in Figure 9 shows that our approach based on the fast marching method is much more accurate than the Dijkstra’s method [19]. Moreover, relative errors of our method decrease as the number of points increase. Furthermore, for the uniform sample cases, the relative error of geodesic distance from north pole to south pole has first order convergence with respect to the density of points using our method, while the Dijkstra’s method does not have the same convergence property.

Similarly, we also compute the distance function on a Swiss roll:

[TABLE]

based on incomplete local distance from uniform/non-uniform sampling, in which two diagonal points are fixed and their position can be seen in Figure 9. Using $D_{60\%,20}$ from different sample sizes vary from $1002$ to $16002$, the bottom Table in Figure 9 shows that our method can approximate the distance function with much less error than the Dijkstra’s method, and approximate first convergence can also be observed from the case of the uniform sampled Swiss roll.

5.4 Global coordinate reconstruction from using patch stitching

In this experiment, we test our proposed global reconstruction model as an application of using global information from solving the LB eigenproblem based on incomplete distance. We assume the input manifold $\mathcal{M}\subset\mathbb{R}^{3}$ sampled as a point cloud and has been separated as partially overlapped $L$ patches, namely, $\mathcal{M}=\bigcup_{j=1}^{L}\Omega_{j}$. We generate distance matrices $\{D_{j}\}_{j=1}^{L}$ from pairwise Euclidean distance on each patch and randomly choose $50\%$ information of each distance matrix. Based on these incomplete $\{D_{j}\}_{j=1}^{L}$, we first compute the first 100 LB eigenfunctions of $\mathcal{M}$ using the proposed method discussed in Section 3. After that, we reconstruct coordinates of each $\Omega_{j}$ from its incomplete distance using our matrix completion method. Finally, global reconstruction of $\mathcal{M}$ is obtained by the manifold stitching model discussed in Section 4, which can handle the coordinate inconsistency from each local reconstruction. In Figure 10, we report numerical results for reconstructing an armadillo surface (with 16519 points, 28 patches) and a kitten surface (with 2884 points, 30 patches). The two left images in Figure 10 illustrate coordinate reconstruction for each patch, the two middle images show the corresponding global reconstruction and the two right energy curves indicate the convergence of the proposed global reconstruction algorithm. It is clear to see that the proposed method can successfully reconstruct a manifold from a set of incomplete distance matrices from separated patches, while the global matrix completion model can not handle this case as the missing distance information is quite coherent. Moreover, compared to direct global coordinate reconstruction, the method of local distance reconstruction and patch stitching saves more computational time. To validate this claim, we measured that the $50\%$ of local distance $\{D_{j}\}_{j=1}^{L}$ actually contains $3.59\%$ of full distance $D$ for armadillo surface and $3.45\%$ for kitten surface. Then we also implement the direct global coordinate reconstruction for these two surfaces using the same rate of random missing global distance $D$. Table 3 shows that the computation time is much smaller using the stitching scheme especially for the armadillo surface with a large number of points.

We further conduct computation based on an incomplete geodesic distance matrix provided by a set of uniformly sampled points on the Swiss roll $S$ used in the previous example. It is straightforward to check that the geodesic distance between any two points $((t_{1}+0.1)\cos(t_{1}),(t_{1}+0.1)\sin(t_{1}),8\pi s_{1})$ and $((t_{2}+0.1)\cos(t_{2}),(t_{2}+0.1)\sin(t_{2}),8\pi s_{2})$ is given by $\sqrt{(8\pi(s_{1}-s_{2}))^{2}+(\mathcal{G}(t_{2})-\mathcal{G}(t_{1}))^{2}}$, where $\mathcal{G}(t)=\frac{1}{2}[(t+0.1)\sqrt{1+(t+0.1)^{2}}+\log(|\sqrt{1+(t+0.1)^{2}}+(t+0.1)|)]$. In fact, this Swiss roll $S$ is isometrical to a 2D flat domain $[\mathcal{G}(0),\mathcal{G}(4\pi)]\times[0,8\pi]$. We uniformly sample 2048 points on $S$ and construct a distance matrix $D$ using pairwise geodesic distance. We test two ways of reconstruct the Swiss roll based on incomplete information of $D$. Firstly, similar as patch stitching simulations, we assume the Swiss roll is separated as $23$ patches with partial overlap. Then we chose only $50\%$ of local geodesic distance $\{D_{j}\}_{j=1}^{23}$, which is essentially $3.61\%$ coherent sampling of totally geodesic distance $D$. The bottom row of figure 11 shows that the global stitched coordinate is a 2D rectangle, which is identical to the theoretical ground truth of the dimension reduced Swiss roll. Secondly, we assume $3.61\%$ random sampled global geodesic distance is available and apply our coordinate reconstruction algorithm to have the global coordinate reconstruction for the Swiss roll. Top right image in figure 11 shows that direct global coordinate reconstruction from incoherent random missing distance is also applicable. Table 3 shows that with same sampling rate from the whole distance matrix, the direct global coordinate reconstruction takes much more computational time. The reason can be explained similarly as computation complexity discussion in subsections 3.1 and  5.2. The step of solving PDEs from local distance has time complexity $\mathcal{O}(\min(l^{2}mn,l^{3}n))$ while direct global coordinate reconstruction has time complexity $\mathcal{O}(mn^{2})$ which is larger if total number of points $n$ is large enough. Although the eigen-decomposition of Laplace-Beltrami matrix also involves global computation but it needs to be performed only once while eigenvalue thresholding is required in every iteration for direct global coordinate reconstruction. Therefore, it is not uprising that the stitching scheme saves a lot of time. We remark that our method of using geodesic distance to reconstructing manifolds is highly related to an important topic called nonlinear dimension reduction in machine learning and statistics community [42]. In our future work, we will investigate more along this direction based on our methods.

6 Conclusion

In this paper, we proposed a framework for discretizing PDEs on manifolds represented as incomplete inter-point distance information. Our approach is to conduct PDE discretization point-wisely from a local coordinate reconstruction. This strategy successfully avoids a time-consuming step of global coordinate reconstruction and leads to a method with complexity linearly scaling to the number of sample size. Our local reconstruction model is inspired from the recent advances of low-rank matrix completion theory, which guarantee accuracy local coordinate reconstruction by only requiring a very small portion of distance information. Our method can be viewed as natural extensions of the moving least square method and the local method method for solving PDEs on point clouds to manifold-structured data represented as incomplete inter-point distance information. As an application of using solutions of PDEs, we propose a new manifold reconstruction model by stitching local patches on the spectrum domain. Intense numerical experiments indicate effectiveness and efficiency of our methods for solving LB eigenvalue problem, Eikonal equations and manifold reconstruction using stitching methods based on global information from the LB eigensystem.

7 Acknowledgement

We would like to thank Prof. Xiangxiaong Zhang’s suggestions about optimizing distance matrix directly as we discussed in Remark 2.3.