Hamiltonian operator for spectral shape analysis

by   Yoni Choukroun, et al.

Many shape analysis methods treat the geometry of an object as a metric space that can be captured by the Laplace-Beltrami operator. In this paper, we propose to adapt the classical Hamiltonian operator from quantum mechanics to the field of shape analysis. To this end we study the addition of a potential function to the Laplacian as a generator for dual spaces in which shape processing is performed. We present a general optimization approach for solving variational problems involving the basis defined by the Hamiltonian using perturbation theory for its eigenvectors. The suggested operator is shown to produce better functional spaces to operate with, as demonstrated on different shape analysis tasks.


Sparse Approximation of 3D Meshes using the Spectral Geometry of the Hamiltonian Operator

The discrete Laplace operator is ubiquitous in spectral shape analysis, ...

Steklov Spectral Geometry for Extrinsic Shape Analysis

We propose using the Dirichlet-to-Neumann operator as an extrinsic alter...

Shape-from-intrinsic operator

Shape-from-X is an important class of problems in the fields of geometry...

Steklov Geometry Processing: An Extrinsic Approach to Spectral Shape Analysis

We propose Steklov geometry processing, an extrinsic approach to spectra...

A preconditioned deepest descent algorithm for a class of optimization problems involving the p(x)-Laplacian operator

In this paper we are concerned with a class of optimization problems inv...

Why you should learn functional basis

Efficient and practical representation of geometric data is a ubiquitous...

Statistical analysis of random objects via metric measure Laplacians

In this paper, we consider a certain convolutional Laplacian for metric ...

1 Introduction

The field of shape analysis has been evolving rapidly during the last decades. The constant increase in computing power allowed image and shape understanding algorithms to efficiently handle difficult problems that could not have been practically addressed in the past. A large set of theoretical tools from metric and differential geometry, and spectral analysis has been imported and translated into action within the shape understanding arena. Among the numerous ways of analyzing shapes, a common one is to embed them into a different space where they can be processed more efficiently.

1.1 Related efforts

[Elad and Kimmel (2003)] introduced a method for analyzing surfaces based on embedding the intrinsic geometry of a given shape into a Euclidean space, extending previous efforts of [Schwartz et al. (1989), Zigelman et al. (2002), Grossmann et al. (2002)]. Their key idea was to consider a shape as a metric space, whose metric structure is defined by geodesic distances between pairs of points on the shape. Two non-rigid shapes are compared by first having their respective geometric structures mapped into a low-dimensional Euclidean space using multidimensional scaling (MDS) [Cox and Cox (2008)], and then comparing rigidly the resulting images, also called canonical forms.

[Mémoli and Sapiro (2005)] proposed a metric framework for non-rigid shape comparison based on the Gromov-Hausdorff distance that was suggested by Gromov as a theoretical tool to quantify disimilarity between metric spaces. Using the Gromov-Hausdorff formalism, the distance between two shapes is defined by matching pairwise distances on the shapes. However, the Gromov-Hausdorff distance is difficult to compute when treated in a straightforward manner. To overcome this difficulty [Bronstein et al. (2006a), Bronstein et al. (2006b)] proposed an efficient numerical solver based on a continuous optimization problem, known as Generalized MDS (GMDS). Recently, other relaxation schemes have been proposed, see for example [Chen and Koltun (2015), Aflalo et al. (2016)].

In the past decade, the Laplace-Beltrami operator (LBO) – the extension of the Laplacian to non-Euclidean manifolds, has become growingly popular. Its properties have been well studied in differential geometry and it was used extensively in computer graphics. The LBO can be found in countless applications such as mesh filtering [Vallet and Levy (2008)], mesh compression [Karni and Gotsman (2000)], shape retrieval [Bronstein et al. (2011)], to name just a few. It has been widely used in shape matching where several approaches treat the correspondence problem by comparing isometric invariant pointwise descriptors between the two shapes. For example, the Global Point Signature (GPS) [Rustamov (2007)], the Heat Kernel Signature (HKS) [Sun et al. (2009)] and the Wave Kernel Signature (WKS) [Aubry et al. (2011)]

, all use the eigenfunctions and eigenvalues of the LBO to compute local shape descriptors. Matching only signatures at a small set of points, the correspondence between the points on the two shapes can be found. These points can serve as anchors and interpolated for the entire shape

[Ovsjanikov et al. (2010)] where refinement of the basis can be performed to produce precise dense correspondence [Ovsjanikov et al. (2012), Pokrass et al. (2013), Shtern and Kimmel (2014a)].

Recently, learning based approaches [Litman and Bronstein (2014), Wei et al. (2016), Boscaini et al. (2016)] have also become highly popular in the shape matching arena.

The use of the basis defined by the LBO is in many senses a natural choice for surfaces analysis. It was chosen in the functional map framework [Ovsjanikov et al. (2012)] because of its compactness, stability, and invariance to isometries. Subsequently, it was proven to be optimal [Aflalo et al. (2015)] for representing smooth functions on the surface. In an attempt to overcome the topological sensitivity of the LBO and the non-local support of its eigenfunctions, compressed eigenfunctions have been adapted from mathematical physics to shape analysis [Neumann et al. (2014), Bronstein et al. (2016)]. Here, we try to find a richer family of basis functions that are based on intrinsic properties that can go beyond the geometry of the shape. Exploring a similar goal, [Kovnatsky et al. (2011)] combined geometric and photometric information within a unified metric for shape retrieval. [Iglesias and Kimmel (2012)] used artificial surface textures on shapes to define elliptic operators that give birth to a new family of diffusion distances. Along the same line of thought, [Hildebrandt et al. (2012)] designed a new family of eigenvibrations using extrinsic curvatures and deformation energies.

We suggest to further explore those ideas and construct the so-called potential operator that is added to the Laplace Beltrami operator. Here, a designed perturbation to the Laplacian permits a supervised control of the vibrational modes on the manifold.

1.2 Contributions

