We consider the problem of non-rigid shape matching using the functional map
framework. Specifically, we analyze a commonly used approach for regularizing functional maps,
which consists in penalizing the failure of the unknown map to commute with the Laplace-Beltrami
operators on the source and target shapes. We show that this approach has certain undesirable
fundamental theoretical limitations, and can be undefined even for trivial maps
in the smooth setting. Instead we propose a novel, theoretically well-justified approach for
regularizing functional maps, by using the notion of the resolvent of the Laplacian
operator. In addition, we provide a natural one-parameter family of regularizers, that can be
easily tuned depending on the expected approximate isometry of the input shape pair. We show on a
wide range of shape correspondence scenarios that our novel regularization leads to an improvement in the quality of the estimated functional, and ultimately pointwise
correspondences before and after commonly-used refinement techniques.

<concept_desc>Computing methodologies Shape analysis</concept_desc>

[500]Computing methodologies Shape analysis

§ INTRODUCTION
Shape matching is a fundamental problem in geometry processing and computer graphics, with
applications ranging from shape interpolation [15] to shape exploration [16] and statistical shape analysis [6].

In this paper, we concentrate on the functional maps framework for computing correspondences between shapes, which has proven to be especially useful when dealing with non-rigid and especially near-isometric shape pairs. This approach, originally introduced in [28] and subsequently extended in multiple follow-up works, e.g. [17, 31, 16], is based on optimizing for a linear mapping between function spaces defined on the shapes, which can be conveniently encoded as a matrix. One attractive feature of this approach is that many desirable geometric objectives for the unknown pointwise correspondence, such as preservation of geodesic distances, can be conveniently encoded as constraints on the matrices representing functional maps, and often lead to simple convex optimization problems.

The majority of objectives when optimizing for a functional map between a pair of shapes consist of a) preservation of some pre-computed descriptors and b) a functional map regularization term, based on promoting some desired global map properties. This latter, regularization term is especially important for ensuring the overall global consistency of the computed map.

Common strategies for regularization include commutativity with the Laplacian operator, exploited in the original article [28], a mask promoting a slanted diagonal in the case of partial correspondences [31], or orthonormality of the functional map, which corresponds to local area preservation by the underlying pointwise map [28, 17, 33].

[trim=0cm 0cm 0cm 0cm,clip,width=1, grid=false]figures_new/teaser_v2.png

Comparison of different masks. Standard: the Laplacian commutativity operator can be equivalently formulated as a penalty or regularizing mask. Slanted: the weight mask proposed in [31] designed to promote a slanted structure. Ours: the mask proposed in this paper based on the resolvent of the Laplacian operator. The penalty masks of these three methods are visualized in the first row for an example pair from the SHREC dataset, and the corresponding optimized point-wise maps from the functional map pipeline are shown in the second row.

Unfortunately, although  functional map regularization has been instrumental in obtaining high-quality results on some shape correspondence benchmarks, these regularization terms often lack a theoretical foundation and indeed in some cases result in optimization objectives which would not be well-defined in the case of smooth surfaces.

In this paper, we introduce a novel way of regularizing functional maps, based on the concept of the resolvent of the Laplacian operator. Unlike the original approach, based on commutativity with the Laplacians, our method has a natural theoretically well-defined analogue in the case of smooth surfaces, and is guaranteed to be bounded even in the full (infinite dimensional) basis case. Moreover, our approach provides a simple one-parameter family of regularizers, that can be tuned depending on the expected approximate deviation from isometry in a given shape pair. Finally, we demonstrate on a wide range of challenging settings that our approach  leads to a quantitative and qualitative improvement for the computed functional and eventually pointwise maps, without incurring any additional computational complexity (see Fig. <ref> as an example).

§ RELATED WORK

Shape matching in its full generality is an extremely well-studied area of geometry processing and
computer graphics, and its full overview is beyond the scope of our paper. Therefore, below we only
review the most closely related methods based primarily on the functional maps framework. We refer
the interested readers to recent surveys including
[41, 40, 4] for an in-depth treatment of other
shape matching approaches.

Functional Maps Our approach fits within the functional map framework, which was
originally introduced in [28] for solving non-rigid shape matching
problems, and extended significantly in follow-up works, including
[17, 1, 18, 31, 10, 5] among many others
(see also [29] for an overview). The key observation in these techniques
is that it is often easier to estimate correspondences between real-valued functions, rather than
points on the shapes. This is because functions have a natural vector space structure, and moreover
functional transformations are most often linear, which means that functional maps can be encoded,
and thus optimized for, as small matrices in a reduced functional basis.

As observed by several works in this domain,
[17, 33, 31, 5] many natural properties on the
underlying pointwise correspondences can be expressed as objectives on functional maps. Most notably,
this includes: orthonormality of functional maps, which corresponds to the local area-preservation
nature of pointwise correspondences [28, 17, 33]; preservation
of inner products of gradients of functions, which corresponds to conformal maps [33, 5, 45];
preservation of pointwise products of functions, which corresponds to functional maps arising
from point-to-point correspondences [26, 25]; slanted diagonal structure of
functional maps, which corresponds to correspondences between partial shapes
[31, 21]. Similarly, several other regularizers have been proposed, including using
robust norms and matrix completion techniques [17, 18], exploiting the relation between functional maps in different directions [11], the map adjoint  [14], and powerful
cycle-consistency constraints [16] in the context of shape collections, among many
others. More recently constraints on functional maps have been introduced to promote
continuity of the recovered pointwise correspondence [30], maps between curves defined on
shapes [12], kernel-based techniques aimed at extracting more information
from given descriptor constraints [43], and an approach for incorporating
orientation information into the functional map infererence pipeline [34] among

Among all of these, perhaps the most widely-used building block for regularizing functional maps,
introduced in [28] and extended in follow-up works, including
[44, 31, 20, 21], is based on the commutativity with the Laplacian operators, which implies
a diagonal (or slanted diagonal in the case of partial correspondence) structure for functional
maps. To promote this structure, the most common method (see also Chapter 2.4 in
[29]) consists in adding an energy to the functional map estimation
pipeline, which penalizes the failure of the unknown functional map to commute with the
Laplace-Beltrami operators on the source and target shapes. Conveniently, this term still leads to a
convex optimization problem during functional map estimation. Moreover, as observed in
[28], for functional maps arising from pointwise ones, this term is zero
if and only if the map preserves geodesic distances exactly.

Unfortunately, although functional map regularization via commutativity with the Laplace-Beltrami
operators has been instrumental in obtaining high quality results in challenging benchmarks, the
exact properties of this regularization are not well-understood. In particular, as we show below,
the commonly-used energy is not bounded in the full (infinite-dimensional) basis in general. Instead, our novel regularizer, although based on a
similar underlying principle, overcomes this limitation, and both leads to a theoretically better
justified energy, and a practical improvement on a range of challenging datasets.

Optimal Transport We also note briefly that other commonly-used relaxations for
matching problems include those based on optimal transport,
e.g. [39, 22, 42].  These techniques have
recently gained prominence especially due to the computational advances for adressing large-scale
transport problems, primarily using the Sinkhorn method with entropic regularization
[8, 38]. Furthermore, other techniques that exploit the
formalism of optimal transport, for solving assignment problems include the recent variants
of Product Manifold Filter using Gaussian and Heat Kernels
[22, 42]. Interestingly, these latter methods argue that
near-isometric shape matching can be better addressed via preservation of general functional
kernels rather than preservation of e.g. geodesic distances. Our modification of the
regularization term of functional maps can also be seen through a similar light, as we show that a
more theoretically justified functional operator leads to an improvement of the overall structural
properties of the map, which eventually result in more accurate functional and pointwise maps.

§ BACKGROUND

