A Unified Definition and Computation of Laplacian Spectral Distances

TL;DR

This paper introduces a unified framework for defining and computing Laplacian spectral distances in both continuous and discrete settings, enabling efficient approximation and analysis of spectral kernels on graphs and manifolds.

Contribution

It extends recent discrete results to the continuous case, providing a novel definition and efficient computation methods for Laplacian spectral kernels and distances.

Findings

Discrete spectral kernels are approximated by sparse linear systems.

The proposed methods are efficient and well-conditioned.

The paper discusses the stability and optimality of spectral distances.

Abstract

Laplacian spectral kernels and distances (e.g., biharmonic, heat diffusion, wave kernel distances) are easily defined through a filtering of the Laplacian eigenpairs. They play a central role in several applications, such as dimensionality reduction with spectral embeddings, diffusion geometry, image smoothing, geometric characterisations and embeddings of graphs. Extending the results recently derived in the discrete setting~\citep{PATANE-STAR2016,PATANE-CGF2017} to the continuous case, we propose a novel definition of the Laplacian spectral kernels and distances, whose approximation requires the solution of a set of inhomogeneous Laplace equations. Their discrete counterparts are equivalent to a set of sparse, symmetric, and well-conditioned linear systems, which are efficiently solved with iterative methods. Finally, we discuss the optimality of the Laplacian spectrum for the approximation of the spectral kernels, the relation between the spectral and Green kernels, and the stability of the spectral distances with respect to the evaluation of the Laplacian spectrum and to multiple Laplacian eigenvalues.

Equations