The main contribution of this paper is the exploration of the Hamiltonian operator on manifolds. We study spectral properties of the operator and the impact of an additional potential function to the Laplacian for shape analysis applications. The properties of the Hamiltonian allow it to be efficiently utilized by many spectral-based methods. The potential part can lead to a more descriptive operator when treated as a truncated basis generator. Modulated harmonics on the surface are obtained by treating different regions of interest as different values of the potential. We show that by simply plugging the resulting basis into existing spectral shape analysis pipelines could improve their performance.

The rest of the paper is organized as follows: in Section 2

, we propose to study the Hamiltonian on manifolds from the variational calculus point of view with motivation from quantum mechanics. We prove optimality properties of its eigenspace, characterize the associated diffusion process, the resulting nodal sets, introduce a discretization method and analyze the robustness of the operator.

In Section 3, we propose a global optimization framework for variational problems involving the basis defined by the Hamiltonian. We provide an approach for computing derivatives with respect to the potential based on eigenvectors perturbation theory. We demonstrate the effectiveness of the framework on the task of data representation.

In Section 4, we review recent improvement of the computation of the compressed modes [Ozolins et al. (2013)] that make use of the decomposition of the Hamiltonian. Finally, in Section 5 we present properties of the proposed basis that make it a better alternative for the task of shape matching where priors can be inserted through the potential in order to improve performance.

2 Hamiltonian operator

2.1 The Laplace Beltrami Operator

Consider a parametrized surface

with a metric tensor

. The space of square-integrable functions on is denoted by with the standard inner product , where is the area element induced by the Riemannian metric . The Laplace Beltrami Operator acting on a scalar function is defined as


where is the determinant of the metric matrix and is the inverse metric. If

is a domain in the Euclidean plane, the metric matrix is generally the identity matrix and the LBO reduces to the well-known Laplacian


The LBO is self-adjoint and thus admits a spectral decomposition , where and , such that,


with the Kronecker delta. In case has boundary, we add homogeneous Neumann boundary condition



is the normal vector to the boundary


The LBO eigendecomposition can be extracted from the Euler Lagrange solution to Dirichlet energy minimization


Here, each ordered eigenfunction composing the basis on the manifold corresponds to the function with the smallest possible energy that is orthogonal to all the previous ones. Therefore, the LBO eigenfunctions can be seen as an extension of the Fourier harmonics in Euclidean spaces to manifolds and are often referred to as Manifold Harmonics.

2.2 Hamiltonian

A Hamiltonian operator on a manifold acting on a scalar function , is an elliptic operator of the form


where is a real-valued scalar function. It plays a fundamental role in the field of quantum mechanics appearing in the famous Schrödinger equation that describes the wave motion of a particle with mass under potential ,


where is the Planck’s constant and represents the wave function of the particle such that

is interpreted as the probability distribution of finding the particle at a given position

at time .
The Schrödinger equation can be analyzed via perturbation theory by solving the spectral decomposition of the Hamiltonian


also known as the time-independent Schrödinger equation, where is the eigenenergy of a particle at the stationary eigenstate .

Since the potential is a diagonal operator, the Hamiltonian is self-adjoint as a sum of two self-adjoint operators and its eigenfunctions form a complete orthonormal basis on the manifold . As a generalization of the regular Laplacian, its spectral theory can be derived almost straightforwardly from that of the latter. Classical examples of the influence of potential functions in a one-dimensional Euclidean domain are depicted in Figure 1.

[width=1bb=0 0 1280 960,trim=150 50 110 150]potentials_euc.png

Figure 1: Influence of different potentials on the harmonics in one dimension.

2.3 Variational principle

Let us consider the following variational problem


whose the Euler-Lagrange equation defines the eigendecomposition of the Hamiltonian defined in (8).

The basis defined by the Hamiltonian operator corresponds to the orthogonal harmonics modulated by the potential function. The potential defines the trade-off between the orientation and the compactness of the basis and its global support. Larger values of the potential will enforce smooth solutions that concentrate on the low potential regions, while smaller ones will give solutions that better minimize the total energy at the expense of more extended wave functions.

2.4 Finite step potential

The time-independent Schrödinger equation can yield a rather complicated problem to solve analytically, even in one dimension. Let us consider a system with an ideal step potential in one dimension [Griffiths (2005)]. We need to solve the differential equation , with denoting the energy of the particle, and the Heaviside function with step of magnitude , at point , given by


The step divides the space in two constant-potential regions. At the zero potential region, the particle is free to move and the harmonic solutions are known. In the high potential region, on the other hand, for , the solution is a decaying exponentially, meaning that the particle cannot pass the potential barrier and is reflected according to classical physics. If , the solution is also harmonic, which means there is a probability for the particle to penetrate into the effective potential region with a different energy than that of particles in the zero potential region. We illustrate this effect in Figure 2 by numerically computing the eigenvectors of the LBO and the Hamiltonian with a potential defined on a human body surface.

[width=bb=0 0 1280 960,trim=0 0 0 0,clip]step_potential_2.png

Figure 2: Absolute values of the and eigenfunctions of the LBO (top). Absolute values of the corresponding eigenfunctions of the Hamiltonian with a step-function potential (bottom left), with step value . For this potential, the first eigenstate with energy greater than is the eighth. As analyzed, the eigenfunctions corresponding to lower eigenenergies are restricted to the region with , while the higher ones can have effective values (and oscillate) at the region. An evanescent wave can be observed at the seventh eigenstate.

Therefore, the potential energy can be tuned to enforce localization of the basis at the expense of loss of smoothness.

Let and be the spectral decompositions of the Laplacian, and the Hamiltonian, respectively. Then, everywhere on the manifold implies that the eigenvalues satisfy


Proof. According to the Courant-Fischer min-max theorem, we have




Since the family of eigenvalues of the Helmholtz equation (3) consist of a diverging sequence ( [Weyl et al. (1950)]), there exists an such that

and the trade-off between local-compact and global support of the basis elements can be controlled by the potential energy. Then, we can estimate the magnitude of the potential required in order to allow for oscillations outside the regions where the potential vanishes.