Our work is based on the functional map representation and the estimation pipeline. Below we review
the basic notions and the main steps for estimating a map between a pair of shapes using this
framework, and refer the interested reader to a recent set of course notes
[29] for a more in-depth discussion.

Basic Pipeline Given a pair of shapes, $S_1,S_2$ represented as triangle meshes, and
containing, respectively, $n_1$ and $n_2$ vertices, the general pipeline for computing a map between
them using the functional map representation, consists of the following main steps:

* Compute a small set of $k_1 << n_1$ and $k_2 << n_2$ basis functions on each shape. The most common choice consists in using the first $k$ eigenfunctions of the
Laplace-Beltrami operator of each shape, although other bases derived from  the Hamiltonian operator [7] and more localized basis functions [27, 24] have also been used.

* Compute a set of descriptor functions on each shape, that are expected
to be approximately preserved by the unknown map. Store their coefficients in the corresponding
bases as columns of matrices $A_1, A_2$.

* Compute the optimal functional map $C$ by solving the following optimization problem:

\begin{align}
C_{\text{opt}} = \argmin_{\matr{C}_{12}} E_{\text{desc}}\big(\matr{C}_{12}\big) + \alpha E_{\text{reg}}\big(\matr{C}_{12}\big)
\end{align}

where the first term aims at the descriptor preservation: $E_{\text{desc}}\big(\matr{C}_{12}\big)
= \big\Vert \matr{C}_{12} \matr{A}_1 - \matr{A}_2\big\Vert^2$, whereas the second term regularizes
the map by promoting the correctness of its overall structural properties. As mentioned above, the
most common approach consists of penalizing the failure of the unknown
functional map to commute with the Laplace-Beltrami operators, which can be written as:

\begin{align}
\label{eq:energy:laplacian_comm}
E_{\text{reg}}(C_{12}) = E_{\text{comm}}\big(\matr{C}_{12}\big) = \big\Vert \matr{C}_{12}\Delta_1 - \Delta_2\matr{C}_{12} \big\Vert^2
\end{align}

where $\Delta_1$ and $\Delta_2$ are the Laplace-Beltrami operators of the two shapes expressed in the respective bases. Here, and throughout the rest of our paper, unless stated otherwise $\Vert \,\cdot \,\Vert$ denotes the matrix Frobenius norm.

* Convert the functional map $C$ to a point-to-point map, for example using an iterative approach,
such as the Iterative Closest Point (ICP) in the spectral embedding, or using other more advanced
techniques [32, 10].

One of the attractive properties of this pipeline is that the functional map computation in step 3 leads to a simple (convex) least squares optimization with a relatively small number of unknowns, independent
of the number of points on the shapes. This step has been further extended e.g. using more powerful descriptor preservation constraints via commutativity [26], or using manifold optimization [19] among many others (see also Chapter 3 in [29]).

§ FUNCTIONAL MAP REGULARIZATION

[trim=0cm 0cm 0cm 0cm,clip,width=0.95, grid=false]figures_new/eg_methodology.png
(2,60)Given descriptor
(2,30)Ground-truth map

(45, 68.5)Mask
(64, 70)Functional
(68, 67)map
(88, 70)Pointwise
(92, 67) map

Given a single pair of WKS descriptors, we optimize a 100-by-100 functional map using the
standard mask and our resolvent mask and compare to ground truth. The converted point-wise maps are shown on the right. We can see that the functional map stemming from our resolvent mask has less noise than the functional map computed with the standard mask. Also, our resolvent mask leads to a
point-wise map with better quality.

Our main goals are to analyze the commonly-used functional map regularization term, to bring attention to some of its
theoretical limitations, and to propose a novel regularizer with better theoretical properties, which lead to practical
improvements. Therefore, in this in this section, we first consider the standard Laplacian-commutativity term, and then
introduce our new regularizer based on the resolvent operator of the Laplacian.

§.§ Reformulation of the Laplacian-Commutativity term

As mentioned above, the Laplacian-commutativity term given in Eq. (<ref>) was first introduced
to promote approximately isometric correspondences. If the functional map is expressed in the basis given by the
eigenfunctions of the Laplace-Beltrami operator, and letting $\Lambda_1$ and $\Lambda_2$ represent vectors that store
the eigenvalues of the Laplacians of shape $S_1$ and $S_2$ respectively, this term can be equivalently reformulated as:

\begin{equation} \label{eq:comm_as_mask}
\begin{split}
&= \big\Vert \matr{C}_{12}\Delta_1 - \Delta_2\matr{C}_{12} \big\Vert^2 \\
&= \big\Vert\matr{C}_{12}\text{diag}\big(\Lambda_1\big) - \text{diag}\big(\Lambda_2\big)\matr{C}_{12}\big\Vert^2\\
&= \big\Vert \matr{C}_{12} \odot \big(\mathbf{1}_{k_2}\Lambda_1\tran\big) - \big(\Lambda_2\mathbf{1}_{k_1}\tran\big)\odot \matr{C}_{12} \big\Vert^2\\
&= \big\Vert \big(\mathbf{1}_{k_2}\Lambda_1\tran - \Lambda_2\mathbf{1}_{k_1}\tran \big) \odot \matr{C}_{12}\big\Vert^2\\
%&\triangleq E_{\text{mask}}\big(\matr{C}_{12}\big) =  \big\Vert\matr{M}_{LB} \odot \matr{C}_{12} \big\Vert^2 \\ %Old , wrong, version
&\triangleq E_{\text{mask}}\big(\matr{C}_{12}\big) = \sum_{ij} \left[\matr{M}_{LB} \right]_{ij} \left[ \matr{C}_{12} \right]_{ij}^2~,
\end{split}\end{equation}

where $\odot$ is the entry-wise matrix multiplication, $\mathbf{1}_{k}$ is a $k$-dim all-ones vector and $[A]_{ij}$ denotes the $(i,j)$ entry of a matrix $A$.

In other words, the regularization term  $E_\text{comm}$ can be written as a product between the matrix $M_{LB}$, which
we call the penalty mask matrix and the squares of the entries of the unknown functional map $C_{12}$. In the case of $E_\text{comm}$, the
entries of $M_{LB}$ are given as $M_{LB}(i,j) = (\Lambda_2(i)-\Lambda_1(j))^2$. Fig. <ref> shows a heat
map of this matrix (see the first row in the “Mask” column).

Unfortunately, this simple formulation has certain fundamental undesirable properties. In particular, in the case of smooth surfaces, Laplace-Beltrami operators are unbounded operators on square-integrable functions [23].
Consequently, in general, the energy term $\| C_{12}\Delta_{1} - \Delta_{2}C_{12} \|^{2}$ is undefined on smooth surfaces.
As a simple example, consider the situation where $C_{12} = Id$, the identity operator, and the two surfaces are scaled versions of each other (i.e., $\Delta_2 = c\Delta_1$ for some constant $c \neq 1$). In this case, $\| C_{12}\Delta_{1} - \Delta_{2}C_{12} \|^{2} = |1-c|^2 \| \Delta_1 \|^2$ is infinite.

That being said, the ill-definiteness of $\| C_{12}\Delta_{1} - \Delta_{2}C_{12} \|^{2}$ is not a mere question of scale. Recall that, by Weyl's law [9], large Laplacian eigenvalues of surfaces can be estimated in terms of the surface area $S$ as follows:

\begin{equation}
\lambda_k \sim \frac{4 \pi}{S}k~.
\end{equation}