$$\left\{\begin{array}[]{ll}f=\sum_{n=0}^{+\infty}\alpha_{n}\phi_{n},&\alpha_{n}:=\langle f,\phi_{n}\rangle_{2};\\ \|f\|_{2}^{2}=\sum_{n=0}^{+\infty}\alpha_{n}^{2},&\|\nabla f\|_{2}^{2}=\sum_{n=0}^{+\infty}\alpha_{n}^{2}\lambda_{n}.\end{array}\right.$$

$$\Vert r_n{\Vert }_2^2\le \frac{\Vert \nabla f{\Vert }_2^2}{{\lambda }_{n+1}}.$$

$${\Vert f-\sum _{i=1}^n⟨f,{\psi }_i{⟩}_2{\psi }_i\Vert }_2\le \alpha \frac{\Vert \nabla f{\Vert }_2}{{\lambda }_{n+1}},\forall f.$$

$$\begin{array}{r}{\Delta }^{i}f=\sum _{n=0}^{+\infty }{\lambda }_n^i⟨f,{\varphi }_n{⟩}_2{\varphi }_n,({\Delta }^{†}{)}^{i}f=\sum _{n=1}^{+\infty }\frac{1}{{\lambda }_n^i}⟨f,{\varphi }_n{⟩}_2{\varphi }_n.\end{array}$$

$${\Phi }_{\rho }f=\sum _{n=0}^{+\infty }{\alpha }_n{\Delta }^{n}f=\sum _{n=0}^{+\infty }\rho ({\lambda }_n)⟨f,{\varphi }_n{⟩}_2{\varphi }_n,$$

$$K_{\rho }(\mathbf{p},\mathbf{q})=\sum _{n=0}^{+\infty }\rho ({\lambda }_n){\varphi }_n(\mathbf{p}){\varphi }_n(\mathbf{q})$$

$$\begin{array}[]{l}\Phi_{\rho}\Phi_{1/\rho}\Phi_{\rho}=\Phi_{\rho},\quad\Phi_{1/\rho}\Phi_{\rho}\Phi_{1/\rho}=\Phi_{1/\rho},\\ \langle\Phi_{1/\rho}\Phi_{\rho}f,g\rangle_{2}=\langle f,\Phi_{1/\rho}\Phi_{\rho}g\rangle_{2}.\end{array}$$

$$\left\{\begin{array}[]{ll}\langle f,g\rangle:=\langle\Phi_{1/\rho}f,\Phi_{1/\rho}g\rangle_{2}=\sum_{n=0}^{+\infty}\frac{\langle f,\phi_{n}\rangle_{2}\langle g,\phi_{n}\rangle_{2}}{\rho^{2}(\lambda_{n})},&(a)\\ d^{2}(f,g):=\|f-g\|^{2}=\sum_{n=0}^{+\infty}\frac{|\langle f-g,\phi_{n}\rangle_{2}|^{2}}{\rho^{2}(\lambda_{n})}.&(b)\end{array}\right.$$

$$\begin{array}{rl}{d}^2(\mathbf{p},\mathbf{q})& :=\Vert {\delta }_{\mathbf{p}}-{\delta }_{\mathbf{q}}{\Vert }^2\\ & {=}_{\text{Eq. }(\text{ref}eq:FILT-SCAL-PRODb)}\sum _{n=0}^{+\infty }\frac{|{\varphi }_n(\mathbf{p})-{\varphi }_n(\mathbf{q}){|}^2}{{\rho }^2({\lambda }_n)}\\ & {=}_{\text{Eq. }(\text{ref}eq:FILT-SCAL-PRODa)}\Vert {\Phi }_{1/\rho }({\delta }_{\mathbf{p}})-{\Phi }_{1/\rho }({\delta }_{\mathbf{q}}){\Vert }_2^2\\ & =\Vert K_{1/\rho }(\mathbf{p},\cdot )-K_{1/\rho }(\mathbf{q},\cdot ){\Vert }_2^2\\ & =K_{1/\rho }(\mathbf{p},\mathbf{p})-2K_{1/\rho }(\mathbf{p},\mathbf{q})+K_{1/\rho }(\mathbf{q},\mathbf{q});\end{array}$$

$$\begin{array}{rl}{d}^2(\mathbf{p},\mathbf{q})& =\Vert K_{1/\rho }(\mathbf{p},\cdot )-K_{1/\rho }(\mathbf{q},\cdot ){\Vert }_2^2\\ & =\sum _{n=0}^{+\infty }{\rho }^{-2}({\lambda }_n){\varphi }_n(\mathbf{p}){\varphi }_n(\mathbf{p})-2\sum _{n=0}^{+\infty }{\rho }^{-2}({\lambda }_n){\varphi }_n(\mathbf{p}){\varphi }_n(\mathbf{q})\\ & +\sum _{n=0}^{+\infty }{\rho }^{-2}({\lambda }_n){\varphi }_n(\mathbf{q}){\varphi }_n(\mathbf{q}).\end{array}$$

$$\left({\Phi }_{\rho }+\delta \mathcal{E}\right)\varphi (\delta )=\mu (\delta )\varphi (\delta ),\varphi (0)=\varphi ,\mu (0)=\mu .$$

$$\mathcal{E}\varphi +{\Phi }_{\rho }{\varphi }^{\prime }(0)={\mu }^{\prime }(0)\varphi +\mu {\varphi }^{\prime }(0).$$

$$|{\mu }^{\prime }(0)|=|⟨\varphi ,\mathcal{E}\varphi {⟩}_2|\le \Vert \mathcal{E}{\Vert }_2\Vert \varphi {\Vert }_2^2=\Vert \mathcal{E}{\Vert }_2;$$

$$(s-\mu {)}^m=\frac{p(s)}{q(s)}\approx \frac{O(\delta )}{q(s)},\delta \to 0,\text{as }p(s)\to 0,\text{ when }s\to \mu ;$$

$$d_k(\mathbf{p},\mathbf{q}):=\sum _{n=0}^k\frac{|{\varphi }_n(\mathbf{p})-{\varphi }_n(\mathbf{q}){|}^2}{{\rho }^2({\lambda }_n)}$$

$$K_{\Delta }(\mathbf{p},\mathbf{q})=\sum _{n=1}^{+\infty }\frac{1}{{\lambda }_n}{\varphi }_n(\mathbf{p}){\varphi }_n(\mathbf{q});$$

$$K_{{\Delta }^2}(\mathbf{p},\mathbf{q}):=\sum _{n=1}^{+\infty }\frac{1}{{\lambda }_n^2}{\varphi }_n(\mathbf{p}){\varphi }_n(\mathbf{q});$$

$$K_t(\mathbf{p},\mathbf{q})=\sum _{n=0}^{+\infty }\exp (-t{\lambda }_n){\varphi }_n(\mathbf{p}){\varphi }_n(\mathbf{q});$$

$$d(f,g)=\Vert u{\Vert }_2,u={\Phi }_{1/\rho }(f-g)⟺{\Phi }_{\rho }u=f-g,$$

$$\mathcal{K}(\mathcal{N}):=\{K_{\rho }:\mathcal{N}\times \mathcal{N}\to \mathbb{R},K_{\rho }\text{ spectral kernel in Eq. (refeq:SPECTRAL-KERNEL)}\},$$

$$\begin{gathered} \Vert K_{\rho }-K_{\phi }{\Vert }_2^2=\sum _{n=0}^{+\infty }|\rho ({\lambda }_n)-\phi ({\lambda }_n){|}^2 \\ \le \Vert \rho -\phi {\Vert }_2^2, \end{gathered}$$

$$c_{rl}(s):=\frac{p_l(s)}{q_r(s)}=\frac{{\beta }_0+{\beta }_{1}s+\dots +{\beta }_l{s}^l}{{\alpha }_0+{\alpha }_{1}s+\dots +{\alpha }_r{s}^r},s\in [a,b].$$

$$\sum _{i=0}^l{a}_i{s}^i+\sum _{i=1}^r{b}_i{s}^{-i}=c_{rl}(s).$$

$$\Delta g_1=\Delta {\Delta }^{†}f=f-⟨f,{\varphi }_0{⟩}_2{\varphi }_0,$$

$$g_i:=({\Delta }^{†}{)}^{i}f={\Delta }^{†}({\Delta }^{†}{)}^{i-1}f={\Delta }^{†}{g}_{i-1},i\ge 2,$$

$$({\Phi }_{1/\rho }^{(r)}f{)}_{r=0}^{+\infty },{\Phi }_{1/\rho }^{(r)}f:=\sum _{n=0}^{+\infty }{c}_{rl}({\lambda }_n)⟨f,{\varphi }_n{⟩}_2{\varphi }_n,$$

$$\begin{array}{rl}{\Vert {\Phi }_{1/\rho }^{(r)}f-{\Phi }_{1/\rho }f\Vert }_2^2& \le \Vert c_{rl}-1/\rho {\Vert }_{\infty }^2\sum _{n=0}^{+\infty }|⟨f,{\varphi }_n{⟩}_2{|}^2\\ & ={\sigma }_{rl}^2\Vert f{\Vert }_2^2,{\sigma }_{rl}\approx \mathcal{O}(s^{r+l+1}),s\to 0;\end{array}$$

$${\rho }^{-1}({\lambda }_i){=}_{\text{Eq. }(\text{ref}eq:FUNCT-OPER)}⟨{\Phi }_{1/\rho }{\varphi }_i,{\varphi }_j{⟩}_2={\mathbf{x}}_i^{\top }{\mathbf{K}}_{1/\rho }^{\top }\mathbf{B}{\mathbf{x}}_j,\forall i=1,\dots ,n.$$

$$d^2({\mathbf{p}}_i,{\mathbf{p}}_j)=\Vert {\mathbf{K}}_{1/\rho }({\mathbf{e}}_i-{\mathbf{e}}_j){\Vert }_{\mathbf{B}}^2=\sum _{l=1}^n\frac{|⟨{\mathbf{x}}_l,{\mathbf{e}}_i-{\mathbf{e}}_j{⟩}_{\mathbf{B}}{|}^2}{{\rho }^2({\lambda }_l)},$$

$$d({\mathbf{p}}_i,{\mathbf{p}}_j)=\Vert {\mathbf{K}}_{1/\phi }({\mathbf{e}}_i-{\mathbf{e}}_j){\Vert }_{\mathbf{B}}.$$

Full text

1 Introduction

Spectral kernels and distances are easily defined through a filtering of the Laplacian eigenpairs and include random walks [41], biharmonic [28, 45], heat diffusion [10, 11, 17, 25], wave [12] kernels and distances. Laplacian spectral distances have been applied to shape segmentation [19] and comparison [11, 17, 35, 50] with multi-scale and isometry-invariant signatures [33, 32, 31]. Main applications of the heat kernel include dimensionality reduction [7, 53, 16, 4], diffusion geometry [6, 47], graphs’ embeddings [30] and analysis [18], shape comparison [9, 13, 43], data representation [54] and classification [34, 46, 49].

Overview and contribution

Extending the results recently derived in the discrete setting [38, 39] to the continuous case, we propose a novel definition of the Laplacian spectral kernels and distances, whose approximation requires the solution of a set of inhomogeneous Laplace equations. Their discrete counterparts are equivalent to a set of sparse linear systems, which are efficiently solved with iterative methods. To this end, the Laplacian spectral distances $d^{2}(\mathbf{p}_{i},\mathbf{p}_{j}):=\sum_{n=0}^{+\infty}|\phi_{n}(\mathbf{p}_{i})-\phi_{n}(\mathbf{p}_{j})|^{2}/\rho^{2}(\lambda_{n})$ are defined by filtering the Laplacian spectrum $(\lambda_{n},\phi_{n})_{n=0}^{+\infty}$, where $\rho:\mathbb{R}^{+}\rightarrow\mathbb{R}$ is a proper filter function.

The proposed spectral kernels and distances are independent of the data dimensionality (e.g., surfaces, volumes, $n$-dimensional data), of the discretisation of the input domain, and of the Laplace-Beltrami operator. The Laplacian spectral distances generalise the bi-harmonic [28] ($\rho(s):=s^{2}$), wave [12] ($\rho_{t}(s):=\exp(ist)$), and diffusion [14, 21, 36] ($\rho_{t}(s):=\exp(st)$) distances, Mexican hat wavelets [23] ($\rho(s):=s^{1/2}\exp(s^{2})$), and random walks [41]. The behaviour of the filter $\rho$, and in particular the convergence of $1/\rho$ to zero, determines the main properties of the resulting spectral distances, such as smoothness, encoding of local and global shape properties (e.g., Gaussian curvature, geodesic distance), localisation in time and frequency [21]. Through a proper filter design and a learning process [2, 8], we can obtain spectral kernels and distances that are also invariant with respect to isometric transformations or discriminative of a given shape class.

Recalling the definition of the spectral operator in [39, 37, 38], we review equivalent representations of the spectral kernels and distances in terms of the spectral norm of the $\delta$-functions, the Laplacian spectrum, the spectral operator, and the spectral kernel (Sect. 2). Through these relations, the Laplacian spectral kernels and distances can be computed by approximating the filter with a rational function and converting their evaluation to the solution of a set of differential equations that involve only the Laplace-Beltrami operator.

The proposed discrete approximation of the spectral kernels (Sect. 3) is equivalent to solving $r$ sparse, symmetric, and well-conditioned linear systems, where $r$ is the degree of the rational polynomial approximation of the filter. Applying an iterative linear solver, the computation of the spectral kernel takes $\mathcal{O}(r\tau(n))$ time, where $\tau(n)$ typically varies from $\tau(n)=n$ to $\tau(n)=n\log n$ and depends on the number $n$ of shape samples and on the sparsity of the coefficient matrix. In a similar way, the computation of the spectral distance between two points has the same order of complexity of the evaluation of the corresponding spectral kernel at the same input points.

Through the spectrum-free approximation, the computed spectral kernel is not affected by the Gibbs phenomenon, as a consequence of the high accuracy of the rational approximation with respect to the low-pass filter associated with the truncated spectral approximation. Furthermore, the spectrum-free approximation is free of user-defined parameters (e.g., the number of Laplacian eigenpairs).

Finally (Sects. 4, 5), we experimentally verify the accuracy, efficiency, and numerical robustness of the definition and spectrum-free computation of the spectral kernels and distances.

Novelties with respect to previous work

With respect to previous work and our recent results on the definition and computation of discrete spectral distances [39, 38], the main novelties of this paper are

the discussion of the optimality of the Laplacian spectrum for the approximation of the spectral kernel (Sect. 2.1);

• 2.

the study of the relation between the spectral and Green kernels (Sect. 2.4), associated with differential operators (Sect. 2.2);

• 3.

the analysis of the stability of the spectral distances with respect to the evaluation of the Laplacian spectrum and to multiple Laplacian eigenvalues (Sect. 2.3). In fact, the evaluation of the Laplacian spectrum is generally affected by small perturbations of the Laplace-Beltrami operator;

• 4.

a novel spectrum-free computation of the spectral distances (Sect. 2.5), which is achieved by properly approximating the associated filter with a rational polynomial and is expressed in terms of the canonical basis;

• 5.

the computation of the discrete spectral distances (Sect. 3.2), which reduces to the solution of a set of sparse, symmetric, and well-conditioned linear systems. This analogy between the continuous and discrete cases shows the generality of the proposed approach and it has not been addressed by previous work. Furthermore, these results are complementary and more general than the spectrum-free approaches presented in [39, 38];

• 6.

a new upper bound to the conditioning number of the linear systems associated with the spectrum-free approximation, which confirms the numerical stability and well-conditioned computation of the class of the spectrum-free approaches (Sect. 3.1);

• 7.

new experiments (Sect. 4), which enrich the tests initially presented in [39, 38] and address the robustness of the computation of the Laplacian spectral distances and kernels with respect to the data resolution, partial sampling, geometric or topological noise, and deformations.

2 Laplacian spectral distances

Firstly, we present the Laplace-Beltrami operator, its spectrum, and the optimality of the Laplacian eigenbasis (Sect. 2.1). Recalling the definition of the spectral operator and kernel introduced in [39, 37, 38], we review equivalent representations of the spectral distances in terms of the spectral norm of the $\delta$-functions, the Laplacian spectrum, the spectral operator, and the spectral kernel (Sect. 2.2). As novel contribution, we discuss the stability of the spectral distances with respect to the evaluation of the Laplacian spectrum and to multiple Laplacian eigenvalues, which are generally affected by small perturbations of the Laplace-Beltrami operator (Sect. 2.3). Finally, we introduce a novel relation between the spectral and Green kernels (Sect. 2.4) and propose a novel spectrum-free computation of the spectral distances in the continuous case, which is presented in a way analogous to the discrete case [38, 39] and shows the generality of the proposed approach (Sect. 2.5).

2.1 Laplace-Beltrami operator

Let $\mathcal{N}$ be a smooth surface, possibly with boundary, equipped with a Riemannian metrics and let us consider the inner product $\langle f,g\rangle_{2}:=\int_{\mathcal{N}}fg\textrm{d}\mu$ defined on the space $\mathcal{L}^{2}(\mathcal{N})$ of square integrable functions on $\mathcal{N}$ and the corresponding norm $\|\cdot\|_{2}$. The Laplace-Beltrami operator $\Delta$ is self-adjoint $\langle\Delta f,g\rangle_{2}=\langle f,\Delta g\rangle_{2}$, and positive semi-definite $\langle\Delta f,f\rangle_{2}\geq 0$, $\forall f,g$, [44].

Laplacian eigenbasis and its optimality

Recalling that the Laplace-Beltrami operator is self-adjoint and positive semi-definite, it has an orthonormal eigensystem $(\lambda_{n},\phi_{n})_{n=0}^{+\infty}$, $\Delta\phi_{n}=\lambda_{n}\phi_{n}$, $\lambda_{0}=0$, $\lambda_{n}\leq\lambda_{n+1}$, in $\mathcal{L}^{2}(\mathcal{N})$. Expressing a function $f$ in $\mathcal{L}^{2}(\mathcal{N})$ in terms of the Laplacian eigenfunctions, we have that

[TABLE]

According to [1], the Laplacian eigenfunctions are an optimal basis for the representation of signals with bounded gradient magnitude. In fact, the spectral decomposition $f_{n}=\sum_{i=0}^{n}\langle f,\phi_{i}\rangle_{2}\phi_{i}$ is optimal for the approximation of functions with $\mathcal{L}^{2}$ bounded gradient magnitude; i.e., the residual error $r_{n}:=f-f_{n}$ is bounded as

[TABLE]

The spectral decomposition is also optimal for the approximation of functions with respect to the error estimation in (1). In fact, for any $0\leq\alpha<1$ there is no integer $n$ and no sequence $(\psi_{i})_{i=0}^{n}$ of linearly independent functions in $\mathcal{L}^{2}$ such that

[TABLE]

Indeed, the Laplacian eigenfunctions are an optimal basis for the definition of scalar functions on a given domain, such as the Laplacian spectral kernels and distances on 3D shapes (Sect. 3). From the computational point of view, they have two main drawbacks: a high computational cost and storage overhead, which prevent the evaluation of a large number of Laplacian eigenpairs, and numerical instabilities with respect to the surface discretisation [38, 39].

2.2 Theoretical background: spectral kernels and distances

We briefly review the equivalent definitions of the spectral kernels (Sect. 2.2.1) and distances (Sect. 2.2.2), which will be useful to introduce their spectrum-free computation and main properties.

2.2.1 Spectral operator and kernel

According to [38, 39], let $\rho:\mathbb{R}^{+}\rightarrow\mathbb{R}$ be a positive filter map that is square integrable and admits the power series $\rho(s)=\sum_{n=0}^{+\infty}\alpha_{n}s^{n}$. Considering the orthonormal Laplacian eigensystem $(\lambda_{n},\phi_{n})_{n=0}^{+\infty}$, the spectral representation of the functions $\Delta^{i}f$ and $(\Delta^{{\dagger}})^{i}f$ is

[TABLE]

If $\lambda_{k}=0$, then we neglect the corresponding entry in $(\Delta^{{\dagger}})^{i}f$. In $\mathcal{L}^{2}(\mathcal{N})$, we define the spectral operator

[TABLE]

which is linear, continuous, self-adjoint, and $\Phi_{\rho}f=\langle K_{\rho},f\rangle_{2}$, with

[TABLE]

spectral kernel. The spectral kernel operator and kernel are well-posed if the filter is bounded or square-integrable. The spectral kernel is also symmetric and self-adjoint, as a consequence of its definition. From these relations, it follows that

the spectral kernel $K_{\rho}(\mathbf{p},\cdot)=\Phi_{\rho}\delta_{\mathbf{p}}$ is achieved as the action of the spectral operator on the $\delta$-function at $\mathbf{p}$;

• 2.

the spectral kernel is the solution to the differential equation $\Phi_{1/\rho}K_{\rho}(\mathbf{p},\cdot)=\delta_{\mathbf{p}}$. In Sect. 2.4, this property is further investigated in terms of the relation between the spectral and Green kernels.

To express the spectral distances in terms of the spectral operator (Sect. 2.5), we show that the pseudo-inverse of $\Phi_{\rho}$ is induced by the filter function $1/\rho$; i.e., $\Phi_{\rho}^{{\dagger}}=\Phi_{1/\rho}$. In fact, from the spectral representation (3) of $\Phi_{1/\rho}$ and $\Phi_{\rho}$, we get that

[TABLE]

In the following, we assume that both $\Phi_{\rho}$ and $\Phi_{1/\rho}$ are well-defined; for instance, this hypothesis is satisfied if $\rho$ is not null only on a compact interval of $\mathbb{R}^{+}$ and it is valid in the discrete case (Sect. 3), where the Laplacian spectrum belongs to the interval $\mathcal{I}:=[0,\lambda_{\max}(\tilde{\mathbf{L}})]$, with $\lambda_{\max}(\tilde{\mathbf{L}})$ maximum Laplacian eigenvalue.

2.2.2 Spectral distances

Through the spectral operator, we introduce the scalar product and the corresponding distance as [38, 39]

[TABLE]

Indicating with $\delta_{\mathbf{p}}$ the map that takes value $1$ at $\mathbf{p}$ and [math] otherwise, and selecting $f:=\delta_{\mathbf{p}}$, $g:=\delta_{\mathbf{q}}$ in Eq. (6b), the spectral distance (Fig. 1) on $\mathcal{N}$ is defined as

[TABLE]

i.e., these equivalent formulations involve the Laplacian spectrum, the spectral operator and kernel. The third equality follows from the identity $\Phi_{\rho}(\delta_{\mathbf{p}})=K_{\rho}(\mathbf{p},\cdot)$ and will be applied to the computation of the spectral distances (Sect. 3.2); in fact, it is independent of the evaluation of the Laplacian spectrum. The last identity is achieved by applying the relation

[TABLE]

2.3 Stability of the eigenpairs of the spectral operator

Generalising the results in [38, 39], we show that the computation of a single eigenvalue of the spectral operator is numerically stable and instabilities are generally due to repeated or close eigenvalues. Firstly, we notice that if $\lambda$ is a Laplacian eigenvalue of multiplicity $m$ then $\mu:=\rho(\lambda)$ is an eigenvalue of $\Phi_{\rho}$ and its multiplicity is equal to or greater than $m$. The corresponding eigenfunction is the Laplacian eigenfunction associated with the eigenvalue $\lambda$; i.e., $\Delta\phi=\lambda\phi$ and $\Phi_{\rho}\phi=\rho(\lambda)\phi$.

We perturb the spectral operator by $\delta\mathcal{E}$, $\delta\rightarrow 0$, and compute the eigenpair $(\mu(\delta),\phi(\delta))$ of the corresponding operator $\Phi_{\rho}+\delta\mathcal{E}$; i.e.,

[TABLE]

Recalling that the derivative of $\mu(\delta)$ measures the variation of the eigenvalue, deriving (7) with respect to $\delta$, and evaluating the resulting relation at [math], we have that

[TABLE]

From (8) and the self-adjointness of $\Phi_{\rho}$, it follows that $\langle\phi,\Phi_{\rho}\phi^{\prime}(0)\rangle_{2}=\mu\langle\phi,\phi^{\prime}(0)\rangle_{2}$. Multiplying both sides of (8) by $\phi$ and applying the previous identity we get that

[TABLE]

i.e., the computation of the eigenvalue of $\Phi_{\rho}$ with multiplicity one is stable.

Assuming that $\mu:=\rho(\lambda)$ is an eigenvalue of $\Phi_{\rho}$ with multiplicity $m$ and rewriting the characteristic polynomial as $p(s)=(s-\mu)^{m}q(s)$, where $q(\cdot)$ is a polynomial of degree $n-m$ and $q(\mu)\neq 0$, we get that

[TABLE]

i.e., $s\approx\mu+O(\delta^{\frac{1}{m}})$. Indeed, modifying the Laplacian matrix in such a way that the filtered eigenvalues are perturbed by $\delta:=10^{-m}$ corresponds to a change of order $0.1$ in $\mu$ (i.e., $s\approx\mu+0.1$) and this amplification becomes larger as the multiplicity of the eigenvalue increases.

While multiple eigenvalues are typically associated with symmetric shapes, numerically close or switched eigenvalues are present regardless of the surface regularity. This unstable computation of multiple eigenpairs of the spectral operator generally affects the accuracy of the truncated spectral approximation

[TABLE]

of the corresponding distances. In fact, each filtered eigenvalue $\rho(\lambda_{n})$, $n=1,\ldots,k$, which appears at the denominator of $d_{k}$, can further accentuate the numerical error of $\lambda_{n}$ in the distance computation. Furthermore, the computation of the Laplacian spectrum is time-consuming and it is difficult to properly select the number of eigenpairs that is necessary to accurately approximate the spectral distance. Indeed, we propose an evaluation of the spectral distance that is independent of the computation of the Laplacian spectrum and is equivalent to a set of differential equations involving only the Laplace-Beltrami operator and its pseudo-inverse. This novel spectrum-free computation is a generalisation of the discrete approach recently presented in [38, 39].

2.4 Relation between spectral and Green kernels

To study the relation between the Green and the spectral kernels, let us consider a linear differential operator $\mathcal{L}$ and the corresponding Green kernel $K:\mathcal{N}\times\mathcal{N}\rightarrow\mathbb{R}$ such that $\mathcal{L}K(\mathbf{p},\mathbf{q})=\delta(\mathbf{p}-\mathbf{q})$. Then, the solution to the differential equation $\mathcal{L}u=f$ is expressed in terms of the Green kernel as $u(\mathbf{p})=\langle K(\mathbf{p},\cdot),f\rangle_{2}$.

Noting that the eigensystem $(\mu_{n},\psi_{n})_{n=0}^{+\infty}$ of the integral operator $\mathcal{A}_{K}f:=\langle K(\cdot,\cdot),f\rangle_{2}$ induced by the Green kernel satisfies the relation $\mathcal{L}\psi_{n}=\mu_{n}^{-1}\psi_{n}$, $\mu_{n}\neq 0$, we have that the eigensystem of $\mathcal{L}$ is $(\mu_{n}^{-1},\psi_{n})_{n=1}^{+\infty}$. Indeed, $\mathcal{L}$ and $\mathcal{A}_{K}$ have the same eigenfunctions and reciprocal eigenvalues. In particular, the spectral representation of the Green kernel $K(\mathbf{p},\mathbf{q})=\sum_{n=0}^{+\infty}\mu_{n}\psi_{n}(\mathbf{p})\psi_{n}(\mathbf{q})$ is uniquely defined by the spectrum of $\mathcal{L}$. Combining the previous results with the properties of the spectral operator and kernel (Sect. 2.2.2), we get that

for the harmonic operator $\mathcal{L}=\Delta$, the Green kernel is the commute-time kernel

[TABLE]

• 2.

for the bi-Laplacian operator $\mathcal{L}:=\Delta^{2}$, the Green kernel is the bi-harmonic kernel

[TABLE]

• 3.

for the diffusion operator $\mathcal{L}=\exp(t\Delta)$, the Green kernel is the diffusion kernel (Fig. 2)

[TABLE]

• 4.

for the spectral operator $\mathcal{L}=\Phi_{1/\rho}$, the Green kernel is the spectral kernel $K_{\rho}$ defined in Eq. (4).

2.5 Spectrum-free approximation of kernels and distances

Recalling the relation $\Phi_{\rho}^{{\dagger}}=\Phi_{1/\rho}$ in Eq. (5) and noting that

[TABLE]

the spectral distance is equivalent to the norm of (i) the function $u=\Phi_{1/\rho}(f-g)$ and (ii) the solution of the differential equation $\Phi_{\rho}u=f-g$. Through these relations, the spectral kernels will be computed by approximating the filter with a rational function (Sect. 2.5.1) and by converting the evaluation of the corresponding distances (Sect. 2.5.2) to the solution of a set of differential equations that involve the Laplace-Beltrami operator. The choice of one of these two equivalent representations depends on the selected filter and the simplicity of evaluating either $\Phi_{\rho}$ or $\Phi_{1/\rho}$.

2.5.1 Kernel approximation, convergence, and accuracy

Let us introduce the space of the Laplacian spectral kernels as

[TABLE]

and equipped with the $\mathcal{L}^{2}(\mathcal{N}\times\mathcal{N})$ scalar product. Since

[TABLE]

the approximation of a given spectral kernel $K_{\rho}$ with a new kernel $K_{\varphi}$ in $\mathcal{K}(\mathcal{N})$ is reduced to the approximation of $\rho$ by $\varphi$ on a proper subspace of functions (e.g., the space of polynomials or rational polynomials). The class of functions used for the approximation of the input filter $\rho$ is selected in such a way that $K_{\varphi}$ provides a good approximation of the input filter $K_{\rho}$ and is easily computable.

2.5.2 Rational approximation based on the canonical basis

We apply the Padé-Chebyshev rational approximation to the filter map $1/\rho$. Let $\mathcal{R}_{r}^{l}$ be the space of all rational functions

[TABLE]

We briefly recall that the representation of the rational approximation is not unique, unless we impose that it is irreducible; for example, by choosing $q_{r}(a)=1$. Given a filter $\rho:\,[a,b]\rightarrow\mathbb{R}$, there exists a unique best approximation of $1/\rho$ in $\mathcal{R}_{r}^{l}$ with respect to the $\ell_{\infty}$ norm [20] (Ch. 9), which is represented in terms of the canonic basis $\mathcal{B}:=\{s^{i}\}_{i=0}^{l}\cup\{s^{-i}\}_{i=1}^{r}$.

Recalling that the dimension of the space $\mathcal{R}_{r}^{l}$ of rational polynomial of degree $(r,l)$ is $l+r+1$, let us express the best rational approximation $c_{rl}$ of $\rho$ in terms of these basis functions. To this end, we can either apply algebraic rules or impose interpolating constraints at $(l+r+1)$ points in order to compute the new coefficients through the identity

[TABLE]

The canonical basis used to represent the rational approximation of the input filter simplifies the computation of the spectral kernel; in fact, we need to evaluate only the functions $\Delta^{i}f$ and $(\Delta^{{\dagger}})^{i}f$. Instead of applying the spectral representation (2), we define a recursive procedure as follows. To compute $g_{1}=\Delta^{{\dagger}}f$, let us multiply both sides by $\Delta$ and notice that

[TABLE]

where $\phi_{0}=\mathbb{1}$ is the constant eigenfunction equal to $1$. By definition of pseudo-inverse, $g_{1}$ is the least-squares solution to the equation $\Delta g_{1}=\tilde{f}$, where $\tilde{f}$ is equal to $f$ minus its mean $\langle f,\phi_{0}\rangle_{2}\phi_{0}$. For the general case, we apply the recursive relation

[TABLE]

which reduces to the previous case. In a similar way, $g_{i}=\Delta^{i}f$ is calculated as $h_{i}=\Delta h_{i-1}$, $i=1,\ldots,r$, with $h_{0}:=f$. Indeed, the proposed approach requires the solution of $(l+r)$ Laplace equation with a different right-hand side.

Convergence and accuracy

To verify that the sequence

[TABLE]

induced by the rational polynomial approximation $c_{rl}$ of $1/\rho$, converges to $\Phi_{1/\rho}f$, we apply the upper bound

[TABLE]

where $\sigma_{rl}$ is the approximation error between $1/\rho$ and $c_{rl}$. For the diffusion operator [52], $\sigma_{rr}=10^{-5}$.

3 Spectrum-free approximation of kernels and distances

We briefly summarise the discretisation of the spectral kernel and distances [38, 39], which is then used to establish a simple relation between the spectral kernel and its pseudo-inverse (Sect. 3.1). Then, we introduce the spectrum-free approximation of the spectral kernels and distances (Sect. 3.2) and discuss the computational cost of the main steps of the proposed approach (Sect. 3.3).

3.1 Discretisation of spectral kernels and distances

Let us consider a (triangular, polygonal, volumetric) mesh $\mathcal{M}:=(\mathcal{P},\mathcal{T})$, which discretises a domain $\mathcal{N}$, where $\mathcal{P}:=\{\mathbf{p}_{i}\}_{i=1}^{n}$ is the set of $n$ vertices and $\mathcal{T}$ is the connectivity graph. On $\mathcal{M}$, a scalar function $f:\mathcal{M}\rightarrow\mathbb{R}$ is identified with the vector $\mathbf{f}:=(f(\mathbf{p}_{i}))_{i=1}^{n}$ of $f$-values at $\mathcal{P}$.

Let us introduce the Laplacian matrix $\tilde{\mathbf{L}}:=\mathbf{B}^{-1}\mathbf{L}$, where $\mathbf{L}$ and $\mathbf{B}$ are the stiffness and mass matrix respectively (e.g., cotangent [40], Voronoi-cotg [15], linear FEM [42] weights), and the corresponding spectral decomposition is $\mathbf{L}\mathbf{X}=\mathbf{B}\mathbf{X}\Lambda$, $\mathbf{X}^{\top}\mathbf{B}\mathbf{X}=\mathbf{I}$, where $\mathbf{X}:=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]$ is the eigenvectors’ matrix and $\Lambda$ is the diagonal matrix of the eigenvalues $(\lambda_{i})_{i=1}^{n}$. Analogous discretisations apply to polygonal [3, 22] and tetrahedral [27, 51] meshes, or point sets [29].

Under these assumptions, we introduce the spectral kernels and distances (Sect. 3.1.1), together with their spectrum-free approximation (Sect. 3.1.2), which is based on a rational approximation of the input filter.

3.1.1 Discrete spectral kernels and distances

The spectral operator $\Phi_{1/\rho}$ is discretised by the spectral kernel matrix $\mathbf{K}_{1/\rho}$ such that

[TABLE]

Indeed, $\mathbf{K}_{1/\rho}=\mathbf{X}\rho^{{\dagger}}(\Lambda)\mathbf{X}^{\top}\mathbf{B}$ is the pseudo-inverse of the spectral kernel $\mathbf{K}_{\rho}=\rho(\tilde{\mathbf{L}})$, which is a filtered version of the Laplacian matrix. Here, $\rho^{{\dagger}}(\Lambda)=\textrm{diag}(1/\rho(\lambda_{i}))_{i=1}^{n}$ and its entry is null if $\rho(\lambda_{i})=0$. Then, the discrete spectral distances [38, 39] are

[TABLE]

where $\mathbf{e}_{i}$ is the vector of the canonical basis of $\mathbb{R}^{n}$, $\langle\mathbf{f},\mathbf{g}\rangle_{\mathbf{B}}:=\mathbf{f}^{\top}\mathbf{B}\mathbf{g}$ and $\|\mathbf{f}\|_{\mathbf{B}}^{2}:=\mathbf{f}^{\top}\mathbf{B}\mathbf{f}$ are the scalar product and the norm induced by the mass matrix, respectively. In this case, we have derived the spectral distances by applying the continuous expression of the spectral operator instead of the spectral kernel, as done in [38, 39].

Analogously to the definitions in Sect. 2.2, the spectral distance is equal to the norm $d(\mathbf{f},\mathbf{g}):=\|\mathbf{u}\|_{\mathbf{B}}$ of the solution to the linear system $\mathbf{K}_{\rho}\mathbf{u}=\mathbf{f}-\mathbf{g}$ or to the norm of the vector $\mathbf{u}=\mathbf{K}_{1/\rho}(\mathbf{f}-\mathbf{g})$. In a similar way, the spectral distance between two points reduces to

[TABLE]

Indeed, this approximation of the spectral distances involves the spectral kernel only and allows us to bypass numerical inaccuracies due to repeated or close Laplacian eigenvalues [20] ($\S$ $7$).

3.1.2 Spectral kernel approximation and conditioning

Since the eigenvalues of the spectral kernel matrix $\mathbf{K}_{\rho}$ are $(\rho(\lambda_{i}))_{i=1}^{n}$, the approximation of $\mathbf{K}_{\rho}$ with a new kernel $\mathbf{K}_{\varphi}$ reduces to the approximation of the corresponding filters with respect to the $\ell_{\infty}$ norm; in fact,

[TABLE]

The approximation $\varphi$ of $\rho$ is computed on the interval $[0,\lambda_{\max}(\tilde{\mathbf{L}})]$, where the maximum Laplacian eigenvalue is evaluated by the Arnoldi method [20], or is set equal to the upper bound $\lambda_{\max}(\tilde{\mathbf{L}})\leq\min\{\max_{i}\{\sum_{j}\tilde{L}(i,j)\},\max_{j}\{\sum_{i}\tilde{L}(i,j)\}\}$ [26, 48].

Conditioning of the spectral kernel

We now analyse the conditioning of the spectral kernel. Assuming that $\rho$ is an increasing function (i.e., $1/\rho$ is a low pass filter), the conditioning number of the filtered Laplacian matrix is bounded as

[TABLE]

and it is ill-conditioned when $\rho(0)$ is close to zero or $\rho$ is unbounded. If $\rho$ is bounded and $\rho(0)$ is not too close to [math], then the filtered Laplacian matrix is well-conditioned. If $\rho(0)$ is null, then we consider the smallest and not null filtered Laplacian eigenvalue at the denominator of the previous relation.

3.2 Spectrum-free approximation of distances: canonical basis

Analogously to the discussion in Sect. 2.5.2, the basis $\mathcal{B}$ used to represent the rational approximation of the input filter simplifies the computation of the spectral kernel; in fact, we need to compute only the following vectors $\tilde{\mathbf{L}}^{i}\mathbf{f}$ and $\tilde{(\mathbf{L}}^{i})^{{\dagger}}\mathbf{f}$. Applying the spectral representation of the Laplacian matrix

[TABLE]

we can represent its powers and the corresponding pseudo-inverse matrices as

[TABLE]

To compute $\mathbf{g}_{1}=\tilde{\mathbf{L}}^{{\dagger}}\mathbf{f}$, let us multiply both sides by $\tilde{\mathbf{L}}$ and notice that

[TABLE]

i.e., $\mathbf{g}_{1}$ is the least-squares solution to the sparse and symmetric linear system

[TABLE]

where $\mathbf{g}_{0}$ is achieved by subtracting to $\mathbf{f}$ its mean value $\langle\mathbf{f},\mathbf{1}\rangle_{\mathbf{B}}$ (c.f., Eq. (10)). For the general case, we apply the recursive relation

[TABLE]

which reduces to the previous case (Fig. 3). In a similar way, $\mathbf{g}_{i}=\tilde{\mathbf{L}}^{i}\mathbf{f}$ is calculated by recursively solving the sparse, symmetric, and positive-definite linear systems $\mathbf{B}\mathbf{g}_{i+1}=\mathbf{L}\mathbf{g}_{i}$, $i=1,\ldots,r-1$, with $\mathbf{g}_{0}:=\mathbf{f}$.

3.3 Computational cost

The truncated spectral approximation, whose computational cost depends on the sparsity degree of the Laplacian matrix, takes from $\mathcal{O}(kn\log n)$ to $\mathcal{O}(kn^{2})$ time, where $k$ is the number of selected eigenpairs. Selecting a rational approximation of the input filter of degree $(r,l)$, the evaluation of the corresponding spectral kernel is reduced to solve $l$ linear systems whose coefficient matrix is $\mathbf{B}$ and $r$ linear systems whose coefficient matrix is $\mathbf{L}$. Through iterative solvers, the computational cost is $\mathcal{O}((r+l)\tau(n))$, where $\tau(n)$ is the cost for the solution of a sparse linear system, which varies from $\mathcal{O}(n)$ to $\mathcal{O}(n^{2})$, according to the sparsity of the coefficient matrix, and it is $\mathcal{O}(n\log n)$ in the average case.

In a similar way, the computation of the spectral distance between two points has the same order of complexity of the evaluation of the corresponding spectral kernel at the same input points. For the computation of the one-to-all spectral distance, we lump or pre-factorise $\mathbf{B}$ (if not already diagonal) and pre-factorise $\mathbf{L}$. Then, the overall computational cost varies from $\mathcal{O}(rn)$ ($\mathbf{B}$ diagonal) to $\mathcal{O}(n\log n+rn)$ ($\mathbf{B}$ not diagonal), where $\mathcal{O}(n\log n)$ is the time for the factorisation of $\mathbf{L}$.

In spite of this different computational cost, the truncated spectral approximation (Fig. 4) is affected by small geometric undulations (especially at small scales), the use of heuristics for the selection of the number of Laplacian eigenpairs with respect to the target approximation accuracy, and the scale of features of the input shape. The spectrum-free computation (Fig. 5) generally provides better results in terms of smoothness, regularity, and accuracy of the computed spectral basis. For further comparison examples, we refer the reader to Sect. 4.

4 Discussion

We now discuss the selection of the filter map and main examples on the computation of the Laplacian spectral distances.

The spectral kernels and distances depend only on the behaviour of the filter in the spectral domain and on the Laplacian spectrum: increasing or decreasing the decay of the filter to zero encodes global or local shape details, respectively. Recalling that an arbitrary filter $1/\rho$ can be approximated in $\mathcal{R}^{l}_{r}$ with an accuracy of order $\mathcal{O}(s^{l+r+1})$ with respect to the $\ell_{\infty}$-norm (Sect. 3.2), we can use the space $\mathcal{R}^{l}_{r}$ to define any filter and the corresponding spectral distances. In this way, we reduce the degree of freedom in the definition of the filter to the selection of $(l+r+1)$ coefficients and without losing the richness of the resulting spectral distances. For the selection of the coefficients of the rational filter and the filter frequencies, we can apply the rules proposed in [24] for Laplacian spectral smoothing. In particular, the filter frequencies are derived from the dimension of a bounding box placed around the chosen feature $\mathcal{F}$ and whose axis are aligned with the eigenvectors of the covariance matrix of $\mathcal{F}$.

The truncated spectral approximation (c.f., Eq. (9)) of the diffusion, and more generally of the spectral distances, is affected by (i) small undulations (especially at small scales, Fig. 6), (ii) heuristics for the selection of the number of Laplacian eigenpairs with respect to the target approximation accuracy and the scale of the shape features, and (iii) the overall computational cost and storage overhead for the computation of a subpart of the Laplacian spectrum. The comparison of these results with the ones induced by the spectrum-free computation (Fig. 7) shows the improvement of the proposed approach in terms of smoothness, regularity, and accuracy of the computed distances.

At small scales (Fig. 8(e,f)), the truncated spectral approximation of the diffusion kernel generally have small negative values as a matter of the slow decay of the exponential filter to zero for small eigenvalues. For the Padé-Chebyshev approximation (Fig. 8(a-d)), the kernel values are positive at all the scales, as we approximate the filter with a higher accuracy with respect to the truncated spectral approximation, which applies a low-pass filter.

Finally, the analogous distribution and shape of the level-sets of the heat kernel and diffusion distance confirm the robustness of the Padé-Chebyshev of the approximation with respect to surface resolution (Fig. 7), partial sampling (Fig. 9), geometric (Fig. 10) and topological (Fig. 11) noise, almost isometric deformations (Fig. 12).

5 Conclusions and future work

We have presented a unified approach to the definition, discretisation, and computation of the spectral kernels and distances, which are defined by filtering the Laplacian spectrum and generalise the commute-time, bi-harmonic, diffusion, and wave kernel and distances. Even though analytic filters are easily defined and encode global/local shape properties, the corresponding spectral kernels and distances are not capable of characterising a specific shape class (e.g., humans with respect to four-legs animals, chairs with respect to tables). Indeed, the main topic for future research is a deeper analysis of the constraints on the filter in order to define “optimal” spectral kernels and distances for shape comparison. To achieve this goal, we plan to apply and specialise learning methods, as initially investigated in [2, 8]. Another interesting aspect is the study of spectral kernels for the adaptive hash retrieval [5] and the identification of those PDEs, whose Green kernels are useful for shape analysis, as already demonstrated by the bi-harmonic and diffusion kernels.

Acknowledgements

We thank the Reviewers for their thorough review and constructive comments, which helped us to improve the technical part and presentation of the revised paper. This work is supported by the H2020 ERC Advanced Grant CHANGE, grant agreement No. 694515. Shapes are courtesy of the AIM@SHAPE Repository, the SHREC 2010 and 2016 data sets; tet-meshes have been generated by the TETGEN software http://wias-berlin.de/software/tetgen/.