Given a scalar we can define the Hamiltonian as


where controls the resistance to diffusion induced by the potential. Let and be the -th eigenvalue of the LBO and Hamiltonian, respectively. We seek for a constant such that so the particle can penetrate the high potential region. Considering the potential as small perturbation of the Laplacian, up to first order, the eigenenergies are defined as . In order to contain the basis support at most until the -th eigenfunction, must satisfy


According to its potential energy, the basis can then provide supervised multiresolution analysis on the manifold by containing the first eigenfunctions and allow global analysis for the following.

2.5 Optimality of the Hamiltonian eigenspace

Let us consider a function . We define the representation residual function as


Defining , we know that


Thus, from (20) and (21) we obtain


Recall that for we return to the LBO case. Among the numerous reasons that motivated the selection of the Laplacian for shape analysis, a major one is its efficiency in representing functions with bounded gradient magnitude. This result was subsequently proved to be optimal for representing functions with bounded gradient magnitude over surfaces in [Aflalo et al. (2015)], which says that there exists no other basis with better representation error for all possible functions.

In case of the Hamiltonian, the Dirichlet energy is coupled with the potential energy. Thus the Hamiltonian operator advocates measuring smoothness differently for different regions of the domain where smoothness remains a less important factor than avoiding vibrations in high potential areas. This is a useful property to exploit in different shape analysis scenarios.

Next, we show that the Hamiltonian is optimal in approximating functions with both bounded gradient and low values in high potential areas. Let . There is no integer and no sequence of linearly independent functions in such that


The proof of Theorem 2.5 is given in the Appendix.

2.6 Diffusion process

Let us be given a Riemannian manifold . The heat equation governing the diffusion process on is defined as


with appropriate boundary conditions. A natural extension to the new operator with a potential , can be written as


The solutions of (24) and (25) have the form [Iglesias and Kimmel (2012)]


that represents the diffusion in time of heat on the manifold with potential , where . We refer to as the heat kernel. A standard proof is given in the Appendix.

According to the Feynman-Kac formula [Simon (2005)], the solution of the diffusion process is expressed in terms of the Wiener process,


In the Laplacian case, the initial value is carried over random paths in time, while the expected value of the stochastic process is equal to the solution . For , the diffusion spreads according to the potential on the manifold, when the transported value is modulated exponentially by the potential , diffusing anisotropically to low potential regions, as shown in Figure 3.

[width=trim=200 0 0 0,bb=0 0 1500 960]Heat_Kernel.png V

Figure 3: Heat diffusion with a delta function at the centaur’s head as initial condition. The diffusion is derived from the LBO (top) and the Hamiltonian (bottom) for different values of . The potential used in this example is the geodesic distance from the front left leg. A signature extracted from a diffusion process using the Hamiltonian is more descriptive and in this case allows to resolve ambiguities due to symmetry.

2.7 Nodal sets

An interesting property of the Laplacian is the relation between its eigenfunctions, the number of connected nodal (zero) sets, and the number of complementary regions they define. Given an eigenfunction , a nodal set is defined as the set of points at which the eigenfunction values are zero. That is,


The Nodal Theorem [Courant and Hilbert (1966)] states that the -th eigenfunction of the LBO can split to at most connected sub-domains. In other words, the zero set of the -th eigenfunction can separate the manifold into at most connected components.

Proposition 1

Given the self-adjoint Hamiltonian operator on , with arbitrary boundary conditions; if its eigenfunctions are ordered according to increasing eigenvalues, then, the nodal set of the -th eigenfunction divides the domain into no more than connected sub-domains.

The proof is essentially the same as that of the Laplacian case. See [Courant and Hilbert (1966)] Vol.1 Sec. VI.6 for a proof.

As shown in Figure 2, the Hamiltonian eigenfunctions are tuned by the potential. Thus, shape segmentation can be obtained by separating the surface according to the induced nodal sets as described in [Levy (2006)]. Given a potential defined on the surface, meaningful segmentation can be induced by the nodal domains of the resulting eigenfunctions, as presented in Figure 4.

[width=bb=0 0 1280 960]Nodal_Sets_2.png VV

Figure 4: Nodal domains obtained from the nodal sets of the Hamiltonian (second and third columns) and the LBO (fourth and fifth columns). Two potentials are depicted (first and second rows). One can observe a meaningful segmentation induced by the nodal sets of the Hamiltonian.

2.8 Discretization

In the discrete setting, we consider a triangular mesh in with the associated space of functions that are continuous and linear in every triangle. According to the Finite Element Method (FEM) [Dziuk (1988)], the solution of the Hamiltonian eigenvalue problem (8) can be computed by imposing that the equation is satisfied in a weak sense, that is,


where denote the Lagrange basis of piecewiselinear hat-functions on , that take the value one at a vertex and vanishes at all the other vertices. represents the eigenenergy of the Hamiltonian. On the mesh, the bilinear form is evaluated by splitting the integrals into a sum over the triangles of by


Since the Hamiltonian is a linear operator we have


The matrix representation of and with respect to the Lagrange basis are well known [Pinkall and Polthier (1993)] and define the stiffness matrix and the mass matrix with the entries




where the last equality is obtained by representing the potential function as a diagonal matrix according to the Lagrange basis functions. The discretization of the eigenvalue problem (8) is defined by finding all pairs such that


Efficient solution methods can be found in [Vallet and Levy (2008)]. Among the possible explicit representations of the matrices and , we use here the cotangent formula [Pinkall and Polthier (1993), Meyer et al. (2003)] where the stiffness matrix is defined as


with , where is the set of edges of the triangulated surface interpreted as a graph and denote the angles and of the triangles sharing the edge . The mass matrix is replaced by a diagonal lumped mass matrix of the area of local mixed Voronoi cells about each vertex . The manifold inner product is discretized as . Since only modifies the diagonal of , our operator remains a sparse matrix with the same effective entries, and thus, there is no increase in the computational cost of the generalized eigendecomposition compared to that of the LBO.