Thus, by rescaling the surfaces such that their areas match guarantees that their eigenvalues have comparable
asymptotic growth. This, however, is not sufficient to make $\| C_{12}\Delta_{1} - \Delta_{2}C_{12} \|^{2}$ well-defined. This can be seen from another simple, if slightly artificial, example. Consider a sphere ($\mathcal{M}_1$) and a flat square torus ($\mathcal{M}_2$), both of unit
area. While these surfaces are not diffeomorphic, they are among the few whose Laplace-Beltrami eigenvalues can be explicitly computed (see [36, 37] among many others), which is why we consider them here. Said spectra, as well as the corresponding Weyl estimate are illustrated in Fig. <ref>.
For simplicity's sake, we once again pick a functional map of the form $C_{12} = Id$.
We illustrate the spectra, as well as the corresponding $\| C_{12}\Delta_{1} - \Delta_{2}C_{12} \|^{2}$ on Fig.  <ref>. Notice the divergence of the Laplacian commutator energy. On the same figure, we illustrate its proposed replacement, defined in the next section. Note its rapid convergence.

The spectra of a sphere and a square torus of unit area, as well as the corresponding Weyl estimate.

In addition to being ill-defined in the continuous setting, the Laplacian commutativity mask has another significant
problem. Namely, it penalizes the high frequency entries of functional maps in the same way as the low frequency ones,
despite the greater instability of the former, and in this way fails to exhibit the funnel-like structure observed in
ground-truth functional maps (see e.g. Fig. <ref> above or Figure 4 in
[28]). Recall that a mask serves as a penalty during the optimization of the functional map
matrix.  Thus, the mask and the ground truth functional map should be complementary, in the sense that for the regions
of ground-truth functional maps with smaller (resp. large) values, the mask should add more (resp. less) penalty.

We illustrate this phenomenon in Fig. <ref>, which shows the average of the squared values of functional maps constructed from the ground-truth correspondences using 250 eigenfunctions over 100 FAUST non-isometric shape pairs. We then compare its structure to the Laplacian commutativity mask (labeled “Standard”), averaged over the same shape pairs. We perform the same computation for the heuristic slanted mask introduced in [31] and, finally, for our proposed resolvent-based mask defined in the next section. Notice that our mask better reproduces the funnel-like structure of the ground-truth maps.

The Laplacian and Resolvent commutator energies for $C_{12} = Id$ computed using the lowest $k$ eigenvalues
of the unit area sphere and torus (see Fig. <ref>). Note the rapid convergence of $\|
C_{12}R(\Delta_{1}) - R(\Delta_{2})C_{12} \|^{2}$ and the divergence of $\| C_{12}\Delta_{1} - \Delta_{2}C_{12}
\|^{2}$. The spectra are rescaled with respect to the largest eigenvalue with $k = 100$, as described in
Eq. (<ref>).

§.§ A novel regularization based on the Resolvent

In this section, we propose an alternative form for the Laplacian commutativity regularizer that overcomes the problems identified in the previous section. The resulting regularizer is better from both the theoretical and practical standpoints.

The first issue of the original Laplacian-commutativity term arises from the fact that the Laplace-Beltrami operator is
unbounded. We thus propose to replace it with a meaningful bounded operator, namely its
resolvent. Below, we give a brief overview and refer the reader to [35] for the detailed
functional-analytic underpinnings of the following discussion.

Resolvent Let $A: \mathcal{D} \to \mathcal{H}$ be a densely-defined closed  operator on some Hilbert space $\mathcal{H}$ with domain $\mathcal{D} \subset \mathcal{H}$. Let $\rho(A) \subset \mathbb{C}$ be the set of all complex numbers $\mu$ such that $R_{\mu}(A) = \left( A - \mu Id \right)^{-1}$ is defined and bounded. The set $\rho(A)$ and the operator $R_{\mu}(A)$ are known as the resolvent set of $A$ and the resolvent (operator) of $A$ at $\mu$, respectively. The resolvent set $\rho(A)$ is the complement of the spectrum of $A$ in the complex plane.

The general idea of considering the resolvent of an unbounded operator rather than the operator itself comes from the
notion of norm-resolvent convergence, which is used to study the convergence of unbounded self-adjoint
operators [35]. In that sense, our choice to use the resolvent of the Laplace-Beltrami operator is a natural one. Note that here, closedness is a technical condition used in the definition of the resolvent. In particular, it is satisfied by self-adjoint operators, such as the ones we consider.

Before proceeding further, we slightly generalize our problem. In what follows we will consider powers of the Laplacian
rather than the Laplacian itself. Specifically, we will use $\Delta^{\gamma}$ for some $\gamma > 0$ rather than
$\Delta$. As explained below, the parameter $\gamma$ controls the funnel-like structure of the resulting mask. Thus, this parameter takes on an important role in the numerical tests reported later in this paper.

Now since the Laplace-Beltrami operator $\Delta$ is positive and self-adjoint, its spectrum is contained in the non-negative real line. The same also holds for $\Delta^\gamma$. Thus, we are guaranteed that $R_{\mu}(\Delta)$ is a well-defined bounded operator for any complex $\mu$ not in the non-negative real line.

For our purposes, the resolvent of $\Delta^{\gamma}$ is expressed as

\begin{equation}
R(\Delta^{\gamma}) = \left( \Delta^{\gamma} - (a + i b) Id \right)^{-1}~,
\end{equation}

where $i$ is the imaginary unit and $a,b \in \mathbb{R}$. Our proposal is to use the resolvent of the Laplacian instead of the Laplacian in the commutator term. We thus define a new energy term,

\begin{equation}
\label{eq:res_penalty}
E_{\text{resolvent}}\big(\matr{C}_{12}\big) = \big\Vert \matr{C}_{12} R \left( \Delta_1^{\gamma} \right) - R \left( \Delta_2^{\gamma} \right) \matr{C}_{12} \big\Vert^2~.
\end{equation}

Before proceeding further, recall that we only need to consider bounded functional maps $C_{12}$, since functional maps that arise as pullbacks of diffeomorphisms are bounded. This last fact is shown in appendix <ref>, for completeness.
Our proposal is motivated by the following result.

[trim=0cm 0cm 0cm 0cm,clip,width=1, grid=false]figures_new/eg_gt_and_diff_masks.png
(28,27)Standard mask
(53,27)Slanted mask

(81,27)Our mask

The average squared ground-truth functional map and three penalty masks over 100 FAUST non-isometric shape pairs.

Let $C_{12}$ be a bounded functional map. Then, in the operator norm,

\begin{equation}
\left\| \matr{C}_{12} R \left( \Delta_1^{\gamma} \right) - R \left( \Delta_2^{\gamma} \right) \matr{C}_{12}  \right\|^{2}_{op} < \infty~.
\end{equation}

The result directly follows from the fact that products and linear combinations of bounded operators are bounded. As defined above, $R(\Delta^{\gamma})$ is bounded, which concludes the proof.

As explained above, there is no analogous guarantee for the commonly used Laplacian commutativity, since the Laplace-Beltrami operator itself is unbounded. In contrast, our resolvent-based commutativity is always well-defined, and leads to a commutator with bounded operator norm, even in the case of smooth surfaces.

Notice that Theorem <ref> holds in the operator norm, rather than the more convenient Frobenius norm, which we use to define the energy in Equation (<ref>). Our usage of the Frobenius norm can be justified in two ways. First, a Frobenius norm version of Theorem <ref> holds for $\gamma >1/2$. This is shown in Lemma <ref> of Appendix <ref>. Thus, for $\gamma > 1/2$, the usage of the Frobenius norm is perfectly justified in the smooth setting. For lower values of $\gamma$, we invoke the fact that, in practice, we work on a finite dimensional vector space. Since all norms are equivalent on finite dimensional vector spaces, we are free to replace the operator norm with the Frobenius one. In that case, we lose the guarantee that the energy $E_{resolvent}$ makes sense in the full Laplace-Beltrami basis in the smooth setting. Yet, as shown in the experiments reported below, this does not seem to cause any issue.