2.9 Robustness to noise

As a generalization of the Laplacian, the Hamiltonian exhibits similar robustness to noise. Consider the Hamiltonian matrix with the potential. Then, the perturbed Hamiltonian has the form . Let us define and . Based on perturbation theory, and up to second-order corrections, the -th eigenfunction of has the form


with and being respectively the

-th eigenfunction and eigenvalue of the unperturbed Hamiltonian. Assuming uniformly distributed random noise on the mesh, the eigenfunctions of the regular Laplacian may present smaller distortion to noise than the Hamiltonian since the perturbation is amplified by area and potential distortions. Still, in case of potential with small values the distortion is insignificant. In Figure

5, we present the original surface and its noisy version in which vertex positions have been corrupted by additive Gaussian noise with of the mean edge length. The potential is also modified by adding a Gaussian noise with

of the initial variance of the potential.

[width=1trim=0 350 0 0,clip]hands_hole_Noise

Figure 5: Robustness to noise of the Hamiltonian. First eigenfunctions of the Hamiltonian under potential (top). First eigenfunctions of the Hamiltonian subject to Gaussian noise in positions of the vertices and the potential (middle). First eigenfunctions of the Hamiltonian subject topological noise (bottom).

The construction of the Laplacian depends crucially on the mesh connectivity making it sensitive to topological noise such as holes and part removal that can be found in many depth acquisition scenarios. The compact support of the basis elements of the Hamiltonian makes it robust to noise compared to the basis elements that are generated by the Laplacian. We illustrate the robustness property in Figure 5 where of the surface area was removed due to topological noise in the form of small holes.

3 Optimization of the potential

One natural problem emerging when working with the Hamiltonian is the ability to define an optimal potential function for a specific task. The choice of the potential is application dependent but can be represented through minimization problem generically defined as


where denotes the data term depending on the data matrix and the vector defining the diagonal potential matrix. Regularization terms can be further be added. If the analytical solution remains complex, a common approach is to minimize the goal function with an optimization algorithm involving the gradient of the goal function with respect to the potential. In this section we propose an optimization framework based on perturbation theory of the eigenvectors where optimal potential is obtained. To that end, we need to derive the gradient for a given objective .

Here we will consider the problem of data representation using the discrete basis of the Hamiltonian referred to as representing the eigenvectors of the Hamiltonian such that . The discretized minimization problem is defined as


with . The objective defines the representation error of the data in the subspace spanned by the columns of and in the sense of the Frobenius norm . For a general orthonormal matrix

, the problem is equivalent to Principal Component Analysis (PCA). We can straightforwardly obtain that


Thus, the differential of the loss function

with respect to is obtained by


It remains to derive the differential of the Hamiltonian eigenvectors. Let us consider the full matrix of eigenvectors , the diagonal matrix of eigenenergies and the discrete Hamiltonian operator . The eigenvalue decomposition problem is given by . Thus, the differential of the spectral decomposition problem is given by


Multiplying by on the left side and denoting with , we have


since . We readily obtain that the off diagonal elements of the matrix can be defined by


Here represents the -th column of the matrix of eigenvectors. The diagonal elements of are defined by the following


The diagonal elements are then defined by . Since second order elements are negligible, we have . We obtain that , with denoting the Hadamard product and the matrix defined as


The selection of the first eigenvectors are obtained by multiplying by the truncated identity matrix . The differential is now known and can be plugged into (40) in order to extract , that is:


The passage in the fourth line stems from the equivalence . Since , we obtain finally


Two problems arise from the suggested scheme. First, the high computational cost of a full (sparse) matrix diagonalization. Second, the matrix remains undefined when eigenvectors have non-trivial multiplicities. The first problem can be relaxed by approximating the matrix with less eigenvectors. This is especially justified for distant indices, where the eigenenergies are well separated and the corresponding elements of matrix become negligible. Also, the data can be projected onto the LBO basis so the solution complexity remains constant with the size of the mesh. Even if the second problem has been treated in [Van Der Aa et al. (2007)], it seems that lack of smoothness at isolated points is not critical for computation and convergence by resorting to a sub-gradient approach. The alternative opted for here is to stabilize the matrix in order to avoid exploding gradients. We use the approximation


where the sign function is not vanishing.

In its geometric setting, the solution remains similar and is obtained by defining a new basis , such that , coupled with the consistently discretized Frobenius norm . In the following experiments we allowed negative potential for performance consideration only, since the potential is defined over the whole codomain . Also, for physically interpretable solutions we enforced positive potential by using quadratic function . The extension of the derivation is straightforward but decreased the performance since it is more restrictive.

As a toy experiment, we propose to find the best potential for the representation of a function in the one dimensional Euclidean domain. Given a function , we seek for the best potential minimizing . We compare in Figure 6 the reconstruction performance on a one dimensional linear function with the Laplacian and the Hamiltonian built from the optimized potential.

[width=1trim=90 30 90 0,clip]euc1

Figure 6: Reconstruction of a linear function using the Laplacian and the Hamiltonian constructed with the proposed framework. 15 eigenvectors were used in this experiment. Observe that the potential is high close to the boundary to reduce the representation error.

In Figure 7, we propose to reconstruct the matrix of coordinates of a mesh so the data matrix is defined by . The experiments were implemented using the quasi-Newton method with initial zero potential.

[width=1bb=0 0 1280 960,trim=0 0 0 0,clip]recon_PCA.png LBOHamiltonian

Figure 7: Potential function defined on the original mesh (left), reconstruction of the mesh coordinates with 50 eigenvectors using the LBO (middle) and the Hamiltonian constructed with the proposed method (right). Blue and red colors represent negative and positive values respectively. The Hamiltonian is able to focus on sharp regions of the mesh designated by the blue regions of the potential for a better reconstruction (fingers). The errors are 0.0015 and 0.00061 for the LBO and the Hamiltonian respectively.