Similarly to the usual Laplacian commutativity term $E_{\text{comm}}$, $E_{\text{resolvent}}$ can also be expressed as a mask matrix. This can be done as follows. Notice that $R(\Delta^{\gamma})$ is diagonal in the Laplacian eigenbasis. Moreover, if $\Delta$ has eigenvalues $\{ \lambda_{n} \}_{n=0}^{\infty}$, then, by matrix inversion, the eigenvalues $\{ r_n \}_{n=0}^{\infty}$ of $R(\Delta^{\gamma})$ will be given by:

\begin{equation} \label{eq:mtd:resolvent_eigenvalues}
r_n = \frac{1}{\lambda^{\gamma}_n - a - ib}  = \frac{\lambda^{\gamma}_n - a}{(\lambda^{\gamma}_n - a)^2 + b^2} + \frac{b}{(\lambda^{\gamma}_n - a)^2 + b^2}i~.
\end{equation}

In practice, we choose $a = 0$ and $b = 1$, which yields

\begin{equation}
r_n = \frac{\lambda^{\gamma}_n}{(\lambda^{\gamma}_n )^2 + 1} + \frac{1}{(\lambda^{\gamma}_n )^2 + 1}i~.
\end{equation}

Since the square of the Frobenius norm independently considers the real and imaginary parts of a matrix, we can re-express Equation (<ref>) in terms of two new matrices $M_{Re}$ and $M_{Im}$, defined below.

\begin{equation}
\begin{aligned}
E_{\text{resolvent}}\big(\matr{C}_{12}\big) &= \big\Vert \matr{M_{Re}} \odot \matr{C}_{12}\big\Vert^2 + \big\Vert \matr{M_{Im}} \odot \matr{C}_{12}\big\Vert^2\\
&= \sum_{ij} \left[M_{Re} \right]_{ij}^2 \left[ C_{12} \right]_{ij}^2 + \sum_{ij} \left[M_{Im} \right]_{ij}^2 \left[ C_{12} \right]_{ij}^2\\
&= \sum_{ij} \left( \left[\matr{M}_{Re}\right]_{ij}^2 + \left[M_{Im} \right]_{ij}^2 \right) \left[ C_{12} \right]_{ij}^2
\end{aligned}
\label{eq:mtd:resolvent2masks}
\end{equation}

The matrices $\matr{M}_{Re}$ and $\matr{M}_{Im}$ correspond to the real and imaginary parts of the eigenvalues of the resolvent, respectively. Explicitly, these matrices are given by

\begin{align}
\label{eq:mtd:def:mask_real} &\matr{M}_{Re}(i,j) = \frac{\Lambda_2(i)^{\gamma}}{\Lambda_2(i)^{2\gamma} + 1} - \frac{\Lambda_1(j)^{\gamma}}{\Lambda_1(j)^{2\gamma} + 1}\\
\label{eq:mtd:def:mask_imag} &\matr{M}_{Im}(i,j) =  \frac{1}{\Lambda_2(i)^{2\gamma} + 1} - \frac{1}{\Lambda_1(j)^{2\gamma} + 1}
\end{align}

These matrices are not quite masks in the sense of Eq. (<ref>). As per Equation (<ref>), the above two matrices can be combined into a single mask

\begin{equation}
\label{eq:mtd:res_mask} \matr{M}_{res}(i,j) = \matr{M}_{Re}(i,j)^2 + \matr{M}_{Im}(i,j)^2~.
\end{equation}

The split of $\matr{M}_{res}$ into $\matr{M}_{Re}$ and $\matr{M}_{Im}$ will be revisited later, when we explore beyond the established theory and consider a mask constructed from weighted combinations of these two matrices.

§.§ Rescaling the Spectra

As it stands now, we observe that in practice the mask defined in equation (<ref>) decays too quickly as the Laplacian eigenvalues grow.
In other words, the mask is not sufficiently sensitive to the higher frequencies. This is due to the scale introduced by the parameter $b$ found in the definition of the resolvent. There are a few equivalent ways of addressing this issue.

Our approach is as follows. We begin by computing the $k$ lowest eigenvalues of the considered Laplacians. Then, both spectra are rescaled according to the rule

\begin{equation} \label{eq:mtd:rescaling}
\Lambda_{i} ~\longmapsto~ \frac{\Lambda_{i}}{ \max \big( \max(\Lambda_1), \max(\Lambda_2) \big)}~.
\end{equation}

In other words, the spectra of the source and target shape are rescaled by the same factor, such that the largest considered eigenvalue (over both spectra) becomes equal to $1$. From Equation (<ref>), one can see that this rescaling can be absorbed into a choice of $b$ and a change in the weight of the resolvent energy in the overall energy. Thus, conceptually, the spectral rescaling is equivalent to a choice of the parameter $b$ used in the definition of the resolvent. Alternatively, recall that rescaling the spectrum is equivalent to rescaling the surface.
Consequently, this rescaling does not affect the theoretical guarantees of the previous section.

Returning to the previously discussed example of the sphere and torus of unit area (see Fig. <ref>), we illustrate the resolvent energy computed by the above procedure in Fig. <ref>. Note the rapid convergence of the resolvent commutator energy (in red) and the divergence of the Laplacian commutator energy (in blue).

§.§ Mask Structure as a Function of $\gamma$

The resolvent mask introduced in Eq. (<ref>) with different $\gamma$. The darker the region, the larger the penalty. Here, $\Delta_1$ and $\Delta_2$ are near-isospectral, which explains the approximately symmetric structure of the mask.

Having defined the way to compute the resolvent mask, we are now ready to explore the way in which tuning the parameter $\gamma$ controls its the funnel-like structure.

Fig. <ref> illustrates the resolvent mask for different values of $\gamma$. The darker the region, the more penalized the corresponding entry in the functional map $C_{12}$ would be. Notice that the funnel-like structure of the mask changes with different values of $\gamma$.

From Fig. <ref>, we see that the behaviour of the mask as a function of $\gamma$ can be separated into two regimes.
The first  corresponds to $\gamma \in (0,1]$. There, increasing $\gamma$ results in a narrowing of the funnel-like shape of the mask. Thus, the larger $\gamma$ is, the more the mask penalizes functional maps that take eigenfunctions of $\Delta_{1}$ to eigenfunctions of $\Delta_{2}$ with distant eigenvalues. Correspondingly, small values of $\gamma$ are more lax in that regard, allowing for maps between eigenfunctions with quite different eigenvalues.
The former choice of $\gamma$ is appropriate when the shapes under consideration are approximately isometric, as one then expects the eigenfunctions and eigenvalues of both surfaces to be roughly the same. The latter choice is sound for shape pairs that are further away from isometry.

The second regime corresponds to $\gamma > 1$. There, we observe an inversion of the funnel-like structure, with the reverse funnel shape being more and more pronounced as $\gamma$ increases. This results in a low penalty for maps between low frequency eigenfunctions, which does not respect the shape of the ground-truth maps (see Fig. <ref>).

Later in this paper, we report empirically obtained optimal values for $\gamma$ on a benchmark dataset (see Fig. <ref>).

§.§ Proposed Functional Map Energy

In summary, we propose the following energy to compute a functional map for a pair of shapes, with a set of given descriptors:

\begin{equation}
\label{eq:energy:complete}
E\big(\matr{C}_{12}\big) = \alpha_1E_{\text{desc}} + \alpha_2 E_{\text{mult}} + \alpha_3 E_{\text{orient}} + \alpha_4 E_{\text{resolvent}}
\end{equation}