An important application related data representation is spectral mesh compression. [Karni and Gotsman (2000)] proposed to project the coordinates functions of the mesh onto the LBO eigenfunctions in order to encode the mesh geometry via the first coefficients only. Since most of the function energy is generally contained in the first coefficients, the reconstruction distortion is low, up to fine details related to higher frequencies. Since matrix decomposition is an expensive operation, they suggested to segment the shape into smaller parts that can be processed separately. By sending the mesh topology (triangles) separately, the combinatorial graph Laplacian is built on the decoder side and the signal can be reconstructed with the received coefficients. We suggest to apply this idea to our basis which potential is obtained by the proposed optimization framework. However, one major drawback is that we need to encode the potential as well as the coefficients. Also, some methods use the ordering of the vertices in order to encode information [Touma and Gotsman (1998)]. Here we suggest to reorder the vertices such that the vertex with the smallest potential is be assigned the index and the vertex with the largest potential is be assigned the index . By using a fixed potential defined as


the decoder simply applies in order to obtain the Hamiltonian basis. Here and are the regression coefficients minimizing that are also encoded. We present compression results in Figure 8.

[width=1trim=40 0 20 0,clip]comp_fandisk

Figure 8: Geometry compression performance comparison between the Laplacian (MHB), the proposed projected operator (H) and the optimal Hamiltonian (H opt) using the proposed framework for the sharp Fandisk shape. The optimal Hamiltonian performance are presented with no encoding of the potential itself.

4 Compressed Manifold Modes

[Ozolins et al. (2013)] proposed a a novel method to create a set of localized eigenfunctions in Euclidean domains. To that end, they modified the construction of standard differential operators by adding an regularization term to the variational leading to the decomposition of the operator. The resulting eigenfunctions were called compressed modes and were shown to be compactly supported [Brezis (1974)]. [Neumann et al. (2014)] extended this construction to manifolds, suggesting the following discrete regularization problem


with the parameter that controlling the localization of the basis. Proposed solutions require the use of expensive and unstable optimization techniques ([Neumann et al. (2014), Kovnatsky et al. (2016)]) based on ADMM and proximal operators.

The latter optimization problem (50) can be written as an Hamiltonian eigendecomposition problem [Bronstein et al. (2016)]


where is the diagonal matrix operator defining the potential that corresponds to the -th eigenvector that localizes the support of in low-potential areas. The potential is defined iteratively using a reweighted least squares scheme


ensuring that the minimizers of (50) and (51) coincide. Interestingly, the potential here is defined as a function of the eigenfunction, namely . The potential and the resulting eigenstate are then intrinsically linked, meaning that the potential is influenced by the state of the particle itself. Consequently, a perturbation of the potential enforces perturbation of the eigenfunction and vice versa until reaching steady state.

We formulate the compressed manifold modes problem as


with and where is a sufficiently large constant such that the third term guarantees that the -th mode is -orthogonal to the previously computed modes , . Observe that albeit non-convex, the problem has a closed form global solution, that is the smallest generalized eigenvector of



For small number of compressed modes, is a low rank matrix and finding the smallest generalized eigenvector can be solved efficiently since the involved matrix is the sum of a sparse and a low-rank matrix.

Several numerical eigendecomposition implementations use the Arnoldi iteration algorithm. In our matrix decomposition problem, the core operation is the multiplication by the inverse of the matrix with a vector, operation that cannot be solved straightforwardly. Also, shifting the maximum eigenvalue using power method is too unstable since it depend on gap of the first eigenvalues, generally tight. In our configuration the Woodbury identity

can be used to compute efficiently the vector multiplication with the inverse of the matrix as a cascade of sparse and low-rank systems. Unlike solutions of the inconsistently discretized problem 50, the basis obtained with the proposed Hamiltonian framework is more robust under various discretizations and can be computed at a fraction of the computational cost (Figure 9).

Figure 9: Runtimes of Neumann et al. and the proposed framework on meshes of varying size (number of vertices ) and number of eigenvectors

. Averages and standard deviations are presented over

runs. Same stopping criteria were applied to all methods.

5 Shape matching

The task of matching pairs of shapes lies at the core of many shape analysis tasks and plays a central role in operations such as 3D alignment and shape reconstruction. While rigid shape matching has been well studied in the literature, non-rigid correspondence remains a difficult task even for nearly isometric surfaces. When dealing with rigid objects, it is sufficient to find the rotation and translation that aligns one shape to the other [Chen and Medioni (1992)]

. Therefore, the rigid matching problem amounts to determining only six degrees of freedom. At the other end, non-rigid matching generally requires dealing with many more degrees of freedom. Since the LBO is invariant to isometric deformations, it has been used extensively to aid the solution of correspondence problem. Several properties of the Hamiltonian operator make it a better choice for this task compared to its zero-potential particular case that is the LBO.

Invariance. The Laplace-Beltrami Operator is defined in terms of the metric tensor which is invariant to isometries. For a potential function defined intrinsically, the resulting Hamiltonian is also isometry-invariant.

Compactness. Compactness means that scalar functions on a shape should be well approximated by using only a small number of basis elements. From Theorem 2.5 and as a generalization of the Laplacian, the global support and compactness hold for a bounded (low) potential.

Descriptiveness. The LBO eigenvalues are related to frequency. Similarly, eigenenergies of the Hamiltonian relate to the number of oscillations on the manifold. Theorem 2.4 demonstrates that the modes corresponding to small eigenvalues of the Hamiltonian defined with a positive potential, encapsulate higher frequencies, even when localized, compared to the modes of the regular LBO. At the other end, highly oscillating eigenfunctions can be used to represent fine details of the shape that can be crucial for shape matching. Also, the potential enforces different oscillations in different regions on the manifold, allowing for better discrimination of similar areas and disambiguation of intrinsic symmetries with asymmetric potential.