where $E_{\text{desc}}$ is the descriptor-preserving term defined in Sec. <ref>, $E_{\text{mult}}$ is the descriptor-commutativity term defined as $E_{\text{mult}} = \sum_{i} \big\Vert \matr{C}_{12}\matr{D}_{1i}^{\text{mult}} - \matr{D}_{2i}^{\text{mult}}\matr{C}_{12}\big\Vert^2$ introduced in [26], $E_{\text{orient}}$ is the orientation-preserving term introduced in [34]. $E_{\text{resolvent}}$ is the our new term introduced in Eq. (<ref>) and discussed above. As mentioned above, when functional maps are expressed in the Laplace-Beltrami eigenbasis, this term can be  written via a penalty using the mask matrix given in Eq. (<ref>) .

Before proceeding to a more extensive evaluation of the proposed energy, in Fig. <ref> above, we provide an example of a functional map obtained using this energy and compare it to one obtained using the standard Laplacian commutator regularizer. Fig. <ref> also shows the resolvent mask with $\gamma = 0.5$ (second row), compared to the standard mask (first row).
The rightmost column shows the quality of pointwise maps, recovered from the functional maps (shown in the second rightmost column) using both our and the standard Laplacian regularizers.

Changing the functional map size. We randomly select 50 FAUST non-isometric pairs and 50 TOSCA isometric pairs. For each of the pairs, we only use one pair of WKS [2] descriptors. We then optimize for a functional map with different size ranging from 20 to 250 using different Laplacian mask terms. The solid lines represent the initialization with different masks and the dashed lines are the ICP-refined results. We report the per-vertex and the direct error measure. We can see that the proposed mask is much more stable than the standard and the slanted mask as the functional map size increases.

§ RESULTS

We tested the proposed approach using a MATLAB-based implementation, which we adapt from the state-of-the-art functional map approach in [34]. Here, we first describe the benchmark datasets and the baseline methods.

We use the two datasets introduced in [34]. These consist of shapes from the FAUST [6] and TOSCA [3] datasets, which were remeshed so that the shapes have different triangulations, and are no longer in one-to-one correspondence, making the matching more challenging and realistic. Specifically, we include 300 FAUST shape pairs and 284 TOSCA pairs for evaluation.

To evaluate our new regularizer, we compare to the following two masks:

* "standard": the standard Laplacian-commutativity mask, which is defined as $\matr{M}_{LB}(i,j) = \big(\Delta_2(i) - \Delta_1(j)\big)$ as discussed before.

* "slanted": the heuristic slanted diagonal penalty mask proposed in [31], which is defined as

$$\matr{M}(i, j) = \exp{\left(-\eta \sqrt{i^2 + j^2}\right)}\left\Vert \frac{\vec{n}}{\Vert \vec{n} \Vert} \times \left((i,j)\tran - \vec{p}\right)\right\Vert$$

where $\vec{p} = (1,1)\tran$, and $\vec{n} = (1,r/k)\tran$ is the line direction with slope $r/k$, where $r$ is the estimated rank of the functional map, and $k$ is the size of the square functional map. This weight matrix is originally applied to partial shape matching in [31], and we use it as a mask matrix to regularize the functional map. In our tests, $\eta$ is set to the default value 0.03 as suggested in the original paper. This mask is illustrated on Fig. <ref>. Note that it exhibits the desired funnel-like structure for the lower part of the spectrum, but not the upper part.

We also compare three different settings:

* Initialization. The wave kernel signatures (WKS) [2] are used to construct $E_{\text{desc}}, E_{\text{mult}}$, and $E_{\text{orient}}$ in Eq. (<ref>). Then we optimize the functional map w.r.t. the energy defined in Eq. (<ref>) with three different masks, namely, the standard, slanted, and our resolvent mask.

* ICP refinement. After initialization, we use ICP [28] to refine the computed functional maps.

* BCICP refinement. After initialization, we use the recently proposed BCICP algorithm [34] to refine the computed functional and pointwise maps, using the open-source implementation and default parameters provided by the authors.

In our experiments, we measured the quality of the functional maps and the recovered point-wise maps:

* Point-wise maps. Since most shapes contain left-right symmetries, which are indistinguishable for intrinsic methods, in each dataset, we considered both the ground-truth direct and symmetric correspondences. To measure the accuracy of a computed map, we used the following measures:

* per-vertex error: for each vertex we accept the ground-truth direct and symmetric correspondences and take the minimum as the error of this vertex. This measure reflects the accuracy of the map regardless of the symmetry.

* direct error: we compute the average per-vertex error to the direct ground-truth correspondences only. This measure reflects both the accuracy and the smoothness of the map.

* Functional maps. We can evaluate a functional map by the quality of its recovered point-wise map, or by measuring the penalty from a given mask. Specifically, for a given functional map $\matr{C}$ and a  mask matrix $\matr{M}$, we can measure the total penalty as
$\sum_{ij} \left[\matr{M} \right]_{ij} \left[ \matr{C} \right]_{ij}^2$ .

[trim=0cm 0cm 0cm 0cm,clip,width=1, grid=false]figures_new/eg_diff_fMapSize_v2.png
(5,70)Given descriptor
(3.5,35)Ground-truth map

Given one pair of WKS descriptors, as visualized on the left top, we use different Laplacian mask terms to optimize for a functional map with size $k$. The ground-truth map is visualized on the bottom left. We can see that our mask is much more stable over different size $k$.

For each dataset, we compare the quality of the functional maps between the standard mask (blue curves), the slanted mask (yellow curves), and our resolvent mask (red curves) with the same set of parameters, where a single pair of descriptors is used to optimize a 100-by-100 functional map. The comparison is made in three different settings: comparing the initialization directly (dotted lines), comparing the initialization with ICP refinement (dashed lines), and with the BCICP refinement (solid lines).
Specifically, these curves are measured on 300 FAUST shape pairs and 284 TOSCA pairs. The average direct errors are reported in the legends.

§.§ Effect of functional map size

For the functional maps pipeline, a set of corresponding descriptors is given as input. We then optimize
a $k \times k$ functional map by minimizing an objective function based on Eq.  (<ref>).
Therefore, we have to solve for $k^2$ variables. If $k$ is a smaller value, e.g., $k < 50$, the optimization problem
is easier to solve since the number of variables is small. However, if $k$ is too small, the information that is encoded
into the optimized functional map is limited to the low-resolution of the spectrum of the shapes, and hence it will be hard
to transfer detailed information. On the other hand, if $k$ is too large, solving the optimization problem is potentially more time-consuming. Even more importantly,
this optimization problem can become under-constrained when the number of variables exceeds the constraints stemming from the input descriptors.
In this case, we need effective regularizers to regularize the functional maps. In real experiments, the choice of the parameter $k$ is a key hyper parameter.

To quantify the stability and the effectiveness of the proposed resolvent mask, we conduct the following test: for each of the test pairs, we only use one pair of corresponding WKS descriptors. We then fix this descriptor pair and optimize for functional maps with different sizes ranging from 20 to 250. We randomly select 50 pairs of FAUST and 50 pairs of TOSCA, and report the average error over the tested shape pairs w.r.t. different functional map size in Fig. <ref>. We can see that the standard mask fails to regularize the functional map with a large size: the average per-vertex error is three or four times larger than ours. At the same time, the slanted mask has a better performance than the standard one in the per-vertex measure. However, the slanted mask has large direct errors, which suggests that the smoothness is not well preserved.
We believe this is due to the fact that the orientation-preserving regularizer starts to fail to disambiguate the symmetry as the functional map size $k$ increases. In this case, the slanted mask cannot help the orientation-preserving operator, while our mask can strengthen the functionality of the orientation-preserving operator and leads to maps with much lower per-vertex and direct error.
In summary, even with limited constraints from a single pair of WKS descriptors, increasing the number of variables does not significantly affect the performance of our mask. Fig. <ref> shows an illustrative example.

This test shows that our mask is much more stable and can better regularize larger functional maps even in very challenging cases with little input information or constraints. Also, it suggests that with this new mask, we no longer need to tune the parameter $k$ as much as needed by the standard mask to achieve a better result.

§.§ Evaluation on shape matching

Changing the number of input descriptors. We use different numbers of descriptors to optimize a 100-by-100 functional map with different masks. The results are reported for 300 FAUST shape pairs. As the number of descriptors increases, the results improve for all the masks.

[trim=0cm 0cm 0cm 0cm,clip,width=1, grid=false]figures_new/res_face_expression_withICP.png

(96,17)-90+ ICP

Example. Comparing the quality of the initial maps and after ICP refinement via texture transfer. (One pair of descriptors is used to optimize a 50-by-50 functional map)

[trim=0cm 0cm 0cm 0cm,clip,width=1, grid=false]figures_new/res_icp.png

Example. Comparing the quality of the ICP refined maps via color transfer. (One pair of descriptors is used to optimize a 100-by-100 functional map)

In this experiment, we compare our resolvent mask to the standard and the slanted mask on a larger FAUST and TOSCA dataset in three settings: the initialization, with ICP refinement, and with BCICP refinement.
To make a fair comparison between different masks, the weights for different terms (the $\alpha_i$ in Eq. (<ref>)) are fixed across different test pairs and different test masks.

In this test, for each test pair, we use one pair of WKS descriptors to optimize a 100-by-100 functional map.
As reported in Fig. <ref>, our mask leads to 40.4%, 43.1%, and 16.3% improvement w.r.t the settings of the initialization, with ICP refinement, and with BCICP refinement respectively, over the best of the standard and the slanted mask on the FAUST dataset regarding the average direct error. Similarly, for the TOSCA dataset, ours achieves 34.2%, 42.3%, and 24.8% improvement respectively over the best of the standard and the slanted mask.

Fig. <ref> shows a test on the FAUST dataset, where we use different number of carefully curated input
descriptors based on WKS [2] using parameters from [34] to compute 100-by-100
functional maps. As the number of input descriptors increases, more constraints are added to regularize the functional
map. Therefore, for all the masks, and all the test settings, the results improve.  Observe that when the number of
variables is small (as shown in Fig. <ref>  for small $k$) or the number of
descriptor constraints is large (as shown in Fig. <ref> for large descriptor number), all
masks perform well since the problem is well-constrained. However, when this is not the case, our resolvent mask can still
regularize the functional map better than the other two.  For completeness,  in
Appendix <ref> we also include results with the BCICP refinement. We remark that this refinement is very computationally and memory intensive due, in part, to requiring all-pairs geodesic distances, but can, as such, improve upon even very poor quality maps.

[trim=0cm 0cm 0cm 0cm,clip,width=1, grid=false]figures_new/res_bcicp_lion_cat.png

Example. Here we show a challenging pair of a lion and a cat and the results are refined by BCICP. We can see that, when the initial maps are in a low quality and the BCICP refinement fails to improve the initial maps, our new mask still gives a more reasonable map. (Twenty pairs of descriptors are used to optimize a 50-by-50 functional map)

Fig. <ref> shows a qualitative example of using a single descriptor pair to optimize for a 100$\times$100 functional map. The first row shows the quality of the initial maps with different masks, and the second rows the quality of the corresponding maps refined by ICP. We can see that our resolvent mask outperforms the other two masks. Also the quality of our initial map is close to the ICP refined map, which shows that with our mask, we do not rely on the post-processing refinement as much as the other two. Fig. <ref> shows another example of the results refined by ICP. Fig. <ref> shows a challenging non-isometric example with results refined after BCICP.

§.§ Analysis of the complex resolvent mask

As shown in Fig. <ref>, the funnel-shape of our resolvent mask aligns well with the ground-truth functional map, which leads to a better performance over the standard and the slanted mask. To further analyze the properties of our mask, we also conduct the following experiments: We test the range of the parameter $\gamma$, and the relative weight between the complex and the real part to construct the resolvent mask. Moreover, we investigate the correlation between the mask penalty added on a functional map and the corresponding recovered point-wise map.

§.§.§ Different parameters for the complex resolvent mask

Left: changing the parameter $\gamma$ (as defined in Eq. (<ref>)); Right: changing the relative weight $w$ between the imaginary part and the real part (see the original definition in Eq. (<ref>)). The results are based on 50 TOSCA pairs.

In this section, we empirically tune the parameters in our resolvent based mask. First, we explore different values of $\gamma$. Recall that $\gamma$ controls the funnel-like structure of the mask. Thus, it is expected that tuning $\gamma$ can influence the functional map quality.

On Fig. <ref>, we report results for 100 FAUST pairs, where the result of the standard mask is colored blue, and ours is colored red. As is shown on the left, when $\gamma$ lies between 0 and 1, our mask always outperforms the standard mask over both the initialization (solid lines) and after the ICP refinement (the dashed lines).

Note that the case $\gamma = 1$ corresponds to the resolvent of the Laplacian. Thus, the fact that our mask outperforms the standard one for $\gamma = 1$ experimentally justifies the usage of the resolvent of the Laplacian, rather than the Laplacian itself. Note also that, as suggested in Sec. <ref>, the resolvent mask performs poorly for $\gamma > 1$.

Finally, we explore the relative contribution of the real and imaginary parts of the resolvent mask in order to analyze the utility of the two components of our mask construction.  For this, we analyze the accuracy of the computed pointwise maps when allowing $M_{Re}$ and $M_{Im}$ to have different weights, $(1-w)$ and $w$ for $w \in [0,1]$, respectively. As shown in Fig. <ref> (right), the convex shape of the red curve (where the weight $w$ changes) demonstrates that both the real part and the imaginary part contribute to the improvement over the standard mask. Note that our mask with any convex combination between the real part and the imaginary parts outperforms the standard mask. In practice we always use the equal weight $w=0.5$.

§.§.§ Correlation between the mask penalty and map accuracy

Correlation between the mask penalty on a functional map and the accuracy of the corresponding point-wise map. For a TOSCA isometric pair, we sample 700 functional maps of size 50-by-50 with different quality (normalized to the same scale). We measure the penalty of different masks, as an indicator of the quality of the functional map. We also measure the average geodesic error (w.r.t. the direct measurement) of the recovered point-wise map of the corresponding functional map. The correlation between the functional map penalty and the point-wise map quality is visualized as a scatter plot with 700 samples. Compared to the standard and the slanted mask, our resolvent mask applies a smaller penalty to functional maps with good underlying point-wise maps and more heavily penalizes functional maps with bad underlying point-wise maps.

To justify that our resolvent mask is a better regularizer than the standard and the slanted ones, we also measure how the mask penalty relates to the accuracy of the corresponding point-wise map shown in Fig. <ref>. In this experiment, we generate 700 different point-wise maps with different levels of accuracy, then convert them to a functional map representation and measure the penalty added by different masks w.r.t. the average geodesic error computed from the pointwise maps. Each scatter point in Fig. <ref> shows such a test sample.

We can observe that, compared to the standard and the slanted mask, the new mask induces a lower penalty on functional maps with a good quality (i.e., smaller average geodesic error), and penalizes a functional map with larger error more heavily.
This further confirms that using our resolvent mask is more likely to produce a better functional, and ultimately
better pointwise map.

§.§ Parameters

Changing the relative weight $\alpha_4$ of the mask term. We choose 68 different values for $\alpha_4$ in the range of 0 and 10 for this test.

Average error on 284 TOSCA pairs of different masks, each mask using its own optimal weight from Fig. <ref>. Specifically, we set $\alpha_4$ to $10^{-2}$, $10^{-4}$, $10^{-3}$ for the standard, slanted, and our mask respectively.

type     6cAverage error $(\times 10^{-3})$
3c|per-vertex measure     3cdirect measure
Ini     + ICP     1c|+ BCICP    Ini    + ICP    + BCICP
standard     87.26     65.97     1c|41.22     178.8     130.8     76.87