Stability. Deformations of non-rigid shapes and articulated objects can stretch the surface. In such cases, the LBO eigendecomposition of the two shapes will be different. We could compensate for such local metric distortions by carefully designing a potential. Assigning high potential to strongly distorted regions would lead to lower values of the eigenfunctions in those areas (9). Such a potential will reduce the discrepancy between corresponding eigenfunctions at least for the lower eigenergies, as shown via the functional maps representations [Ovsjanikov et al. (2012)] in Figure 10. Let define and , the area at vertex on mesh and on the second mesh respectively and a bijection between two (discretized) surfaces and . Then, we define the potential at vertex as


(a) Nearly isometric shapes

\put(-55.0,75.0){\scriptsize{{ }}}
(b) LBO

\put(-55.0,75.0){\scriptsize{{ }}}
(c) Hamiltonian
Figure 10: Two nearly isometric meshes with high potential (hot colors) in large distortion regions (a), functional maps matrix of the LBO (b) and the Hamiltonian (c).

Among the few stable intrinsic invariants that can be extracted from the geometry, we will use the stable first eigenfunctions of the LBO and geodesic distances. Additional non necessarily intrinsic information such as photometric properties or even extrinsic shape properties such as principal curvatures [Hildebrandt et al. (2012)] can also be integrated into the potential field.

5.1 Experimental Evaluation

We tested the proposed basis and compared its matching performances to that the LBO basis as applied to pairs of triangulated meshes of shapes from the TOSCA dataset [Bronstein et al. (2008)] and the SCAPE dataset [Anguelov et al. (2004)]. The TOSCA data set contains densely sampled synthetic human and animal surfaces, divided into several classes with given ground-truth point-to-point correspondences between the shapes within each class. The SCAPE data set contains scans of real human bodies in different poses. The evaluation method used is described in [Kim et al. (2011)] where the distortion curves describe the percentage of surface points falling within a relative geodesic distance from what is assumed to be their true locations. Symmetries were not allowed in all evaluations. Note that we assume that the sign ambiguity of the first eigenfunctions generating the potential is resolved [Shtern and Kimmel (2014b)].

Figure 11 compares the two operators by matching diffusion kernel descriptors derived from the corresponding eigenfunctions. The diffusion on the shape using the Hamiltonian as the diffusion operator is more descriptive than regular diffusion that cannot resolve the symmetries. Also, it would be natural to compute the WKS signature when the Schrödinger equation is governed by a given effective potential. As intrinsic positive potential we use the normalized sum of the four first nontrivial eigenfunctions of the LBO on each shape, adding a constant of minimal value in order to obtain a non-negative potential. This way only the intrinsic unstable geometry of the shape is involved in defining the Hamiltonian operator.

\put(-75.0,80.0){\scriptsize{{ }}}

\put(-75.0,80.0){\scriptsize{{ }}}
Figure 11: Evaluation of the diffusion kernels signatures matches on the TOSCA and SCAPE datasets.

In case we know which regions are prone to elastic distortions, like joints and stretchable skin in articulated objects, we could suppress the effect of those regions in our matching procedures by using an appropriate potential as a selective mask. Figure 12, compares the operator with and without potential by matching the spectral signatures computed by the framework of [Shtern and Kimmel (2015)]. The potential we used is the local area distortion when comparing the meshes of two corresponding objects, as in (55). The descriptiveness of the potential and the localization of the harmonics lead to more accurate matching results.


Figure 12: Evaluation of the spectral signature matches on the TOSCA and SCAPE data-sets.

To investigate the performances of the Hamiltonian with photometric textures used as potential, we present in Figure 13 the results of different signatures matching with a dalmatian texture defined for the ”Dogs” shapes from the TOSCA data set.

[scale=0.2,trim=250 60 250 0,clip]dog0RealDalmatianTextureYoniBlender

(a) Photometric data

[scale=0.21,trim=0 0 0 0,clip]dog_dalmatian

(b) Signatures
Figure 13: Evaluation of the descriptors matches on the ”Dogs” benchmark from the TOSCA dataset with dalmatian texture.

Iterative refinement of functional representations have been proven to be powerful in shape matching [Ovsjanikov et al. (2012)]. Given an initial partial or dense map, it tries to recover iteratively dense and accurate matching between two given shapes. Here we use a similar refinement framework dubbed as Iterative Closest Spectral Kernel Maps (ICSKM) [Shtern and Kimmel (2014a)] for performance comparison between the two bases. Figure 14 compares the regular ICSKM algorithm working with the Laplacian eigenspace and the Hamiltonian method when we provided one, two, or three landmark points, that were randomly selected from the ground-truth mapping. The potential used in these examples is the geodesic distance from the landmark points. Note that again we use only the geometry of the shapes in order to refine the match between them using the new basis.


\put(13.0,120.0){\scriptsize{{ }}}

Figure 14: Evaluation of the ICSKM algorithm with different landmark initialization matches on the TOSCA dataset. We used geodesic distances from given landmark points as intrinsic geometric potential on the shapes.

6 Conclusion

A classical operator was adopted from the field of quantum mechanics and adapted to shape analysis problems. Functional and spectral properties of the Hamiltonian operator were presented and compared to the popular Laplacian operator often used in many shape analysis procedures. A general optimization method for solving variational problems involving the Hamiltonian operator have been proposed and employed to the task of mesh compression. Features and texture properties can be incorporated into the new operator to obtain a descriptive and stable basis that provides a powerful domain of operation for shape matching. Various directions for future research include exploration of the operator on other shape analysis tasks such as partial shape matching where occluded areas could be refined via the potential. This work has been supported by Grant agreement No. 267414 of the European Community’s FP7-ERC program.


  • Aflalo et al. (2015) Aflalo, Y., Brezis, H., and Kimmel, R. 2015. On the optimality of shape and data representation in the spectral domain. SIAM J. Imaging Sciences 8, 2, 1141–1160.
  • Aflalo et al. (2016) Aflalo, Y., Dubrovina, A., and Kimmel, R. 2016. Spectral generalized multi-dimensional scaling.

    International Journal of Computer Vision (IJCV)

     118, 3, 380–392.
  • Anguelov et al. (2004) Anguelov, D., Srinivasan, P., Pang, H.-C., Koller, D., Thrun, S., and Davis, J. 2004. The correlated correspondence algorithm for unsupervised registration of nonrigid surfaces. In Proceedings of the 17th International Conference on Neural Information Processing Systems. NIPS’04. MIT Press, Cambridge, MA, USA, 33–40.
  • Aubry et al. (2011) Aubry, M., Schlickewei, U., and Cremers, D. 2011. The wave kernel signature: A quantum mechanical approach to shape analysis. In Computer Vision Workshops (ICCV Workshops), 2011 IEEE International Conference on. 1626–1633.
  • Boscaini et al. (2016) Boscaini, D., Masci, J., Rodolà, E., and Bronstein, M. M. 2016. Learning shape correspondence with anisotropic convolutional neural networks. Tech. Rep. arXiv:1605.06437.
  • Brezis (1974) Brezis, H. 1974. Solutions with compact support of variational inequalities. Russian Mathematical Surveys 29, 2, 103–108.
  • Brezis (2010) Brezis, H. 2010.

    Functional analysis, Sobolev spaces and partial differential equations

    Springer Science & Business Media.
  • Bronstein et al. (2008) Bronstein, A., Bronstein, M., and Kimmel, R. 2008. Numerical geometry of non-rigid shapes, 1 ed. Springer Publishing Company, Incorporated.
  • Bronstein et al. (2016) Bronstein, A., Choukroun, Y., Kimmel, R., and Sela, M. 2016. ”Consistent discretization and minimization of the L1 norm on manifolds”. In 3D Vision (3DV), 2016 4nd International Conference on. IEEE.
  • Bronstein et al. (2011) Bronstein, A. M., Bronstein, M. M., Guibas, L. J., and Ovsjanikov, M. 2011. Shape google: Geometric words and expressions for invariant shape retrieval. ACM Trans. Graph. 30, 1 (Feb.), 1:1–1:20.
  • Bronstein et al. (2006a) Bronstein, A. M., Bronstein, M. M., and Kimmel, R. 2006a. Efficient computation of isometry-invariant distances between surfaces. SIAM J. Scientific Computing 28, 5, 1812–1836.
  • Bronstein et al. (2006b) Bronstein, A. M., Bronstein, M. M., and Kimmel, R. 2006b. Generalized multidimensional scaling: A framework for isometry-invariant partial surface matching. Proceedings of the National Academy of Sciences 103, 5, 1168–1172.
  • Chen and Koltun (2015) Chen, Q. and Koltun, V. 2015. Robust nonrigid registration by convex optimization. In 2015 IEEE International Conference on Computer Vision (ICCV). 2039–2047.
  • Chen and Medioni (1992) Chen, Y. and Medioni, G. 1992. Object modelling by registration of multiple range images. Image Vision Comput. 10, 3 (Apr.), 145–155.
  • Courant and Hilbert (1966) Courant, R. and Hilbert, D. 1966. Methods of mathematical physics. Vol. 1. CUP Archive.
  • Cox and Cox (2008) Cox, A. M. A. and Cox, F. T. 2008. Multidimensional Scaling. Springer Berlin Heidelberg, Berlin, Heidelberg, 315–347.
  • Dziuk (1988) Dziuk, G. 1988. Finite elements for the beltrami operator on arbitrary surfaces. In Partial differential equations and calculus of variations. Springer, 142–155.
  • Elad and Kimmel (2003) Elad, A. and Kimmel, R. 2003. On bending invariant signatures for surfaces. IEEE Trans. Pattern Anal. Mach. Intell. 25, 10 (Oct.), 1285–1295.
  • Griffiths (2005) Griffiths, D. 2005. Introduction to Quantum Mechanics. Pearson international edition. Pearson Prentice Hall.
  • Grossmann et al. (2002) Grossmann, R., Kiryati, N., and Kimmel, R. 2002. Computational surface flattening: a voxel-based approach. IEEE Transactions on Pattern Analysis and Machine Intelligence 24, 4 (Apr), 433–441.
  • Hildebrandt et al. (2012) Hildebrandt, K., Schulz, C., von Tycowicz, C., and Polthier, K. 2012. Modal shape analysis beyond laplacian. Computer Aided Geometric Design 29, 5, 204–218.
  • Iglesias and Kimmel (2012) Iglesias, J. A. and Kimmel, R. 2012. Schrödinger Diffusion for Shape Analysis with Texture. Springer Berlin Heidelberg, Berlin, Heidelberg, 123–132.
  • Karni and Gotsman (2000) Karni, Z. and Gotsman, C. 2000. Spectral compression of mesh geometry. In Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques. SIGGRAPH ’00. ACM Press/Addison-Wesley Publishing Co., New York, NY, USA, 279–286.
  • Kim et al. (2011) Kim, V. G., Lipman, Y., and Funkhouser, T. 2011. Blended intrinsic maps. ACM Trans. Graph. 30, 4 (July), 79:1–79:12.
  • Kovnatsky et al. (2011) Kovnatsky, A., Bronstein, M. M., Bronstein, A. M., and Kimmel, R. 2011. Photometric heat kernel signatures. In Scale Space and Variational Methods in Computer Vision - Third International Conference, SSVM 2011, Ein-Gedi, Israel, May 29 - June 2, 2011, Revised Selected Papers. 616–627.
  • Kovnatsky et al. (2016) Kovnatsky, A., Glashoff, K., and Bronstein, M. M. 2016. MADMM: a generic algorithm for non-smooth optimization on manifolds. In European Conference on Computer Vision. Springer, 680–696.
  • Levy (2006) Levy, B. 2006. Laplace-beltrami eigenfunctions towards an algorithm that ”understands” geometry. In Proceedings of the IEEE International Conference on Shape Modeling and Applications 2006. SMI ’06. IEEE Computer Society, Washington, DC, USA, 13–.
  • Litman and Bronstein (2014) Litman, R. and Bronstein, A. M. 2014. Learning spectral descriptors for deformable shape correspondence. IEEE transactions on pattern analysis and machine intelligence 36, 1, 171–180.
  • Mémoli and Sapiro (2005) Mémoli, F. and Sapiro, G. 2005. A theoretical and computational framework for isometry invariant recognition of point cloud data. Foundations of Computational Mathematics 5, 3, 313–347.
  • Meyer et al. (2003) Meyer, M., Desbrun, M., Schröder, P., and Barr, A. H. 2003. Discrete differential-geometry operators for triangulated 2-manifolds. Springer Berlin Heidelberg, Berlin, Heidelberg, 35–57.
  • Neumann et al. (2014) Neumann, T., Varanasi, K., Theobalt, C., Magnor, M., and Wacker, M. 2014. Compressed manifold modes for mesh processing. Comput. Graph. Forum 33, 5 (Aug.), 35–44.
  • Ovsjanikov et al. (2012) Ovsjanikov, M., Ben-Chen, M., Solomon, J., Butscher, A., and Guibas, L. 2012. Functional maps: A flexible representation of maps between shapes. ACM Trans. Graph. 31, 4 (July), 30:1–30:11.
  • Ovsjanikov et al. (2010) Ovsjanikov, M., Mérigot, Q., Mémoli, F., and Guibas, L. 2010. One point isometric matching with the heat kernel. Computer Graphics Forum 29, 5, 1555–1564.
  • Ozolins et al. (2013) Ozolins, V., Lai, R., Caflisch, R., and Osher, S. 2013. Compressed modes for variational problems in mathematics and physics. Proceedings of the National Academy of Sciences 110, 46, 18368–18373.
  • Pinkall and Polthier (1993) Pinkall, U. and Polthier, K. 1993. Computing discrete minimal surfaces and their conjugates. Experiment. Math. 2, 1, 15–36.
  • Pokrass et al. (2013) Pokrass, J., Bronstein, A. M., Bronstein, M. M., Sprechmann, P., and Sapiro, G. 2013. Sparse modeling of intrinsic correspondences. Computer Graphics Forum 32, 2pt4, 459–468.
  • Rustamov (2007) Rustamov, R. M. 2007. Laplace-beltrami eigenfunctions for deformation invariant shape representation. In Proceedings of the Fifth Eurographics Symposium on Geometry Processing. SGP ’07. Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, 225–233.
  • Schwartz et al. (1989) Schwartz, E. L., Shaw, A., and Wolfson, E. 1989. A numerical solution to the generalized mapmaker’s problem: Flattening nonconvex polyhedral surfaces. IEEE Trans. Pattern Anal. Mach. Intell. 11, 9 (Sept.), 1005–1008.
  • Shtern and Kimmel (2014a) Shtern, A. and Kimmel, R. 2014a. Iterative closest spectral kernel maps. In 3D Vision (3DV), 2014 2nd International Conference on. Vol. 1. IEEE, 499–505.
  • Shtern and Kimmel (2014b) Shtern, A. and Kimmel, R. 2014b. Matching the LBO eigenspace of non-rigid shapes via high order statistics. Axioms 3, 3, 300–319.
  • Shtern and Kimmel (2015) Shtern, A. and Kimmel, R. 2015. Spectral gradient fields embedding for nonrigid shape matching. Computer Vision and Image Understanding 140, 21 – 29.
  • Simon (2005) Simon, B. 2005. Functional integration and quantum physics. 2nd ed., 2nd ed. ed. Providence, RI: AMS Chelsea Publishing.
  • Sun et al. (2009) Sun, J., Ovsjanikov, M., and Guibas, L. 2009. A concise and provably informative multi-scale signature based on heat diffusion. In Proceedings of the Symposium on Geometry Processing. SGP ’09. Eurographics Association, Aire-la-Ville, Switzerland, Switzerland, 1383–1392.
  • Touma and Gotsman (1998) Touma, C. and Gotsman, C. 1998. Triangle mesh compression. PROC GRAPHICS INTERFACE. pp. 26-34. 1998.
  • Vallet and Levy (2008) Vallet, B. and Levy, B. 2008. Spectral Geometry Processing with Manifold Harmonics. Computer Graphics Forum.
  • Van Der Aa et al. (2007) Van Der Aa, N., Ter Morsche, H., and Mattheij, R. 2007. Computation of eigenvalue and eigenvector derivatives for a general complex-valued eigensystem. Electronic Journal of Linear Algebra 16, 1, 300–314.
  • Wei et al. (2016) Wei, L., Huang, Q., Ceylan, D., Vouga, E., and Li, H. 2016. Dense human body correspondences using convolutional networks. In

    Computer Vision and Pattern Recognition (CVPR)

  • Weyl et al. (1950) Weyl, H. et al. 1950. Ramifications, old and new, of the eigenvalue problem. Bulletin of the American Mathematical Society 56, 2, 115–139.
  • Zigelman et al. (2002) Zigelman, G., Kimmel, R., and Kiryati, N. 2002. Texture mapping using surface flattening via multidimensional scaling. IEEE Transactions on Visualization and Computer Graphics 8, 2 (Apr), 198–207.

Appendix A Proof of Theorem 2.

Let us be given the Hamiltonian operator .

Recall the Courant–Fischer min-max principle; see also [Aflalo et al. (2015)] and [Brezis (2010)] Problems 37 and 49. We have for every ,


That is, the min is taken over a linear subspace with is the Sobolev space of co-dimension and the max is taken over all such subspaces.
Set , so that is a subspace of co-dimension .
By 56 we have that for all


and thus


On the other hand, by 56,


Combining 58 and 59 yields .                                              

Appendix B Diffusion kernel of the Hamiltonian

In order to solve the diffusion equation, we first need to find the fundamental solution kernel to the Dirichlet problem that yields the heat equation


Recall that for we return to the regular LBO diffusion case.
Suppose that has a eigendecomposition . In that case, we can write


and from the linearity of we have