slanted     64.66     50.76     1c|37.45     166.9     131.5     100.0

ours     55.52     43.82     1c|32.48     124.5     90.6     62.33

In our tests, we use $\gamma = 0.5$ and $w = 0.5$ to construct our resolvent mask, and use the default value $\eta = 0.03$ to construct the slanted mask. When comparing the three masks in the initialization setting, i.e., to optimize the energy defined in (<ref>), the weights $\alpha_i$ are set to the same values as reported in [34] but in a relative way (see Appendix <ref> for more details). As for the comparison of the ICP and BCICP refinement settings, the same default parameters and the same number of iterations are used for different masks.

Fig. <ref> shows the results of the changing the weight of the mask term, i.e., $\alpha_4$, while keeping the rest weights fixed on 50 TOSCA shape pairs. We can see that, for all different choices of the weight, our resolvent mask is always better than the standard mask. When the weight is small enough, it seems the slanted mask can achieve a faster decrease of the error than ours. However, the range of the effective weight is smaller than ours besides that the construction is purely heuristic and lacks a theoretical justification.

[trim=0cm 0cm 0cm 0cm,clip,width=1, grid=false]figures_new/eg_centaur.png

(96,17)-90+ ICP

Example. Comparing the quality of the initial maps computed with different masks and after ICP refinement via texture transfer. Each mask uses its own optimal weight from Fig. <ref>.

On the other hand, Fig. <ref> also suggests that different masks have their own preferred choice of the weight $\alpha_4$. Therefore, instead of fixing $\alpha_4$, we set it independently w.r.t. the optimal values on a subset reflected in Fig. <ref>. Specifically, for the standard mask, we set $\alpha_4 = 10^{-2}$, for the slanted mask, we set $\alpha_4 = 10^{-4}$, and for ours, we set $\alpha_4 = 10^{-3}$. The corresponding average errors on the complete TOSCA dataset are reported in Table <ref>. We can see that, even when the parameters are carefully tuned for the other two masks, our resolvent mask still achieves the best accuracy. Fig. <ref> shows a qualitative example.

§.§ Application to non-isometric shape pairs

[trim=0cm 0cm 0cm 0cm,clip,width=1, grid=false]figures_new/eg_shrec.png

Example. Comparing the quality of the initial maps computed with different masks of four SHREC shape pairs via texture transfer.

For 380 SHREC shape pairs, we plot the average error for each pair with different masks, where the result of the standard, slanted and our mask is colored in blue, yellow, and red respectively. Therefore, in the case where the red curve is below the blue or the yellow points, our mask leads to a lower error for these pairs.

[trim=0cm 1cm 0cm 0cm,clip,width=1, grid=false]figures_new/eg_failure.png

Failure case. Here we show a challenging shape pair from SHREC, where all the three masks fail to produce a good map.

To show the usefulness of our resolvent mask on non-isometric shape pairs, we test the 20 FourLeg shapes from the SHREC 2007 dataset [13].
Specifically, we use 10 pairs of descriptors constructed from 4 landmarks to optimize a 120-by-120 functional map using the standard, the slanted, and our resolvent mask.
For a fair comparison, we use the same parameters as the previous tests: we set $\gamma = 0.5$ for our resolvent mask, and set $\eta = 0.03$ for the slanted mask.
Fig. <ref> shows a qualitative example.
We can see that our resolvent mask gives the best initialization.
We also measure the accuracy of the 380 point-wise maps among the 20 shapes on the given landmarks (21 ground-truth landmarks are given for each shape).
The average error for the standard, slanted, and our mask is 0.114, 0.122, and 0.107 respectively.
Fig. <ref> reports the average error for each shape pair.
We observe that resolvent mask gives a better result than the standard mask on 68.4% out of 380 pairs, and outperforms the slanted mask on 69.5% pairs.
The limited quantitative improvement of our mask over the other two masks is due to:
(1) for the shape pairs where our mask significantly outperforms the other two, e.g., as shown in Fig. <ref>, the map quality is measured on only 21 landmarks, where the average error does not fully reflect our improvement.
(2) for some extremely challenging pairs, e.g., a failure case shown in Fig. <ref>, all the maps from different masks have a poor quality, and thus the “relative improvement” on the average error is not informative.
Moreover, we also computed the average error in the case where the mask term is removed from the total energy, i.e., set $\alpha_4 = 0$. In this case, the average error over the complete dataset is 0.112. We can see that, the standard and the slanted mask can have a negative effective in this case, while our resolvent mask still works to improve the map quality to some extent.

§ CONCLUSION, LIMITATIONS AND FUTURE WORK

In this paper, we proposed a new regularizer, the resolvent Laplacian commutativity, for the functional map pipeline.  We first analyzed the limitations of the original Laplacian commutativity term and theoretically justified the effectiveness of our proposed new term. This new regularizer can significantly improve the quality of the computed functional maps and the corresponding recovered point-wise maps before or after refinement.

However, our method also has several limitations that we would like to overcome in future work.  First, our proposed regularizer is well justified on isometric and non-isometric shape pairs, but not on partial shapes, where the ground-truth functional maps can have a different structure. Therefore, it would be interesting to extend the analysis to partial shape pairs.
Second, besides the funnel pattern, in our experiments, we also observed the slanted-diagonal structure of the ground-truth functional map of some non-isometric shape pairs as discussed in [31]. It would be interesting to consider this feature into the mask constructions for further improvement for non-isometric datasets.

Thirdly, the role of the (a, b) introduced in Eq. (5) is not well studied, and we would like to leave this exploration as future work.

Finally, we believe that it would be interesting to study and potentially apply data-driven techniques to
learn the optimal operators for enforcing commutativity across diverse shape pairs.

The authors would like to thank the anonymous reviewers for their valuable comments and helpful suggestions.
Parts of this work were supported by the KAUST OSR Awards No. CRG2017-3426 and CRG2018-3730, a gift from the NVIDIA Corporation, and the ERC Starting Grant No. 758800 (EXPROTEA).

[1]
Aflalo Y., Dubrovina A., Kimmel R.:
Spectral generalized multi-dimensional scaling.
International Journal of Computer Vision 118, 3 (2016),

[2]
Aubry M., Schlickewei U., Cremers D.:
The wave kernel signature: A quantum mechanical approach to shape
In 2011 IEEE international conference on computer vision
workshops (ICCV workshops) (2011), IEEE, pp. 1626–1633.

[3]
Bronstein A. M., Bronstein M. M., Kimmel R.:
Numerical Geometry of Non-Rigid Shapes.
Springer Science & Business Media, 2008.

[4]
Biasotti S., Cerri A., Bronstein A., Bronstein M.:
Recent trends, applications, and perspectives in 3d shape similarity
In Computer Graphics Forum (2016), vol. 35, pp. 87–119.

[5]
Burghard O., Dieckmann A., Klein R.:
Embedding shapes with Green's functions for global shape matching.
Computers & Graphics 68 (2017), 1–10.

[6]
Bogo F., Romero J., Loper M., Black M. J.:
FAUST: Dataset and evaluation for 3D mesh registration.
In Proceedings IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR) (Piscataway, NJ, USA, June 2014), IEEE.

[7]
Choukroun Y., Shtern A., Bronstein A. M., Kimmel R.:
Hamiltonian operator for spectral shape analysis.
IEEE transactions on visualization and computer graphics

[8]
Cuturi M.:
Sinkhorn distances: Lightspeed computation of optimal transport.
In Advances in neural information processing systems (2013),
pp. 2292–2300.

[9]
Dodziuk J.:
Eigenvalues of the Laplacian and the heat equation.
The American Mathematical Monthly 88, 9 (1981), 686–695.

[10]
Ezuz D., Ben-Chen M.:
Deblurring and denoising of maps between shapes.
In Computer Graphics Forum (2017), vol. 36, Wiley Online
Library, pp. 165–174.

[11]
Eynard D., Rodola E., Glashoff K., Bronstein M. M.:
Coupled functional maps.
In 3D Vision (3DV), 2016 Fourth International Conference on
(2016), IEEE, pp. 399–407.

[12]
Gehre A., Bronstein M., Kobbelt L., Solomon J.:
Interactive curve constrained functional maps.
In Computer Graphics Forum (2018), vol. 37, Wiley Online
Library, pp. 1–12.

[13]
Giorgi D., Biasotti S., Paraboschi L.:
Shape retrieval contest 2007: Watertight models track.

[14]
Huang R., Ovsjanikov M.:
Adjoint map representation for shape analysis and matching.
In Computer Graphics Forum (2017), vol. 36, Wiley Online
Library, pp. 151–163.

[15]
Heeren B., Rumpf M., Schröder P., Wardetzky M., Wirth B.:
Splines in the space of shells.
In Computer Graphics Forum (2016), vol. 35, Wiley Online
Library, pp. 111–120.

[16]
Huang Q., Wang F., Guibas L.:
Functional map networks for analyzing and exploring large shape
ACM Transactions on Graphics (TOG) 33, 4 (2014), 36.

[17]
Kovnatsky A., Bronstein M. M., Bronstein A. M., Glashoff K., Kimmel
Coupled quasi-harmonic bases.
In Computer Graphics Forum (2013), vol. 32, pp. 439–448.

[18]
Kovnatsky A., Bronstein M. M., Bresson X., Vandergheynst P.:
Functional correspondence by matrix completion.
In Proceedings of the IEEE conference on computer vision and
pattern recognition (2015), pp. 905–914.

[19]
Kovnatsky A., Glashoff K., Bronstein M. M.:
MADMM: a generic algorithm for non-smooth optimization on
In Proc. ECCV (2016), Springer, pp. 680–696.

[20]
Litany O., Rodolà E., Bronstein A. M., Bronstein M. M., Cremers
Non-rigid puzzles.
In Computer Graphics Forum (2016), vol. 35, Wiley Online
Library, pp. 135–143.

[21]
Litany O., Rodolà E., Bronstein A. M., Bronstein M. M.:
Fully spectral partial shape matching.
In Computer Graphics Forum (2017), vol. 36, Wiley Online
Library, pp. 247–258.

[22]
Mandad M., Cohen-Steiner D., Kobbelt L., Alliez P., Desbrun M.:
Variance-minimizing transport plans for inter-surface mapping.
ACM Transactions on Graphics (TOG) 36 (2017), 14.

[23]
Minakshisundaram S., Pleijel A.:
Some properties of the eigenfunctions of the Laplace-operator on
Riemannian manifolds.
Canadian J. Math 1, 8 (1949).

[24]
Melzi S., Rodolà E., Castellani U., Bronstein M. M.:
Localized manifold harmonics for spectral shape analysis.
In Computer Graphics Forum (2018), vol. 37, Wiley Online
Library, pp. 20–34.

[25]
Nogneng D., Melzi S., Rodolà E., Castellani U., Bronstein M.,
Ovsjanikov M.:
Improved functional mappings via product preservation.
In Computer Graphics Forum (2018), vol. 37, Wiley Online
Library, pp. 179–190.

[26]
Nogneng D., Ovsjanikov M.:
Informative descriptor preservation via commutativity for shape
Computer Graphics Forum 36, 2 (2017), 259–267.

[27]
Neumann T., Varanasi K., Theobalt C., Magnor M., Wacker M.:
Compressed manifold modes for mesh processing.
In Computer Graphics Forum (2014), vol. 33, Wiley Online
Library, pp. 35–44.

[28]
Ovsjanikov M., Ben-Chen M., Solomon J., Butscher A., Guibas L.:
Functional Maps: A Flexible Representation of Maps
Between Shapes.
ACM Transactions on Graphics (TOG) 31, 4 (2012), 30.

[29]
Ovsjanikov M., Corman E., Bronstein M., Rodolà E., Ben-Chen M.,
Guibas L., Chazal F., Bronstein A.:
Computing and processing correspondences with functional maps.
In ACM SIGGRAPH 2017 Courses (2017), SIGGRAPH '17,
pp. 5:1–5:62.

[30]
Poulenard A., Skraba P., Ovsjanikov M.:
Topological function optimization for continuous shape matching.
In Computer Graphics Forum (2018), vol. 37, Wiley Online
Library, pp. 13–25.

[31]
Rodolà E., Cosmo L., Bronstein M. M., Torsello A., Cremers D.:
Partial functional correspondence.
In Computer Graphics Forum (2017), vol. 36, Wiley Online
Library, pp. 222–236.

[32]
Rodolà E., Moeller M., Cremers D.:
Point-wise map recovery and refinement from functional
In Proc. Vision, Modeling and Visualization (VMV) (2015).

[33]
Rustamov R., Ovsjanikov M., Azencot O., Ben-Chen M., Chazal F., Guibas
Map-based exploration of intrinsic shape differences and variability.
ACM Transactions on Graphics (TOG) 32, 4 (July 2013),

[34]
Ren J., Poulenard A., Wonka P., Ovsjanikov M.:
Continuous and orientation-preserving correspondences via functional
ACM Transactions on Graphics (TOG) 37, 6 (2018).

[35]
Reed M., Simon B.:
Methods of Modern Mathematical Physics I: Functional Analysis,
revised and enlarged ed.
Academic Press, San Diego, 1980.

[36]
Sauvigny F.:
Partial Differential Equations 1: Foundations and Integral
Springer, Berlin, Heidelberg, 2006.

[37]
Sauvigny F.:
Partial Differential Equations 2: Functional Analytic Methods.
Springer, Berlin, Heidelberg, 2006.

[38]
Solomon J., De Goes F., Peyré G., Cuturi M., Butscher A., Nguyen
A., Du T., Guibas L.:
Convolutional wasserstein distances: Efficient optimal transportation
on geometric domains.
ACM Transactions on Graphics (TOG) 34, 4 (2015), 66.

[39]
Solomon J., Peyré G., Kim V. G., Sra S.:
Entropic metric alignment for correspondence problems.
ACM Transactions on Graphics (TOG) 35, 4 (2016), 72.

[40]
Tam G. K., Cheng Z.-Q., Lai Y.-K., Langbein F. C., Liu Y., Marshall D.,
Martin R. R., Sun X.-F., Rosin P. L.:
Registration of 3D point clouds and meshes: a survey from rigid to
IEEE TVCG 19, 7 (2013), 1199–1217.

[41]
Van Kaick O., Zhang H., Hamarneh G., Cohen-Or D.:
A survey on shape correspondence.
In Computer Graphics Forum (2011), vol. 30, pp. 1681–1707.

[42]
Vestner M., Litman R., Rodolà E., Bronstein A., Cremers D.:
Product manifold filter: Non-rigid shape correspondence via kernel
density estimation in the product space.
In Proc. CVPR (2017), pp. 6681–6690.

[43]
Wang L., Gehre A., Bronstein M. M., Solomon J.:
Kernel functional maps.
In Computer Graphics Forum (2018), vol. 37, Wiley Online
Library, pp. 27–36.

[44]
Wang F., Huang Q., Guibas L. J.:
Image co-segmentation via consistent functional maps.
In Proceedings of the IEEE International Conference on Computer
Vision (2013), pp. 849–856.

[45]
Wang Y., Liu B., Zhou K., Tong Y.:
Vector field map representation for near conformal surface
In Computer Graphics Forum (2018), vol. 37, Wiley Online
Library, pp. 72–83.