# Partial Transport for Point-Cloud Registration

Yikun Bai1, Huy Tran1, Steven B. Damelin2, Soheil Kolouri1, Senior Member IEEE

1 Department of Computer Science, Vanderbilt University, Nashville, TN  2
Michigan, Ann Arbor, MI

###### Abstract

Point cloud registration plays a crucial role in various fields, including
robotics, computer graphics, and medical imaging. This process involves
determining spatial relationships between different sets of points, typically
within a 3D space. In real-world scenarios, complexities arise from non-rigid
movements and partial visibility, such as occlusions or sensor noise, making
non-rigid registration a challenging problem. Classic non-rigid registration
methods are often computationally demanding, suffer from unstable performance,
and, importantly, have limited theoretical guarantees. The optimal transport
problem and its unbalanced variations (e.g., the optimal partial transport
problem) have emerged as powerful tools for point-cloud registration,
establishing a strong benchmark in this field. These methods view point clouds
as empirical measures and provide a mathematically rigorous way to quantify
the ‘correspondence’ between (the transformed) source and target points. In
this paper, we approach the point-cloud registration problem through the lens
of optimal transport theory and first propose a comprehensive set of non-rigid
registration methods based on the optimal partial transportation problem.
Subsequently, leveraging the emerging work on efficient solutions to the one-
dimensional optimal partial transport problem, we extend our proposed
algorithms via slicing to gain significant computational efficiency, resulting
in fast and robust non-rigid registration algorithms. We demonstrate the
effectiveness of our proposed methods and compare them against baselines on
various 3D and 2D non-rigid registration problems where the source and target
point clouds are corrupted by random noise.

###### Index Terms:

Optimal transport, point cloud registration

##  I Introduction

Data acquisition and the increasing interest in augmented and virtual reality
have led to an explosion of volumetric data, as point clouds exemplify. Point
cloud data is prevalent in numerous applications, including robotics [2, 117,
93], autonomous driving [57, 89, 94], medical imaging [116, 102] and computer
graphics [34, 132]. In many of these applications, the captured point clouds
correspond to noisy observations of an object/scene undergoing different
deformations. One of the core challenges in these applications is to perform
point cloud registration, which refers to finding a transformation that aligns
or partially aligns the source and target point sets. At a high level, any
point cloud registration algorithm must solve two problems concurrently: 1)
finding accurate correspondences between the points in the source and target
point clouds (implicitly or explicitly), and 2) modeling the deformation to
match the corresponding source and target points. The existing methods then
propose different correspondence estimation algorithms [49, 35, 120, 72, 38,
40, 28] and/or propose novel deformation modeling [10, 133, 85, 23, 127].

The registration/deformation map (i.e., the transformation) could be rigid
[52, 53, 19, 120, 56], only involving translation and rotation, or non-rigid
(e.g., an affine transformation or other nonlinear deformation) [23, 29, 58,
79]. Most existing works in the literature have focused on the rigid
registration of point clouds, as it is a more prevalent problem in classic
computer vision tasks such as Simultaneous Localization And Mapping (SLAM)
[83, 4]. The core innovations in these approaches are often concerned with
finding the right correspondences between the points. For instance, the
classic Iterative Closest Point (ICP) algorithm [52, 53, 19, 139, 36] relies
on nearest-neighbor correspondences as measured via the Euclidean distance
between points. However, the Euclidean distance between points might not be a
reliable measure of proximity if, for instance, the point clouds are not
approximately aligned or the point clouds are noisy. To establish better
correspondences, many have looked into improving the similarity measures by
defining features/descriptors for point cloud data[31, 30, 97, 135, 130, 131,
90] that not only capture the location of a point but also encode the local
geometry of the object in the vicinity of that point (e.g., curvature). These
methods rely on nearest neighbor matching in an explicit or implicit feature
space instead of the raw input space and are closely related to kernel methods
[120, 79, 58, 85]. Moreover, to obtain robustness and avoid false
correspondences, these methods often use the Random Sample Consensus (RANSAC)
algorithm [32, 43] or its variations [119, 24]. Unfortunately, however, RANSAC
significantly increases the computational cost of the registration algorithm
as it requires running the correspondence problem multiple times for different
random subsets of the data.

In many real-world applications, however, the deformation between two sets of
point cloud data is inherently non-rigid. For instance, in medical imaging,
the point cloud data could come from the surface of a tissue, which can
undergo large non-rigid deformations. Such nonlinear deformations are
generally modeled using two categories of approaches, parametric [23, 129,
120] and non-parametric approaches [85, 86, 127]. In the parametric
approaches, the deformation is characterized via a parametric function, e.g.,
parameters of an affine transformation [120] or Thin Plate Spline (TPS)
parameters [23], and they are optimized to minimize the expected distance
between the corresponding source and target points. On the other hand, the
non-parametric approaches directly calculate the displacement (velocity)
between source points and their corresponding target points. The existing
approaches for non-parametric non-rigid point cloud registration vary in how
to estimate the velocity of each point and how to regularize the velocity
vector field for coherency and smoothness[85, 134]. In this paper, we consider
parametric non-rigid registration of point clouds and utilize Optimal Partial
Transport (OPT) [16, 41], an instance of the unbalanced optimal transport [21,
39], as a unifying framework to achieve this task.

Solving the correspondence between points in the source and target point
clouds is closely related to the celebrated optimal transport (OT) problem
[122, 121, 91]. In short, treating point clouds as empirical distributions and
given a transportation cost, the OT problem seeks the optimal assignment
between the samples to minimize the expected transportation cost between the
assigned samples. This principle has led to many OT-based point-cloud
registration algorithms. However, a major limitation of OT is the mass
preservation assumption, which requires all points in the source to be matched
to all the points in the target point cloud. The mass preservation assumption
limits the application of OT to problems where the points must be partially
matched, e.g., noisy or occluded point clouds and, in general, partial
registration problems [28]. Various ideas have been recently developed to
allow the application of OT for partial registration [127, 98]. For instance,
some dynamically infer the mass of each particle (i.e., the importance of a
point in the point cloud) [114] so that the OT problem can ignore the
particles with zero/small mass. Others have looked into defining outlier bins
and solving the OT problem while allowing for matching points to the outlier
bin [28]. This latter idea is deeply rooted in the unbalanced optimal
transport problem, which allows for the creation and destruction of mass. In
this paper, we use OPT and its sliced variation as an excellent match for
partial and robust registration problems and demonstrate their performance in
non-rigid partial registration of point clouds.

Contributions: Our specific contributions are as follows.

  1. 1.

A robust and unifying framework for parametric non-rigid registration using
OPT.

  2. 2.

Proposing a sliced-OPT framework for accelerated non-rigid registration
between large-scale point clouds.

  3. 3.

Demonstrating the performance of the proposed framework on various noisy point
cloud registration problems and benchmarking them against baselines.

In what follows, we first review the related work in the literature in Section
II. Next, in section III, we introduce some basic background of optimal
transport and optimal partial transport. Following this section, we introduce
our OPT-based methods in section IV and demonstrate our methods in 2D and 3D
experiments in V. Finally, we discuss the limitations and future direction in
section VI. In appendix A, we introduce all the notations in section A.

##  II Related Work

The prior work on point cloud registration could be categorized based on the
type of registration, e.g., rigid versus non-rigid, or based on how
correspondences are calculated/found, e.g., feature-based approaches. A large
body of existing work takes an alternating optimization approach, where they
jointly estimate the correspondences and the transformation iteratively .
These methods update the correspondences assuming the transformation is fixed
and then update the transformation assuming that correspondences are fixed.

Alternatively, the recent and emerging work from the deep learning community
on point cloud registration focus on 1) feature learning for calculating
robust correspondences [68, 96], and 2) end-to-end estimation of the
transformation [74].

Our proposed approach in this paper is to directly optimize the transformation
by minimizing an optimal partial transportation metric between the transformed
source and the target point clouds. Hence, it is closely related to the
alternating optimization approaches, with the main difference that the
calculation of correspondences is implicitly done when calculating the OPT
metric. Below we briefly overview the rigid and non-rigid registration
approaches, the emerging deep learning-based registration methods, and lastly,
the existing OT-based registration algorithms. Finally, we differentiate our
approach from the existing work and state our contributions.

###  II-A Rigid/Affine Point Cloud Registration

In rigid registration, one seeks a transformation composed of only a rotation
and a translation to match the source and target point clouds. The dominant
framework on this front is the seminal works of [10] and [137] on Iterative
Closet Point (ICP), and their numerous extensions [103, 136]. For a given
initial transformation and a distance between the source and the target points
(e.g., the Euclidean distance), the classic ICP algorithm finds the
correspondences via the nearest neighbor assignment, and then it updates the
transformation to be the least-squares minimizer. Given the correspondences,
and using the Euclidean distance, finding the rotation and translation
parameters has a closed-form solution. Through these alternating schema, ICP
is guaranteed to converge to a locally optimal transformation. The classic ICP
algorithm, see for example [52, 53, 19, 139, 36] relies on a good initial
transformation that provides nearly aligned point clouds so that the
correspondences calculated with the nearest neighbor assignments are reliable.
Moreover, the classic ICP algorithm fails in the presence of noise and when
dealing with partial registration (i.e., when only a subset of points should
be matched). Lastly, the classic ICP algorithm has a linear convergence rate,
which can be too slow for various registration applications (e.g., for
simultaneous localization and mapping (SLAM) in robotics).

To address the limitations of ICP, a large body of work has been developed in
the past two decades. A critical idea in these works is to replace the binary
assignment in ICP with relaxed soft assignments and provide a generalized
version of the classic ICP algorithm. For instance, the Robust Point Matching
(RPM) algorithm [101, 46] and its extensions. The footprints of the idea of
using soft assignments can also be seen in numerous approaches that pose the
registration problem as a _maximum likelihood_ estimation problem [61, 25, 77,
80]. An interesting connection to the OT theory here is that the entropy
regularized Kantorovich problem [26], also provides a soft assignment (i.e.,
the entropy regularized transport plan) that can be utilized for the iterative
matching problem in registration. This idea is referred to as the Wasserstein
Procrustes problem [47, 100, 59]. For instance, Grave et al. [47] use this
idea to align two high-dimensional point clouds for Natural Language
Processing (NLP) applications. Another common idea for robustifying the
_maximum likelihood_ based registration algorithms against noise, is to
introduce an extra distribution term to account for noise/outliers [101]
(i.e., using dummy outlier bins). An interesting connection to the OT theory,
here, is the idea of introducing dummy source and sink points in the optimal
partial transport theory (for instance, see [17]) to allow for destruction and
creation of mass in the source and target distributions, respectively. These
connections motivated our approach to use OPT and sliced OPT for noisy point
cloud registration.

###  II-B Nonrigid registration

Non-rigid registration techniques align two shapes that undergo non-rigid
deformation, such as bending, stretching, and twisting. Unlike rigid
registration, which only involves a single rotation and translation, non-rigid
registration must often determine a deformation field [29]. A proper
representation of the deformation fields must be chosen to balance sufficient
expressiveness for accurate alignment and computational efficiency. Generally,
the deformation fields can be divided into two categories: non-parametric and
parametric models. In the following, we introduce some widely used
representation techniques for deformation.

Cohrent point drift (CPD) method. The coherent point drift (CPD) method is
proposed by [86, 85]. In this method, the point clouds are regarded as
Gaussian mixture models (GMMs), and they have aligned via minimizing the
discrepancies between them111In different literature, CPD is attributed into
RKHS-based method, non-parametric approach, or probabilistic (GMM) model-based
method. Here, to avoid controversial classifications, we treat it as a
separate category.. In this work, the deformation can be formulated as

$\hat{y}=x+\sum_{k\in\mathcal{I}}\alpha^{k}\phi_{G}(x,x^{k}),$

where $\mathcal{I}\subset[1:n]$ is a subset of (index of) source points
$\phi_{G}$ is the Gaussian kernel, which is introduced in Appendix A. Hirose
et al. [50] proposed a variant of the CPD method in a Bayesian setting. The
new method is proved to have a convergence guarantee under variational
Bayesian inference. Wang et al. [127] use the Gaussian kernel matrix in the
CPD method to regularize the displacement parameters, and the correspondence
is estimated by optimal partial transportation.

In addition, deep learning techniques have been applied in this framework. For
example, [125], proposes the so-called CPD-network method, which can learn a
displacement field function to estimate a certain geometric transformation
from a training dataset. Thus, these techniques can predict the desired
geometric transformation to align previously unseen pairs of data without an
additional optimization process. Similarly, [78] combines CPD with a certain
transformer network technique in this method.

Thin plate spline-based methods. This type of method uses a deformation field
defined by the so-called thin plate spline function [14, 123, 37] in
$\mathbb{R}^{3}$ space. In particular,

$f(x)=\sum_{k=1}^{K}\alpha^{k}\phi_{T}(x^{k},x)+B^{T}x+\beta,$

where $\phi_{T}$ is TPS kernel, $B\in\mathbb{R}^{D\times
D},\beta\in\mathbb{R}^{D}$, we refer appendix A and subsection B in section IV
for details. The above formulation is derived from a square error minimization
problem with a second derivative regularization term.

In particular, [23] develops the popular TPS-RPM (thin plate spline-robust
point-matching) algorithm which uses the TPS as the non-rigid spatial mapping
and the soft assign for the correspondence. [35] traces back the registration
problem to the solution of a system of nonlinear equations that directly gives
the parameters of the TPS function. Thus the proposed method recovers
deformation without established correspondence. In addition, [55] developed a
new technique that chooses the control points based on the feature
correspondence between two surfaces.

Kernel correlation and probability model. This type of method employs kernel
correlation function to ascertain correspondences between point clouds or
features. Typically, these point clouds are often represented using
probability models. The seminal work proposed by [120] set the stage for
subsequent advancements in this type of approach. This study leverages the
kernel correlation to the similarity between the source and target point set
and formalizes the registration challenge as a kernel correlation optimization
issue. Building on this foundation, [58] employed the Gaussian Mixture Model
(GMM) to characterize each data point cloud. Their approach to the
registration problem aimed to minimize the discrepancy, especially the L2
distance, between the two GMMs. Moreover, they integrated Thin Plate Splines
(TPS) and Gaussian radial basis functions for non-rigid transformation models.
The CPD method, cited in [85, 79], is a subsequent development. This method
adopts the GMM model to represent the point cloud and determines
correspondences by solving a maximum likelihood problem. In recent times,
point cloud kernel correlation has been merged with deep learning models. For
instance, [113] enhances the Pointnet \- a stage-of-art neural network model
for semantic learning on 3D point clouds by integrating a kernel correlation
layer, allowing the model to harness the local geometric structures of point
clouds.

ICP-FFD The ICP algorithm is widely used for rigid registration but has
numerous limitations, as discussed earlier. To extend its capabilities,
various techniques have been developed to model local deformations. Abdelmunim
and Farag [1] provide a flexible non-rigid registration method by combining
ICP with Freeform deformation (FFD). FFD defines a lattice of control points,
which can be moved to control the deformation of the underlying point cloud.
The ICP-FFD method iteratively refines the alignment between point clouds by
minimizing the distance between corresponding points while adjusting the
control points’ positions to account for non-rigid transformations. In
practice, this approach has been applied widely, including in medical image
registration and facial expression recognition. Another combination is ICP and
Gaussian Process Regression (GPR) [15]. GPR is a powerpoint-cloud learning
technique that can model the non-rigid deformation between two point clouds as
a smooth, non-linear function. Like ICP-FFD, the ICP-GPR algorithm also
iteratively refines the alignment between the point clouds using ICP while
modeling the non-rigid deformation with GPR. This approach can handle large
deformations while providing an explicit representation of the deformation
function, which helps with interpretability. ICP can also be combined with an
energy function to handle non-rigid registration [72].

The ICP algorithm is then used iteratively to minimize the energy function,
resulting in a registration that considers both the point cloud alignment and
deformation smoothness.

Large deformations. In non-rigid registration, large deformations refer to a
class of transformations that may change the overall structure of the point
cloud significantly. Isometric deformations are a specific type of large
deformations where the intrinsic geometric properties, such as geodesic
distances between points, are preserved, whereas the extrinsic properties,
like Euclidean distances, may change. For instance, in their seminal work,
Adams et al. [54] proposed a large deformation framework based on the theory
of functional maps. Their proposed method establishes correspondences between
two shapes by aligning their Laplace-Beltrami eigenfunction bases, which are
intrinsic geometric operators invariant under isometric deformations. The
method consists of several key steps, including computing the Laplace-Beltrami
operators, obtaining eigenbases for representing shapes, representing shapes
as functions on their respective surfaces, finding a transformation, and
establishing point correspondence. The method has been extended in various
ways with additional geometric or topological information to handle complex
deformations [112, 3]. It has been successfully applied widely in practice for
registration, object recognition[22], and computer animation [71, 118].

###  II-C Deep Learning Based Registration

Deep learning has rapidly emerged in recent years as a powerful tool for
solving many computer vision tasks, including point cloud registration [138].
Recent advances in the field have led to a plethora of neural network-based
approaches that offer improved speed [96], accuracy [124], and robustness [6]
for registration. In fact, these methods demonstrate great promise in
enhancing the pipeline and overcoming challenging scenarios where traditional
methods fall short.

Existing approaches often involve finding correspondences in 3 key steps:
feature selection, matching, and motion estimation [10] [104]. First, for
feature selection, the traditional methods have evolved from using basic
Cartesian coordinates [10] [137] to hand-crafted descriptors capturing complex
geometric properties [60] [51] [90]. These traditional methods are often sub-
optimal for large datasets. In sharp contrast, deep learning approaches like
PointNet [96] and its extensions [97] offer automated, robust feature
extraction by leveraging neural networks [68]. Regarding the second step,
matching, traditional methods, like ICP, struggle with issues such as
sensitivity to initialization, convergence to local minima, and poor
performance in partial matching. while Random Sample Consensus (RANSAC)
requires balancing accuracy and computational cost. Deep learning-based
solutions like Deep Closest Point (DCP) [126] offer robust and efficient
matching by leveraging PointNet-based architectures and attention mechanisms.
For the third step, motion estimation, traditional methods for motion
estimation, such as least square linear regression and Singular Value
Decomposition (SVD), are sensitive to noise and limited to rigid
transformations. In contrast, deep learning-based solutions like FlowNet3D
[74] offer robust and flexible motion estimation by learning point-wise
correspondence flow between point clouds.

Correspondence-free methods There are also deep learning approaches that
completely bypass finding correspondences by regressing motion parameters
using global features. For example, PointNetLK [6] and PCRNet [105] both
utilize deep learning techniques to effectively estimate the relative pose
between two point clouds without explicitly finding point correspondences.
PointNetLK [6] extends the original PointNet [96] architecture by
incorporating a differentiable Lucas-Kanade (LK) layer to estimate the
transformation parameters iteratively. PCRNet [105], on the other hand,
employs a Siamese network architecture to learn global point cloud features
and regress the transformation parameters directly. PCRNet leverages the
robustness of PointNet for feature extraction and incorporates a fully
connected regression network to predict the 6-DoF transformation. Some other
methods are completely end-to-end, such as DeepVCP [76], allowing prediction
of the correspondences on raw data, without pre- or post-processing.

Despite their effectiveness, deep learning-based supervised approaches have
several significant limitations. For instance, they require a large amount of
annotated training data, which could be prohibitively expensive with scale.
Our work in this paper is orthogonal to the research in the deep learning
community, as our method is also compatible with these approaches, and can be
integrated into learning-based pipelines.

###  II-D OT Based Point Set Registration

Optimal transport (OT) is a mathematical approach that can compare two
probability measures. In the context of point clouds, they can be treated as
discrete probability measures, and the OT problem can be utilized to determine
the optimal (soft) correspondences between two sets of points. These
correspondences aim to minimize the expected transportation cost, which is the
cost of moving from one point to another. After obtaining the optimal
correspondences between the points in the two point clouds, it becomes
possible to compute the transformation that aligns the two point clouds [40].
This transformation can be rigid, such as translation, rotation, scaling, or
non-rigid, like general deformations.

For example: [95] proposed a method to estimate the scene flow on point clouds
using optimal transportation and deep features. On synthetic and real-world
datasets, they demonstrate that the new method can perform as well as the
state-of-art existing methods while requiring much fewer parameters and
without multi-scale analysis. [81] proposes a learning framework for
predicting correspondences of 3D point cloud registration by transferring
point-wise matching and structural matching into a Wasserstein distance-based
and a Gromov-Wasserstein distance-based optimizations, respectively. The
proposed framework can accurately predict correspondences of 3D point cloud
registration and achieve state-of-the-art performance on several benchmarks.
In [38], they use entropic regularized optimal transportation and the smooth
shells method to estimate the unsupervised 3D shape correspondence.

In addition to classical OT in the balanced setting, as mentioned above,
unbalanced OT has been widely applied in shape registration problems too.
[114] investigates the use of Hellinger Kantorovich distance (the so-called
robust optimal transport in the paper) for shape matching. The TPS-RPM (thin
plate sprine-robust point-matching) algorithm proposed by [23] essentially
applied optimal partial transport distance to estimate the correspondence for
the non-rigid registration problem. Similarly, [98] utilize optimal partial
transport via a hard marginal constraint for solving non-rigid problem.
Furthermore, sliced-OPT based approaches have been proposed recently [11] [8].

OT in Deep Learning: Wang and Solomon (2019) proposed Deep Closest Point (DCP)
[126], a deep learning-based approach for point cloud registration that
employed the Sinkhorn algorithm to minimize the regularized OT distance
between two point sets. DCP utilized PointNet to learn local geometric
features and incorporated a differentiable module based on the Sinkhorn
algorithm for the alignment process. This end-to-end trainable architecture
outperformed traditional ICP-based methods. Since the introduction of DCP,
several follow-up works have emerged, addressing different challenges and
further improving the performance of point cloud registration methods that
combine OT and deep learning. Yuan et al. (2020) presented DeepGMR [133],
which used deep learning to compute a Gaussian Mixture Representation (GMR) of
point clouds and applied OT to align the GMRs. This approach reduced
computational complexity while maintaining high registration accuracy. Yew and
Lee (2020) introduced RPM-Net [131], a method that combined deep learning and
the Robust Point Matching (RPM) framework. RPM-Net employed a learned affinity
matrix based on deep features and solved the registration problem using a
differentiable OT-based module. These works demonstrate the potential of
combining OT and deep learning for point cloud registration, leading to
improved performance over traditional methods. However, further research is
needed to address challenges such as computational complexity, partial overlap
constraint, and lack of interpretability.

##  III Background

###  III-A Optimal Transport

Optimal transport (OT) theory, pioneered by Monge [84] and Kantorovich [62,
63], studies the most cost-efficient way to move mass from one measure to
another measure. It has attracted abundant attention in data science,
statistics, machine learning, signal processing and computer vision.

####  III-A1 Classic Optimal Transport

The classic optimal transport, known as Kantorovich formulation, is defined as
follows: let $\mathcal{P}(\mathbb{R}^{D})$ denote the set of all Borel
probability measures defined in $\mathbb{R}^{D}$. Given two probability
measures, $\mu,\nu\in\mathcal{P}(\mathbb{R}^{D})$, and
$c:(\mathbb{R}^{D})^{2}\to\mathbb{R}_{+}$ a lower semi-continuous function
denoting the transportation cost, the optimal transport problem is defined as:

$\displaystyle\text{OT}(\mu,\nu):=\inf_{\gamma\in\Pi(\mu,\nu)}c(\hat{y},y)d\gamma(\hat{y},y)$
(1)

where $\Pi(\mu,\nu)\subset\mathcal{P}((\mathbb{R}^{D})^{2})$ is the set of
joint distribution whose marginals are $\mu,\nu$, respectively. We
mathematically denote the two marginals of $\gamma$ as
$(\pi_{1})_{\\#}\gamma=\mu,(\pi_{2})_{\\#}\gamma=\nu$, where $\pi_{1}$,
$\pi_{2}$ are canonical projection maps, and for any (measurable) function
$f:(\mathbb{R}^{D})^{2}\to\mathbb{R}^{D}$, $f_{\\#}\gamma$ is the push-forward
measure defined as $f_{\\#}\gamma(A)=\gamma(f^{-1}(A))$ for any Borel set
$A\subseteq\mathbb{R}^{D}$.

Regarding the point cloud registration scenario, we use $\mu,\nu$ to represent
the transformed source point cloud and target point cloud, that is, we set
$\mu,\nu$ as empirical distributions:
$\mu=\frac{1}{N}\sum_{n=1}^{N}\delta_{\hat{y}^{n}},\nu=\frac{1}{M}\sum_{m=1}^{M}\delta_{y^{m}}$,
where $\delta_{\hat{y}}$ is the Dirac delta function,
$\\{\hat{y}^{n}\in\mathbb{R}^{D}\\}_{n=1}^{N}$ and
$\\{y^{m}\in\mathbb{R}^{D}\\}_{m=1}^{M}$ are distinct points, the above
formulation becomes:

$\displaystyle
OT(\mu,\nu):=\inf_{\gamma\in\Gamma(\frac{1}{N}1_{N},\frac{1}{M}1_{M})}\sum_{n=1}^{N}\sum_{m=1}^{M}c(\hat{y}^{n},y^{m})\gamma_{n,m}$
(2)

where $c(\hat{y}^{n},y^{m})$ denote the distance between points
$(\hat{y}^{n},y^{m})$, $\gamma$ defines an “soft” correspondence between the
two point clouds.

By OT’s theory, the minimizer for problem (1) (and (2)) always exists, and we
can replace $\inf$ by $\min$. Intuitively, for every feasible joint
distribution $\gamma$ in problem (2) (i.e. $\gamma$ satisfies the constraints
in the problem), $\text{supp}(\gamma)$ defines a correspondence between
$\mu,\nu$ and $\gamma$ describes a transportation plan from
$\\{\hat{y}^{n}\\}_{n=1}^{N}$ to $\\{y^{m}\\}_{m=1}^{M}$. Indeed, for each
$(n,m)\in[1:N]\times[1:M]$, the value of $\gamma_{n,m}$ denotes the amount of
mass that will be moved from $\hat{y}^{n}$ to $y^{m}$ and thus
$\gamma_{n,m}=0$ denotes that there is no correspondence between $\hat{y}^{n}$
and $y^{m}$. When a transportation plan transport masses from $\hat{y}^{n}$ to
at least two $y$ points, say $y^{m_{1}},y^{m_{2}}$, (i.e.,
$\gamma_{n,m_{1}},\gamma_{n,m_{2}}>0$), we say there is mass splitting in this
transportation plan $\gamma$ or the correspondence between $\mu$ and $\nu$ is
soft. Mass splitting is allowed in Kantorovich’s formulation of the OT
problem; however, in some particular cases, for example, when $N=M$, there is
no need to consider such a plan. In particular, when $M=N$ and all samples
have uniform mass $1/N$, by OT’s theory (see page 5 in [122]), the optimal
$\gamma$ is an $N\times N$ permutation matrix (i.e.,
$\gamma_{n,m}\in\\{0,1\\},\forall(n,m)$ and for each row and column, there
exists only one $\gamma_{n,m}=1$), which implies that a one-to-one mapping
induces the OT plan. Therefore, the above Kantorovich’s formulation is
equivalent to the following so-called Monge formulation:

$\displaystyle MOT(\mu,\nu):=\inf_{\begin{subarray}{c}L:[1:n]\to[1:n]\\\
L\text{ is }1-1\end{subarray}}\sum_{n=1}^{N}c(\hat{y}^{n},y^{L(n)}))$ (3)

and the optimal $L$ is called Monge’s mapping. Regarding the registration
scenario, $L$ defines a one-by-one correspondence between the source and
target point clouds.

When $N\neq M$, (for convenience, say $N<M$), feasible $\gamma$s would have
mass splitting, and Monge’s mapping does not exist.

In shape registration tasks, we still prefer hard correspondences; we can use
barycentric projection to get an approximation of the Monge mapping. In
particular, we find the optimal transportation plan $\gamma^{*}$ for
$OT(\mu,\nu)$ and we call
$\hat{\nu}:=\frac{1}{N}\sum_{n=1}^{N}\delta_{((M\gamma Y)[n,:])^{T}}$ (where
$Y=[y^{1},\ldots y^{M}]^{T}$) as the barycentric projection [5, 7] of $\nu$
with respect to $\mu$, which is an $N-$points empirical distribution and can
be regarded as the “representation” of $\nu$. Furthermore, if $N=M$ (i.e.
Monge’s mapping exists), then

$OT(\mu,\nu)=\sum_{n=1}^{N}c(\hat{y}^{n},(M\gamma Y)[n,:]^{T}),$

where $\gamma Y\in\mathbb{R}^{N\times D}$ is matrix multiplication, $M$ is a
scalar, and $(M\gamma Y)[n:]$ is $n^{th}$ row of matrix $M\gamma Y$. Thus, in
practice, the Monge mapping between $\mu$ and $\hat{\nu}$ can be used to
approximate the optimal transportation plan $\gamma$ for $OT(\mu,\nu)$.

Figure 1: In the first column $(N=M)$, the optimal transportation plan
provides a hard assignment where every point has at most one correspondence.
In the second column $(N<M)$, the optimal transportation plan induces a soft
assignment, which can be turned into a hard assignment with the barycentric
projection of $\nu$, as shown in the third column.

####  III-A2 Optimal Partial Transport

The OPT problem, in addition to mass transportation, allows mass destruction
at the source and mass creation at the target. Here the mass destruction and
creation penalty will be linear. Let $\mathcal{M}_{+}(\mathbb{R}^{D})$ denote
the set of all positive Radon measures defined on $\mathbb{R}^{D}$, suppose
$\mu,\nu\in\mathcal{M}_{+}(\mathbb{R}^{D})$, and $\lambda\geq 0$, the OPT
problem is then defined as:

$\displaystyle\text{OPT}_{\lambda}(\mu,\nu):=\inf_{\gamma\in\Gamma_{\leq}(\mu,\nu)}\int
c(\hat{y},y)\,d\gamma+\lambda(|\mu|+|\nu|-2|\gamma|),$ (4)

for
$\Gamma_{\leq}(\mu,\nu):=\\{\gamma\in\mathcal{M}_{+}((\mathbb{R}^{D})^{2})|\pi_{1\\#}\gamma\leq\mu,\pi_{2\\#}\gamma\leq\nu\\}$
where $\pi_{1\\#}\gamma\leq\mu$ denotes that for any Borel set
$A\subseteq\mathbb{R}^{D}$, $\pi_{1\\#}\gamma(A)\leq\mu(A)$, and we say
$\pi_{1\\#}\gamma$ is dominated by $\mu$, analogously for
$\pi_{2\\#}\gamma\leq\nu$; the notation $|\mu|$ denotes the total mass of
measure $\mu$. When $c(\hat{y},y)$ is a metric,
$\text{OPT}_{\lambda}(\cdot,\cdot)$ is a metric on
$\mathcal{M}_{+}(\mathbb{R}^{D})$ (see [21, Proposition 2.10], [92,
Proposition 5], [69, Section 2.1] and [18, Theorem 4]).

In point cloud registration problem, we set $\mu,\nu$ to be empirical
measures, $\mu=\sum_{n=1}^{N}\delta_{\hat{y}^{n}}$ and
$\nu=\sum_{m=1}^{M}\delta_{y^{m}}$, with distinct
$\hat{y}^{n},y^{m}\in\mathbb{R}^{D},\forall n,m$ to denote the estimated and
target point clouds. The OPT problem (4), denoted as
$\text{OPT}(\\{\hat{y}^{n}\\}_{n=1}^{N},\\{y^{m}\\}_{m=1}^{M})$, can be
rewritten as:

$\displaystyle\text{OPT}(\\{\hat{y}^{n}\\}_{n=1}^{N},\\{y^{m}\\}_{m=1}^{M})$
$\displaystyle:=$
$\displaystyle\min_{\pi\in\Gamma_{\leq}(1_{N},1_{M})}\sum_{n,m}c(\hat{y}^{n},y^{m})\gamma_{n,m}+\lambda(N+M-2\sum_{n,m}\gamma_{n,m})$
(5)

where

$\Gamma_{\leq}(1_{N},1_{M}):=\\{\gamma\in\mathbb{R}_{+}^{N\times M}:\gamma
1_{M}\leq 1_{N},\gamma^{T}1_{N}\leq 1_{M}\\},$

and $1_{N}$ denotes the $N\times 1$ vector whose entries are $1$ and
analogously for $1_{M}$.

It has been shown that the optimal plan $\gamma$ for the empirical OPT problem
is induced by a 1-1 mapping [8]. Thus, the above problem can be further
simplified as follows:

$\displaystyle\text{OPT}(\\{\hat{y}^{n}\\}_{n=1}^{N},\\{y^{m}\\}_{m=1}^{M})$
$\displaystyle=\min_{L}\sum_{n\in\text{Dom}(L)}c(\hat{y}^{n},y^{L(n)})$
$\displaystyle+\lambda(N+M-2|\text{Dom}(L)|)$ (6)

where $L:[1:N]\hookrightarrow[1:M]$ is a partial bijection, i.e.
$\text{Dom}(L)\subset[1:N]$, $L$ is a 1-1 mapping.
$|\text{Dom}(L)|=\\#\text{Dom}(L)$ is the cardinality of the set Dom(L). Note
in the registration problem, $L$ denotes a partial correspondence. That is, we
only define the correspondence for points in the domain and range of $L$.

####  III-A3 Primal-form of OPT

If we replace the penalty term of (4) with a constraint, i.e. we impose the
condition $|\gamma|=\zeta$, 222In the original formulation, this constraint is
$|\gamma|\geq\zeta$. It is straightforward to verify the equivalence. Thus in
this paper, we do not distinguish the conditions “$\geq\zeta$” and
“$=\zeta$”., where $\zeta\in[0,\min(|\mu|,|\nu|)]$ is a constant, then (4) is
closely related to the Lagrangian formulation of the following “primal
problem”:

$\displaystyle\text{Primal-
OPT}(\mu,\nu;\zeta)=\inf_{\gamma\in\Pi_{\leq}(\mu,\nu)}$ $\displaystyle\int
c(\hat{y},y)d\gamma(\hat{y},y)$ (7) $\displaystyle\text{s.t. }|\gamma|=\zeta.$

For empirical measures
$\mu=\sum_{n=1}^{N}\delta_{\hat{y}^{n}},\nu=\sum_{m=1}^{M}\delta_{y^{m}}$, the
above problem can be written as

$\displaystyle\text{Primal-
OPT}(\mu,\nu;\zeta)=\min_{\gamma\in\Pi_{\leq}(1_{N},1_{M})}$
$\displaystyle\sum_{n,m}c(\hat{y}^{n},y^{m})\gamma_{n,m}$ (8)
$\displaystyle\text{s.t. }|\gamma|=\zeta.$

####  III-A4 Solvers and their computational complexities

For empirical OT (2), its objective function is linear, and the feasible space
(for $\gamma$) is a polytope. Thus it can be solved by Linear programming
[64]. Furthermore, as shown in [16], OPT can be formulated as a balanced OT
problem by introducing reservoir points, thus empirical OPT (5) and empirical
primal-OPT (8) 333We refer the Appendix B for the computation of primal-OPT.
could also be solved by linear programming. However, the time complexity of
linear programming is $\mathcal{O}(NM(N+M))$.

One idea to facilitate linear programming is to convert the linear
optimization problem into a strictly convex problem. The most successful
approach in this area is Entropy Regularization, which adds the transport
plan’s entropy to the OT objective function and then applies the Sinkhorn-
Knopp algorithm [115, 26]. The algorithm can be extended to the large-scale
stochastic [45] and unbalanced settings [9, 20]. For moderate regularization,
these algorithms converge fast, however, there is a trade-off between accuracy
versus stability and convergence speed for small regularization parameters. In
practice, under an unsuitable regularization, the time complexity of the
Sinkhorn algorithm could, unfortunately, be even worse than linear
programming.

Other notable approaches include restricting the cost function to specific
metrics [48, 88, 106], or enforcing low-rank constraints on $\gamma$ [108,
107]. In particular, when the transportation cost is a metric, network flow
methods [48, 88] can be applied, or when $c$ is a tree metric, an efficient
algorithm based on dynamic programming with time complexity
$\mathcal{O}(n\log^{2}n)$ is proposed in [106]. However, in high dimensions
existence (and identification) of an appropriate metric tree remains
challenging. Lastly, low-rank approximations [108, 107] are shown to lead to
very fast solvers, however, at the cost of accuracy.

####  III-A5 Sliced Optimal (Partial) Transport

When the space is ordered (e.g. 1 Dimensional Euclidean space) and the cost
function is consistent with respect to the order, i.e. $c(x^{1},x^{2})\leq
c(x^{1},x^{3})$ if $x^{1}\leq x^{2}\leq x^{3}$, then by OT’s theory [122], the
balanced OT problem has a closed form solution, i.e., the increasing re-
arrangement function given by the north-west corner rule. Utilizing this
theory, Sliced OT techniques [99, 66, 12, 67, 75] are proposed with the main
idea to calculate the expected OT distance between 1-dimensional marginal
distributions (i.e., slices) of two $d$-dimensional distributions. The
expectation is numerically approximated via a Monte Carlo integration scheme.
Other notable extensions of these distances include the generalized, the max-
sliced, and the convolutional sliced Wasserstein distances [65, 33, 87].

Sliced OT techniques can be extended into the unbalanced transport setting [8,
111]. In particular, Bai et al. [8] define the sliced optimal partial
transport (SOPT) as follows:

###### Definition 1

In $\mathbb{R}^{D}$ space, given $\mu,\nu\in\mathcal{M}_{+}(\mathbb{R}^{D})$
and $\lambda:\mathbb{S}^{D-1}\to\mathbb{R}_{++}$ is an $L_{p}$ function, we
define the sliced optimal partial transport (SOPT) problem as follows:

$\displaystyle\text{SOPT}_{\lambda}(\mu,\nu)=\int_{\mathbb{S}^{D-1}}\text{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu,\langle\theta,\cdot\rangle_{\\#}\nu)d\sigma(\theta)$
(9)

where $\text{OPT}_{\lambda}(\cdot,\cdot)$ is defined in (4),
$\sigma\in\mathcal{P}(\mathbb{S}^{D-1})$ is a probability measure such that
$\text{supp}(\sigma)=\mathbb{S}^{D-1}$.

The 1D OPT problem
$\text{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu,\langle\theta,\cdot\rangle_{\\#}\nu)$
can be solved efficiently by the algorithm 1 in [8]. Similar to OPT problem,
when cost function $c$ is a ($p$-power of a) metric, SOPT (9) defines a metric
in $\mathcal{M}_{+}(\mathbb{R}^{D})$.

###  III-B Procrustes Alignment Problem

Procrustes analysis [109, 27] is a well-known method to learn a linear
transformation between two sets of matched points (point clouds)
$X\in\mathbb{R}^{N\times D}$ and $Y\in\mathbb{R}^{N\times D}$. If the
correspondences between the two sets are known (i.e., which point of $X$
corresponds to which point of $Y$), and for convenience, suppose the
correspondence is induced by identity matrix (i.e. $x^{n}\mapsto y^{n},\forall
n$), then this linear transformation can be recovered by solving the least
square problem: ${\rm min}_{W\in\mathbb{R}^{D\times D}}\|XW-Y\|_{2}^{2}$. The
case ${\rm min}_{W\in O_{D}}\|XW-Y\|_{2}^{2}$ where $O_{D}$ is the set of
orthogonal matrices has a closed form solution discovered by Schönemann [109]
given by $W^{*}=UV^{T}$ where $USV^{T}$ is the singular value decomposition of
$X^{T}Y$. The orthogonality ensures that the distances between points are
unchanged by the transformation. This problem is called the “Orthogonal
Procrustes Problem.”

Procrustes alignment problem assumes known correspondences between source and
target points. When these correspondences are not known, the optimal transport
(or optimal partial transport) can be used to first obtain correspondences (or
soft correspondences) between the source and target and then solve for the
linear transformation $W$. This problem is referred to as “Wasserstein
Procrustes” [47]. Procrustes Wasserstein is formalized as:

${\rm min}_{W\in O_{D}}{\rm min}_{P\in{\cal P}_{N}}\|XW-PY\|_{2}^{2}$

where ${\cal P}_{N}$ is the set of $N\times N$ permutation matrices. While for
a fixed $P$ the problem is convex in $W$, and vice versa, the problem is,
unfortunately, not jointly convex in $W$ and $P$. Hence, stochastic
alternating optimization techniques are utilized to solve this optimization
problem [47]. The rigid registration problem between source point cloud $X$
and target point cloud $Y$ can be regarded as the above Wasserstein Procrustes
problem, and the transformation model is defined as $x\mapsto W^{T}x.$

##  IV Our Method

In space $\mathbb{R}^{D}$, we consider the registration map,
$f:\mathbb{R}^{D}\to\mathbb{R}^{D}$, as follows:

$\displaystyle f(x):=\sum_{k=1}^{K}\alpha^{k}\phi(x,c^{k})+B^{T}x+\beta$ (10)

where $\alpha^{k}\in\mathbb{R}^{D}$ for each $k$,
$\phi:(\mathbb{R}^{D})^{2}\to\mathbb{R}$ is a smooth kernel function 444We
refer Appendix A for the introduction of kernel function. The TPS kernel will
also be discussed in subsection B., $B\in\mathbb{R}^{D\times D}$ modeling the
linear portion of the deformation, $\beta\in\mathbb{R}^{D}$ is the translation
parameter, $c^{1},...,c^{k}\in\mathbb{R}^{D}$ are called control points. In
the remainder of the paper, we assume that the linear part of the registration
map (10) is restricted to scaling and rotation $B=SR$, where
$S=\text{diag}(s_{1},\ldots,s_{D})$ with each $s_{i}>0$ is the scaling matrix,
and $R\in\mathbb{R}^{D\times D}$ is the rotation matrix. Hence, we use the
following registration map 555In methods based on TPS (Thin-Plate Spline), as
discussed in sub-sections B and D, we continue to apply the traditional TPS-
regression technique. This approach typically recovers a general linear matrix
$B$. As such, these methods do not guarantee a linear mapping in the form of
$SR$ finally, unless the true deformation mapping between $X$ and $Y$
satisfies this assumption.:

$\displaystyle f(x):=\sum_{k=1}^{K}\alpha^{k}\phi(x,c^{k})+R^{T}Sx+\beta.$
(11)

At their core, the proposed methods iteratively update the transformation
parameters in two steps, estimating the correspondence and updating the
transform parameters based on recent correspondence. For the first step, we
will apply sliced optimal partial transport and optimal partial transport; for
the second step, we will apply the RBF regression method and TPS regression
methods. In particular, our methods will be introduced in the following four
subsections. Additionally, as one condition in the problem setup, we assume
$\zeta$ is the number of points of the “clean part,” which is given as our
prior knowledge. We first introduce the main ideas of our methods in the
following four subsections; then, in the last section, we summarize these
methods as pseudo-codes before we proceed to our numerical experiments.

Before we introduce our methods, we first formulate the point cloud
registration with noise as the following problem:

$\displaystyle\min_{f}\min_{L}\sum_{n\in\text{Dom}(L)}\|f(x^{n})-y^{L[n]}\|^{2}+\epsilon\text{reg}(f)$
(12)

where $L$ is over all $[1:N]\hookrightarrow[1:M]$ partial bijections such that
$|\text{Dom}(L)|=\zeta$. And there exists an underlying solution $L^{*}$, that
satisfies the above conditions. The $\text{reg}({f})$ term is a regularization
term. When $f$ is RBF model (11) or TPS model (40), we refer (16) and (55) for
its particular formulation.

###  IV-A OPT-RBF

Our paper’s primary contribution is the application of sliced OPT to the non-
rigid registration of point sets. To provide context and motivation for our
approach, we begin by formulating the problem using the primal form of the
partial transport, as given in Eq. (8).

Given two point sets $X=\\{x^{n}\in\mathbb{R}^{D}\\}_{n=1}^{N}$ and
$Y=\\{y^{m}\in\mathbb{R}^{D}\\}_{m=1}^{M}$, we treat these sets as empirical
measures, $\mu=\sum_{n=1}^{N}\delta_{\hat{y}^{n}}$ and
$\nu=\sum_{m=1}^{M}\delta_{y^{m}}$. The OPT-RBF then can be formalized as the
optimization problem in (12). Note that this optimization problem is similar
to the Wasserstein Procrustes problem [47] but with two significant
differences: 1) $f$ is non-rigid, and 2) we use OPT as opposed to OT. To solve
this optimization problem, we propose an alternating iterative optimization
scheme through the following two steps:

Step 1. For fixed registration parameters, first calculate,

$\hat{y}^{n}=f(x^{n}),\forall n\in[1:N],$

where $f$ is defined by (11). Regarding control point set
$\\{c^{1},\ldots,c^{K}\\}$, we can either set it as source point set
$\\{x^{1},\ldots x^{n}\\}$ or choose a different point set based on prior
model knowledge. Then we solve the primal OPT problem:

$\displaystyle\gamma^{*}=\text{primal-
OPT}(\\{\hat{y}^{n}\\}_{n=1}^{N},\\{y^{n}\\}_{m=1}^{M};\zeta)$ (13)

where $\zeta$ is the number of points (the mass) of the clean data point
cloud, and $\gamma^{*}$ is the optimal partial transportation plan. Let

$\displaystyle\mathcal{D}=\\{n:\sum_{m=1}^{M}(\gamma_{n,m})>0\\},$ (14)

then for each $n\in\mathcal{D}$, we update $\hat{y}^{n}$ via:

$\displaystyle\hat{y}^{n}\leftarrow\left(\frac{1}{\sum_{m=1}^{M}(\gamma_{n,m})}\gamma
Y\right)[n,:].$ (15)

which is the barycentric projection of the partial transport plan [5, 7].
Lastly, we obtain the following correspondence666See appendix A for detailed
explanation of the term “correspondence”.:

$\\{(x^{n},\hat{y}^{n})\\}_{n\in\mathcal{D}}.$

Step 2. Given a correspondence $\\{(x^{n},\hat{y}^{n})\\}_{n\in\mathcal{D}}$,
we will update $\\{\alpha^{k}\\}_{k=1}^{K},S,R,\beta$. For convenience, by
reindexing, we suppose correspondence is
$\\{(x^{n},\hat{y}^{n})\\}_{n=1}^{N_{sub}}$, where $N_{sub}\leq\min(N,M)$,
i.e. $\mathcal{D}=[1:N_{sub}]$.

The parameters $S,R,\beta,\alpha$ are selected via to minimizing the
following:

$\displaystyle\min_{S,R,\beta,\alpha}\sum_{n=1}^{N_{sub}}\left\|\left(\sum_{k=1}^{K}\alpha^{k}\phi(x^{n},c^{k})+R^{T}Sx^{n}+\beta\right)-\hat{y}^{n}\right\|^{2}$
$\displaystyle\quad\quad\quad+\epsilon\sum_{k=1}^{K}\|\alpha_{k}\|^{2}$ (16)

where $\epsilon>0$, and the regularization term
$\epsilon\sum_{i=1}^{K}\|\alpha_{k}\|^{2}$ is applied to restrict the model
flexibility and improve the numerical computational stability.

Let

$\displaystyle X_{sub}=\left[\begin{array}[]{c}(x^{1})^{T}\\\ \ldots\\\
(x^{N_{sub}})^{T}\end{array}\right],\hat{Y}_{sub}=\left[\begin{array}[]{c}(\hat{y}^{1})^{T}\\\
\ldots\\\
(\hat{y}^{N_{sub}})^{T}\end{array}\right],\alpha=\left[\begin{array}[]{c}(\alpha^{1})^{T}\\\
\ldots\\\ (\alpha^{K})^{T}\end{array}\right],$ (26)
$\displaystyle\Phi=\left[\begin{array}[]{ccc}\phi(x^{1},c^{1})&\ldots&\phi(x^{1},c^{K})\\\
\ldots&\\\
\phi(x^{N},c^{1})&\ldots&\phi(x^{N},c^{K})\end{array}\right],\Phi_{sub}=\Phi[\mathcal{D},:].$
(30)

Thus, (16) can be written in the following matrix form:

$\displaystyle\min_{S,R,\beta,\alpha}\|\Phi_{sub}\alpha+X_{sub}SR+\beta^{T}1_{N_{sub}}-\hat{Y}_{sub}\|^{2}+\epsilon\cdot\text{tr}(\alpha^{T}\alpha)$
(31)

We will update the transformation parameters in two steps:

Step 2.1: Fix $\alpha$, and let
$\hat{Y}^{\prime}=\hat{Y}_{sub}-\Phi_{sub}\alpha$, that is, each entry, noted
as $\hat{y}^{\prime n}=\hat{Y}^{\prime}[n,:]^{T}$, is defined as

$\hat{y}^{\prime
n}=\hat{y}^{n}-\sum_{k=1}^{K}\alpha^{k}\phi(x^{n},c^{k}),\forall
n\in[1:N_{sub}].$

We aim to solve

$\min_{S,R,\beta}\sum_{n=1}^{N_{sub}}\|f(x^{n})-\hat{y}^{n}\|^{2}=\|X_{sub}SR+\beta^{T}1_{N_{sub}}-Y^{\prime}\|^{2}.$

By [36], we can obtain the optimal $R,S,\beta$ as follows:

Let $Y_{c}=\hat{Y}^{\prime}-\left[\begin{array}[]{c}\bar{y}^{T}\\\ \ldots\\\
\bar{y}^{T}\end{array}\right]$ where
$\bar{y}=\frac{1}{N_{sub}}\sum_{n=1}^{N_{sub}}(\hat{y}^{\prime n})$ and
$y_{c}^{n}=Y_{c}[n,:]^{T}$ denote $n^{th}$ entry of $Y_{c}$, for each $n$.
Similarly, we define $X_{c}$ and $x_{c}^{n}$.

Given a scaling matrix $S$ (initially, $S$ is set to be identity matrix
$I_{D}$), let

$\displaystyle H_{c}=(X_{c}S)^{T}Y_{c}=U_{H}\Sigma_{H}V_{H}^{T}.$ (32)

The optimal $R$ (for fixed $S$) is given by

$\displaystyle R=U_{H}\left[\begin{array}[]{cc}I_{D-1}&0\\\
0&\text{det}(V_{H}U_{H}^{T})\end{array}\right]V_{H}^{T}$ (35)

Next, given $R$, the optimal $S=\text{diag}(s_{1},\ldots,s_{d})$ is computed
by

$\displaystyle
s_{d}=\frac{\sum_{i=1}^{N_{sub}}y^{iT}_{c}R^{T}E_{d}x^{i}_{c}}{\sum_{i=1}^{N_{sub}}x^{iT}_{c}E_{d}x^{i}_{c}},\forall
d$ (36)

where $E_{d}=[0,\ldots,e_{d},\ldots,0]$ and $e_{1},...e_{D}$ are the Canonical
vectors in $\mathbb{R}^{D}$. We iterate through steps (35) and (36) until $R$
and $S$ converge. Note, in the case of uniform scaling, where $S=sI_{D}$ for
some $s>0$, it is sufficient to compute (35) just once to find the optimal
values for $R$. Additionally, the scaling factor $S=sI_{D}$ can be calculated
in two ways. One approach is to use the average value of $\\{s_{1},\ldots
s_{d}\\}$, as specified by Equation (36). Alternatively, $s$ can be determined
using the formula: $s=\sqrt{\frac{\text{cov}(Y_{c})}{\text{cov}(X_{c})}}.$

Finally, the optimal translation $\beta$ is computed by:

$\displaystyle\beta$
$\displaystyle=\frac{1}{N_{sub}}\sum_{n=1}^{N_{sub}}[(\hat{y}^{\prime
n}-R^{T}Sx^{n})]$
$\displaystyle=\frac{1}{N_{sub}}[1_{N}^{T}\hat{Y}^{\prime}-1_{N}^{T}X_{sub}SR]$
(37)

Step 2.2: Given $S,R,\beta$, we find the optimal $\alpha$. Let
$\hat{Y}^{\prime\prime}=Y-X_{sub}SR-\beta^{T}1_{N}$, we have

$\displaystyle\arg\min_{\alpha}\|\Phi_{sub}\alpha+X_{sub}SR+\beta^{T}1_{N}-\hat{Y}_{sub}\|^{2}+\epsilon\text{tr}(\alpha^{T}\alpha)$
$\displaystyle=\arg\min_{\alpha}\|\Phi_{sub}\alpha-\hat{Y}^{\prime\prime}\|^{2}+\epsilon\text{tr}(\alpha^{T}\alpha)$
$\displaystyle=\arg\min_{\alpha}\sum_{d=1}^{D}\|\Phi_{sub}\alpha[:,d]-\hat{Y}^{\prime\prime}[:,d]\|^{2}+\epsilon\alpha[:,d]^{T}\alpha[:,d]$

It suffices to solve the following 1D regression problem:

$\displaystyle\min_{\alpha[:,d]}\|\Phi_{sub}\alpha[:,d]-\hat{Y}^{\prime\prime}[:,d]\|^{2}+\epsilon\alpha[:,d]^{T}\alpha[:,d]$
(38)

for each $d\in[1:D]$, whose solution is given by

$\displaystyle\alpha[:,d]=(\Phi^{T}_{sub}\Phi_{sub}+\epsilon
I_{K})^{-1}\Phi^{T}_{sub}\hat{Y}^{\prime\prime}[:d]$

Therefore, optimal $\alpha$ is given by

$\displaystyle\alpha=(\Phi^{T}_{sub}\Phi_{sub}+\epsilon
I_{K})^{-1}\Phi^{T}_{sub}\hat{Y}^{\prime\prime}.$ (39)

We repeat these two steps until convergence.

###  IV-B OPT-TPS

The second method we introduced is the combination of OPT and TPS, it can be
regarded as a variant of TPS-RPM 777For further details, please see Appendix
C.. In particular, we use the TPS model:

$f(x)=\sum_{n=1}^{N}\alpha^{n}\phi_{T}(x,x^{n})+B^{T}x+\beta$ (40)

Same to the Step 1 in OPT-RBF, we use primal-OPT to estimate the
correspondence and obtain $\\{(x^{n},\hat{y}^{n})\\}_{n=1}^{N}$. It is
important to note that in this method, we utilize the complete set of $x$ and
$\hat{y}$ points.

Next, we discuss how to update the parameters of the above model based on a
given correspondence. We first review the TPS regression technique.

Introduction of TPS regression in $\mathbb{R}^{D}$.

As we introduced in Section II, the thin-plate spline is derived from the
following minimization problem,

$\inf_{\tilde{f}:\mathbb{R}^{D}\to\mathbb{R}}E[\tilde{f}]:=\inf_{\tilde{f}:\mathbb{R}^{D}\to\mathbb{R}}\sum_{n=1}^{N}\|\tilde{f}(x^{n})-\tilde{y}^{n}\|^{2}+\epsilon\int_{\mathbb{R}^{D}}\|\nabla^{2}\tilde{f}\|^{2}dx.$
(41)

where
$\|\nabla^{2}\tilde{f}\|^{2}=\sum_{i,j=1}^{D}\left(\frac{\partial^{2}\tilde{f}}{\partial
x_{i}\partial x_{j}}\right)^{2}$, the regularization term is called bending
energy; $\\{x^{n}\\}_{n=1}^{N}\subset\mathbb{R}^{D}$ and
$\\{\tilde{y}^{n}\\}_{n=1}^{N}\in\mathbb{R}$ are given fixed data points.

It has been proved that the solution $\tilde{f}$ for above (41) has the
following closed form:

$\displaystyle\tilde{f}(x)$
$\displaystyle=\sum_{n=1}^{N}\alpha^{n}\phi_{T}(x,x^{n})+b^{T}x+\tilde{\beta}$
(42)

where $\tilde{\beta},\alpha^{n}\in\mathbb{R},b\in\mathbb{R}^{n}$ and

$\displaystyle\phi_{T}(x,x^{n}):=\begin{cases}\frac{1}{2}\|x-x^{n}\|^{2}\ln(\|x-x^{n}\|^{2})&\text{if
}D=2\\\ \ln(\|x-x_{0}\|)&\text{if }D=4\\\
\|x-x_{0}\|^{4-D}&\text{otherwise}\end{cases}.$ (43)

###### Remark 1

In practice, when $D\geq 4$, $\phi_{T}$ has singularity at $x=x^{n}$ and thus
$\phi_{T}(x,x^{n})$ is set to be $\|x-x^{n}\|^{2}$ or
$\|x-x^{n}\|^{2}\ln(\|x-x^{n}\|^{2})$. Therefore the constructed $f(x)$ is
generally not the minimizer for (41).

Let

$\displaystyle\bar{X}=\left[\begin{array}[]{c}1,x^{1}\\\ \cdots\\\
1,x^{N}\end{array}\right]=[1_{N},X],\tilde{Y}=\begin{bmatrix}\tilde{y}_{1}\\\
\ldots\\\
\tilde{y}^{N}\end{bmatrix},\bar{b}=\left[\begin{array}[]{c}\tilde{\beta}\\\
b\end{array}\right].$ (49)

The optimal parameter $\alpha,\bar{b}$ can be constructed by solving the
following linear system

$\displaystyle\begin{cases}&\tilde{Y}=(\Phi+\epsilon
I_{N})\alpha+\bar{X}\bar{b}\\\ &\bar{X}^{T}\alpha=0_{N}\end{cases}$ (50)

where the second equation follows from the functional analysis in [82], which
force the $\alpha$ be in the null space of $\bar{X}$. Problem (50) can be
solved by the following linear system

$\begin{bmatrix}\Phi+\epsilon I_{N}&\bar{X}\\\
\bar{X}^{T}&0\end{bmatrix}\begin{bmatrix}\alpha\\\
\bar{b}\end{bmatrix}=\begin{bmatrix}\tilde{Y}\\\ 0\end{bmatrix}$ (51)

Thus

$\displaystyle\begin{bmatrix}\alpha\\\
\bar{b}\end{bmatrix}=\begin{bmatrix}\Phi+\epsilon I_{N}&\bar{X}\\\
\bar{X}^{T}&0\end{bmatrix}^{-1}\begin{bmatrix}\tilde{Y}\\\ 0\end{bmatrix}.$
(52)

Alternatively, the matrix decomposition technique can be applied, and thus we
have optimal parameter $(\alpha,\bar{b})$ has the following closed form:

$\displaystyle\begin{cases}&\alpha=(\Phi+\epsilon
I_{N})^{-1}(\tilde{Y}-\bar{X}\bar{b}),\\\ &\bar{b}=(\bar{X}^{T}(\Phi+\epsilon
I_{N})^{-1}\bar{X})^{-1}\bar{X}^{T}(\Phi+\epsilon
I_{N})^{-1}\tilde{Y}.\end{cases}$ (53)

However, both (52)(53) are not computational practicable and [123] proposed an
improved computational method:

We first write the QR decomposition of $\bar{X}$:

$\bar{X}=[Q_{1},Q_{2}]\begin{bmatrix}\mathcal{R}\\\ 0\end{bmatrix},$

where $Q_{1},Q_{2}$ are $(D+1)\times N,(N-D-1)\times N$ orthogonal matrices;
and $R\in\mathbb{R}^{(D+1)\times(D+1)}$ is upper triangular.

$\displaystyle\begin{cases}&\alpha=Q_{2}(Q_{2}^{T}(\Phi+\epsilon
I_{N})Q_{2})^{-1}Q_{2}^{T}\tilde{Y}\\\
&\bar{b}=\mathcal{R}^{-1}Q_{1}^{T}(\tilde{Y}-(\Phi+\epsilon
I_{N})\alpha)\end{cases}.$ (54)

Step 2: TPS method regression for point cloud registration.

Based on the above TPS interpolation technique, given a correspondence
$(\\{(x^{n},\hat{y}^{n})\\}_{n=1}^{N})$, the $D-$dimensional TPS regression
can be simplified to be the following:

$\displaystyle\inf_{\alpha,B,\beta}\sum_{n=1}^{N}\|f(x^{n})-\hat{y}^{n}\|^{2}+\epsilon\sum_{d=1}^{D}\int_{\mathbb{R}^{D}}\|\nabla^{2}f[d]\|^{2}dx.$
(55)

The optimal parameters $\alpha,B,\beta$ can be computed by solving the 1D TPS
interpolation problem

$\inf_{\alpha[:,d],B[:,d],\beta[d]}\sum_{n=1}^{N}\|f(x^{n})[d]-\hat{y}^{n}[d]\|^{2}+\epsilon\int_{\mathbb{R}^{D}}\|\nabla^{2}f[d]\|^{2}dx$

for each $d\in[1:D]$. By (52) (53), let

$\bar{B}=\begin{bmatrix}\beta^{T}\\\
B\end{bmatrix},\hat{Y}=\begin{bmatrix}\hat{y}^{1T}\\\ \ldots\\\
\hat{y}^{NT}\end{bmatrix},$

we have

$\displaystyle\begin{cases}&\alpha=Q_{2}(Q_{2}^{T}(\Phi+\epsilon
I_{N})Q_{2})^{-1}Q_{2}^{T}\hat{Y}\\\
&\bar{B}=\mathcal{R}^{-1}Q_{1}^{T}(\hat{Y}-(\Phi+\epsilon
I_{N})\alpha).\end{cases}$ (56)

###  IV-C SOPT-RBF methods

In this method, we use the RBF modle (11). As we discussed in III, computing
OPT can be time consuming. In improve the computational efficiency. We use
Sliced-OPT to replace the primal-OPT in step 1 in method OPT-RBF.

Step 1. We first select
$\\{\theta_{1},\theta_{2},\ldots\theta_{t}\\}\subset\mathbb{S}^{D-1}$. For
$t^{\prime}=1,2,\ldots,t$, we project $Y$ and $\hat{Y}$ into 1D space spanned
by $\theta_{t^{\prime}}$ and obtain

$\hat{Y}_{\theta_{t^{\prime}}}:=\\{\theta_{t^{\prime}}^{T}\hat{y}^{n}\\}_{n=1}^{N},Y_{\theta_{t^{\prime}}}=\\{\theta_{t^{\prime}}^{T}y^{m}\\}_{m=1}^{M},$

and then we solve the 1D OPT problem

$\displaystyle\text{OPT}_{\lambda}(\hat{Y}_{\theta_{t^{\prime}}},Y_{\theta_{t^{\prime}}}).$
(57)

Let $L_{t^{\prime}}$ denote the OPT transportation map,888In OPT-RBF method,
we use transportation matrix $\gamma$ to represent the transportation plan.
While we use Monge mapping $L$ to denote the transportation plan in this
section. These two descriptions are equivalent. For convenience, we use these
two descriptions interchangeably. . then we update $\hat{y}^{n}$ with:

$\displaystyle\hat{y}^{n}\leftarrow\hat{y}^{n}+\theta_{t^{\prime}}(\theta_{t^{\prime}}^{T}y^{L_{t^{\prime}}(n)}-\theta_{t^{\prime}}^{T}\hat{y}^{n}).$
(58)

We repeat the above process for $\theta_{1},\theta_{2},\ldots\theta_{t}$.

Let $\mathcal{D}_{t^{\prime}}:=\text{Dom}(L_{t^{\prime}})$, and
$\mathcal{D}=\bigcup_{t^{\prime}}\mathcal{D}_{t^{\prime}}$, that is,
$\mathcal{D}$ is union of index of $\hat{y}^{n}$ which has been moved in above
process. We have $\\{({x}^{n},\hat{y}^{n})\\}_{n\in\mathcal{D}}$ defining a
correspondence. For step 2, it is same to the step 2 in OPT-RBF.

###  IV-D SOPT-TPS

This section is an SOPT version of the OPT-TPS method. In particular, we use
TPS model (40).

In particular, we apply the Step 1 in SOPT-RBF to estiamte the correspondence
and use Step 2 in OPT-TPS to update the parameter.

###  IV-E Summary and pseudo-code

Input: $X,Y,\zeta,T,\phi,\text{param},\epsilon,C\leftarrow X$

Output: $R,S,\beta,\Phi,\alpha$

1 Initilize $R,S,\beta,\alpha$

2 Initilize $\Phi\leftarrow\phi(X,C^{T},\text{param})$

3 for  _$T^{\prime}=1,2,\ldots T$_ do

/* Step 1 */

4 $\hat{Y}\leftarrow\Phi\alpha+XSR+\beta^{T}1_{N}$

5 compute optimal plan $\gamma$ for Primal-OPT$(\hat{Y},Y;\zeta)$.

6 $\mathcal{D}\leftarrow\\{n:\sum_{m=1}^{M}\gamma_{n,m}>0\\}$

7 update $\hat{y}^{n},\forall n\in\mathcal{D}$ via (15)

8
$(X_{sub},\hat{Y}_{sub},\Phi_{sub})\leftarrow(X[\mathcal{D},:],\hat{Y}[\mathcal{D},:],\Phi[\mathcal{D},:])$

/* Step 2 */

9 $\hat{Y}^{\prime}\leftarrow\hat{Y}_{sub}-\Phi_{sub}\alpha$

10 By $(X_{sub},\hat{Y}^{\prime})$, update $R,S,\beta$ via (35),(36),(37).

11

12 if  _condition for non-rigid is True_ then

13
$\hat{Y}^{\prime\prime}\leftarrow\hat{Y}_{sub}-(X_{sub}SR+\beta^{T}1_{N_{sub}})$

14 By $(\Phi_{sub},\hat{Y}^{\prime\prime})$, update $\alpha$ via (39)

15

Algorithm 1 OPT-RBF

The pseudo-code for OPT-RBF is presented in Algorithm 1. For the inputs, $X,Y$
are source and target point clouds, $\zeta$ is the number of clean data points
of $X$, which is the prior knowledge and has been elaborated upon in Section
IV. The variable ’param’ represents the kernel-specific parameters; for
example, $\sigma^{2}$ in the case of the Gaussian kernel and dimension $D$ for
the TPS (Thin-Plate Spline) kernel. $C=[c^{1},\ldots
c^{K}]^{T}\in\mathbb{R}^{K\times D}$ is the set of control points, where
$K\in\mathbb{N}$. By default, $C$ is set to $X$, but it can be configured to
other sets of points represented as matrices.

Initialization occurs at line 1, where we conventionally set $R=I_{D}$,
$S=I_{D}$, $\beta=0_{D}$, and $\alpha=0_{N\times D}$. Line 11 sets forth a
condition to initiate the nonrigid registration process. Typically, this
condition is established as enough rigid iterations, coupled with the
convergence of the linear parameters, $R,S,\beta$, or $\bar{B}$.

For the pseudo-code of the other three algorithms, please refer to Appendix E.

###### Remark 2

In a balanced sliced OT setting, the steps specified by (58) is formally known
as sliced Wasserstein gradient flow [42, 75]. As shown in [73, Theorem 1.1,
Theorem 4.7], the transported point set $\hat{Y}$ converges to $Y$ as the
number of iterations tends to infinity. However, convergence properties in the
partially sliced OT setting are not yet fully understood. Empirical
observations indicate rapid convergence when the rotation angles between $X$
and $Y$ are smaller than $\pi/2$. A comprehensive theoretical analysis of this
phenomenon is a subject for future research.

###### Remark 3

As we discussed in sub-section A, under the balanced OT setting, if we exclude
the non-rigid part “$\Phi\alpha$”, the problem formulation given by (12) is
known as the Wasserstein Procrustes problem. As per [47], the gradient
technique with full sample points faces challenges related to statistical
convergence and computational expense. The computational expense of solving
OT(or OPT) has been discussed in section III. Regarding convergence, the
Wasserstein Procrustes is a non-convex problem.999To navigate this, one should
reinterpret $L$, equivalent to a permutation matrix, as any bi-stochastic
matrix. This modification makes the search space for the variable $L$ convex.
Thus, the discussion of convexity becomes suitable. In [47], the benefits of
stochastic gradient descent are explored. sliced OT or (OPT) can intuitively
produce similar results. By reducing $X,Y$ to a one-dimensional domain,
partial information is lost. As a result, the correspondence derived is
generally different from the correspondence from the original D-dimensional OT
(or OPT) approach, and randomness has been introduced. Integrating with the
partial OT setting, only the points in $\mathcal{D}$ are transported, which is
similar to the sub-sampling step in gradient descent. This intuition aligns
with our experiment observation, and we aim to rigorously evaluate the
advantages of the sliced unbalanced OT method concerning the Wasserstein
Procrustes problem in future works.

##  V Experiments

Wall clock time (in second) and registration error for all methods
female male neutral registration error $5\%$ $10\%$ $5\%$ $10\%$ $5\%$ $10\%$
$5\%$ $10\%$ baselines CPD [85] $5.44$ $5.86$ $5.37$ $6.82$ $5.35$ $6.74$
0.063(0.022) 0.362(0.424) OT-RBF $10.76$ $11.70$ $10.13$ $11.90$ $9.50$
$10.33$ 0.339(0.112) 0.550(0.040) OT-TPS $13.81$ $14.01$ $13.26$ $14.33$
$12.35$ $14.21$ 0.076(0.012) 0.151(0.044) SOT-RBF $3.47$ $3.50$ $3.43$ $3.49$
$3.36$ $3.49$ 0.121(0.062) 0.219(0.185) SOT-TPS $6.38$ $7.52$ $6.43$ $7.53$
$6.41$ $7.58$ 0.046(0.007) 0.065(0.004) TPS-RPM[23] $15.1$ $16.6$ $15.2$
$17.6$ $14.8$ $16.4$ 0.183(0.011) 0.218(0.018) TPS-RPM(new) $14.5$ $16.7$
$15.5$ $17.7$ $15.3$ $17.8$ 0.039(0.010) 0.041(0.010) ours OPT-RBF $11.43$
$79.1$ $11.53$ $72.9$ $11.11$ $73.8$ $0.037(0.008)$ 0.040(0.011) OPT-TPS
$14.91$ $74.6$ $15.45$ $72.4$ $15.40$ $64.4$ 0.033(0.010) 0.033(0.010) SOPT-
RBF $8.37$ $9.90$ $8.42$ $9.55$ $9.24$ $9.93$ 0.036(0.013) $0.042(0.011)$
SOPT-TPS $11.71$ $12.21$ $11.59$ $12.72$ $11.55$ $12.31$ 0.034(0.009)
0.040(0.009)

TABLE I: In this table, the first six columns display the wall clock time for
each method, measured in seconds per iteration. All methods in our experiment
converge within a range of 45 to 65 iterations. The final two columns show the
registration error for each method. In these cells, the first value represents
the mean registration error, where the error is defined in 59, while the value
in parentheses indicates the standard deviation.

In this experiment, the methods being compared are our methods (see section
IV): SOPT-RBF, SOPT-TPS, OPT-RBF, OPT-TPS; and baseline methods: CPD [85],
TPS-RPM [23, 129],and an improved version of TPS-RPM (See Appendix C for more
details), denoted as TPS-RPM(new). In addition, we extend the OT Procrustes
method[47] into the non-rigid setting by our RBF and TPS models, and denoted
as OT-RBF, OT-TPS. Furthermore, we provide the sliced-OT version, denoted as
SOT-RBF and SOT-TPS. These five methods are baselines for the experiments.

The dataset being used is STAR dataset (https://github.com/ahmedosman/STAR).
This dataset contains a list of point clouds with labels “female”, “male”, and
“neutral”. For each label, we randomly select two datasets, denoted as $X$ and
$Y$ as the source and target datasets, respectively.

The support of all point clouds is $[-1,1]^{3}$. For each pair of source and
target point clouds, the rotation angle in each dimension is set to be in the
range $[-\frac{\pi}{3},\frac{\pi}{3}]$. Additionally, the translation $\beta$
is selected such that$|\beta[d]|\leq 1$ for $d\in[1:3]$. Scaling is set to be
the identity matrix and remains unchanged for all methods.

Both the source and target datasets have uniform noise added to them, which is
distributed in the support $[-1,1]^{3}$. This experiment compares the accuracy
and running time performance of the mentioned methods under these conditions.

Accuracy. Suppose $X,Y$ are the source and target data respectively and
$X_{0}\subset X,Y_{0}\subset Y$ be the clean part of $X,Y$. Say $L^{*},f^{*}$
is the ground true correspondence and deformation function, i.e.

$y^{L^{*}(n)}=f^{*}(x^{n}),\forall x^{n}\in X_{0}.$

Suppose $\zeta$ is the size of $Y_{0}$. Let $\sigma(Y_{0})$ be the standard
deviation of $Y_{0}$, i.e. $\sigma(Y_{0})>0$ such that
$\frac{1}{D}\frac{1}{N}\left\|\frac{Y_{0}-\text{mean}(Y_{0})}{\sigma(Y_{0})}\right\|^{2}=1.$
For each method, let function $f$ be the returned model, the approximation
error is defined as

$\displaystyle\text{error}=\left(\frac{1}{\zeta}\sum_{x^{n}\in
X_{0}}\left\|\frac{y^{L^{*}(n)}-f(x^{n})}{\sigma(Y_{0})}\right\|^{2}\right)^{1/2}$
(59)

For each method, parameters are carefully selected for optimal performance.
Please refer to the following section on Performance Analysis for a detailed
introduction. Table I and Figure 2 reveal that our methods, along with TPS-
RPM(new), consistently yield superior accuracy, thanks in part to the
utilization of prior knowledge $\zeta$. Additionally, sliced-OT methods, SOT-
RBF and SOT-TPS, also contribute to improved accuracy. This aligns with our
intuitive understanding presented in Remark 3, which suggests that the sliced-
OT approach introduces randomness in the corresponding estimation step. This
proves advantageous in the registration process, particularly when the
correspondence established by balanced OT is inaccurate due to random noise.

Figure 2: This figure visually represents the final result obtained from all
methods, all under a noise level of 10%.

Performance analysis. We evaluated the wall-clock time for our methods
alongside the four baseline techniques, the results of which can be found in
Table (I). This table displays the wall-clock time per iteration and the total
number of iterations needed for convergence for each method and dataset.

All the methods are configured to run for a maximum of 100 iterations, and
rigid registrations are performed exclusively within the initial 20.
Convergence is generally observed between 45 and 65 iterations.

In the CPD method, the Gaussian parameter term $\sigma^{2}$ is fixed at 16.0.
The scaling term is set to 1.0, and the weight term used for the non-rigid
M-step in linear regression is fixed at 4.0.

In TPS-RPM and TPS-RPM (new), the Sinkhorn algorithm [26] is accelerated by
using Numba (https://numba.pydata.org/). We’ve set the weight of entropy
regularization to $\sigma(Y_{0})*0.01$ and capped the maximum number of
iterations at 500 to prioritize the performance. As for OPT-RBF and OPT-TPS,
we utilize the OT solver from PythonOT [44], a C++ library, and impose a
maximum of 1e7 iterations. Regarding SOPT-RBF and SOPT-TPS, we set the number
of projections to be 100 and the 1D-OPT solver [8] is accelerated via Numba
https://numba.pydata.org/. For all these methods, steps involving the
computation of optimal non-rigid parameters, particularly those requiring the
inverse of a large-size matrix, are executed on a GPU. The data type of each
dataset is the 64-bit float number, and all the experiments are conducted on a
machine with an AMD EPYC 7713 64-core Processor and four NVIDIA RTX A6000
GPUs. The CPU operates at a maximum of 3720.7 MHz and supports multi-threading
with 128 threads.

##  VI Summary

In this paper, we proposed new non-rigid registration methods for partial
matching scenarios by incorporating optimal partial transport and classical
non-rigid registration models, RBF and TPS. In addition, to improve the
computation efficiency, we also propose sliced-OPT-based methods. We
demonstrate our methods in 3D and 2D datasets. In our experiment, when data is
corrupted by some proportion of noise data, the partial OT-based method
induces better accuracy. For future research, we plan to investigate the
convergence behavior of methods based on OPT and sliced OPT. Moreover, we aim
to explore the potential benefits of using a sliced-OT-based approach in
addressing Wasserstein Procrustes problems, as suggested by our findings in
both 3D and 2D experiments. Additionally, we intend to examine the potential
applications of generalized sliced unbalanced OT to this problem.

## Acknowledgement

SK acknowledges partial support from the Defense Advanced Research Projects
Agency (DARPA) under Contract No. HR00112190135 and HR00112090023, and the
Wellcome LEAP Foundation. SBD thanks the American Mathematical Society and
Princeton University for financial support. He also thanks Charles Fefferman
and Keaton Hamm for many discussions related to ‘slow twists’ and ‘slides.’

## References

  * [1] Hossam Abdelmunim and Aly A Farag.  Elastic shape registration using an incremental free form deformation approach with the icp algorithm.  In 2011 Canadian Conference on Computer and Robot Vision, pages 212–218. IEEE, 2011.
  * [2] Fernando J AGUILAR, Ismael FERNÁNDEZ, Juan A CASANOVA, Francisco J RAMOS, Manuel A AGUILAR, José L BLANCO, and José C MORENO.  3d coastal monitoring from very dense uav-based photogrammetric point clouds.  In Advances on Mechanics, Design Engineering and Manufacturing: Proceedings of the International Joint Conference on Mechanics, Design Engineering & Advanced Manufacturing (JCM 2016), 14-16 September, 2016, Catania, Italy, pages 879–887. Springer, 2017.
  * [3] Ibraheem Alhashim, Kai Xu, Yixin Zhuang, Junjie Cao, Patricio Simari, and Hao Zhang.  Deformation-driven topology-varying 3d shape correspondence.  ACM Transactions on Graphics (TOG), 34(6):1–13, 2015.
  * [4] Hatem Alismail, L Douglas Baker, and Brett Browning.  Continuous trajectory estimation for 3d slam from actuated lidar.  In 2014 IEEE International Conference on Robotics and Automation (ICRA), pages 6096–6101. IEEE, 2014.
  * [5] Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré.  Gradient flows: in metric spaces and in the space of probability measures.  Springer Science & Business Media, 2005.
  * [6] Yasuhiro Aoki, Hunter Goforth, Rangaprasad Arun Srivatsan, and Simon Lucey.  Pointnetlk: Robust & efficient point cloud registration using pointnet.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 7163–7172, 2019.
  * [7] Yikun Bai, Ivan Vladimir Medri, Rocio Diaz Martin, Rana Shahroz, and Soheil Kolouri.  Linear optimal partial transport embedding.  In International Conference on Machine Learning, pages 1492–1520. PMLR, 2023.
  * [8] Yikun Bai, Bernard Schmitzer, Mathew Thorpe, and Soheil Kolouri.  Sliced optimal partial transport.  arXiv preprint arXiv:2212.08049, 2022.
  * [9] Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré.  Iterative bregman projections for regularized transportation problems.  SIAM Journal on Scientific Computing, 37(2):A1111–A1138, 2015.
  * [10] Paul J Besl and Neil D McKay.  Method for registration of 3-d shapes.  In Sensor fusion IV: control paradigms and data structures, volume 1611, pages 586–606. Spie, 1992.
  * [11] Nicolas Bonneel and David Coeurjolly.  SPOT: sliced partial optimal transport.  ACM Transactions on Graphics, 38(4):1–13, 2019.
  * [12] Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister.  Sliced and Radon Wasserstein barycenters of measures.  Journal of Mathematical Imaging and Vision, 51(1):22–45, 2015.
  * [13] Nicolas Bonneel, Michiel Van De Panne, Sylvain Paris, and Wolfgang Heidrich.  Displacement interpolation using lagrangian mass transport.  In Proceedings of the 2011 SIGGRAPH Asia conference, pages 1–12, 2011.
  * [14] Fred L. Bookstein.  Principal warps: Thin-plate splines and the decomposition of deformations.  IEEE Transactions on pattern analysis and machine intelligence, 11(6):567–585, 1989.
  * [15] Benedict J Brown and Szymon Rusinkiewicz.  Global non-rigid alignment of 3-d scans.  In ACM SIGGRAPH 2007 papers, pages 21–es. Association for Computing Machinery, 2007.
  * [16] Luis A Caffarelli and Robert J McCann.  Free boundaries in optimal transport and monge-ampere obstacle problems.  Annals of mathematics, pages 673–730, 2010.
  * [17] Laetitia Chapel, Mokhtar Z Alaya, and Gilles Gasso.  Partial optimal tranport with applications on positive-unlabeled learning.  Advances in Neural Information Processing Systems, 33:2903–2913, 2020.
  * [18] Yongxin Chen, Tryphon T Georgiou, Lipeng Ning, and Allen Tannenbaum.  Matricial wasserstein-1 distance.  IEEE control systems letters, 1(1):14–19, 2017.
  * [19] Yang Chen and Gérard Medioni.  Object modelling by registration of multiple range images.  Image and vision computing, 10(3):145–155, 1992.
  * [20] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  Scaling algorithms for unbalanced optimal transport problems.  Mathematics of Computation, 87(314):2563–2609, 2018.
  * [21] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  Unbalanced optimal transport: Dynamic and Kantorovich formulations.  Journal of Functional Analysis, 274(11):3090–3123, 2018.
  * [22] Minsu Cho, Jian Sun, Olivier Duchenne, and Jean Ponce.  Finding matches in a haystack: A max-pooling strategy for graph matching in the presence of outliers.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2083–2090, 2014.
  * [23] Haili Chui and Anand Rangarajan.  A new point matching algorithm for non-rigid registration.  Computer Vision and Image Understanding, 89(2-3):114–141, 2003\.
  * [24] Ondrej Chum and Jiri Matas.  Matching with prosac-progressive sample consensus.  In 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05), volume 1, pages 220–226. IEEE, 2005.
  * [25] Andrew DJ Cross and Edwin R Hancock.  Graph matching with a dual-step em algorithm.  IEEE transactions on pattern analysis and machine intelligence, 20(11):1236–1253, 1998.
  * [26] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  Advances in neural information processing systems, 26, 2013.
  * [27] Steven B Damelin.  Near Extensions and Alignment of Data in Rn: Whitney extensions of near isometries, shortest paths, equidistribution, clustering and non-rigid alignment of data in Euclidean space.  John Wiley & Sons, 2024.
  * [28] Zheng Dang, Fei Wang, and Mathieu Salzmann.  Learning 3d-3d correspondences for one-shot partial-to-partial registration.  arXiv preprint arXiv:2006.04523, 2020.
  * [29] Bailin Deng, Yuxin Yao, Roberto M Dyke, and Juyong Zhang.  A survey of non-rigid 3d registration.  In Computer Graphics Forum, volume 41, pages 559–589. Wiley Online Library, 2022.
  * [30] Haowen Deng, Tolga Birdal, and Slobodan Ilic.  Ppfnet: Global context aware local features for robust 3d point matching.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 195–205, 2018.
  * [31] Haowen Deng, Tolga Birdal, and Slobodan Ilic.  3d local features for direct pairwise registration.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 3244–3253, 2019.
  * [32] Konstantinos G Derpanis.  Overview of the ransac algorithm.  Image Rochester NY, 4(1):2–3, 2010.
  * [33] Ishan Deshpande, Yuan-Ting Hu, Ruoyu Sun, Ayis Pyrros, Nasir Siddiqui, Sanmi Koyejo, Zhizhen Zhao, David Forsyth, and Alexander G Schwing.  Max-sliced wasserstein distance and its use for gans.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10648–10656, 2019.
  * [34] Yago Diez, Joan Martí, and Joaquim Salvi.  Hierarchical normal space sampling to speed up point cloud coarse matching.  Pattern Recognition Letters, 33(16):2127–2133, 2012.
  * [35] Csaba Domokos, Jozsef Nemeth, and Zoltan Kato.  Nonlinear shape registration without correspondences.  IEEE Transactions on pattern analysis and machine intelligence, 34(5):943–958, 2011.
  * [36] Shaoyi Du, Nanning Zheng, Shihui Ying, Qubo You, and Yang Wu.  An extension of the icp algorithm considering scale factor.  In 2007 IEEE International Conference on Image Processing, volume 5, pages V–193. IEEE, 2007.
  * [37] Jean Duchon.  Splines minimizing rotation-invariant semi-norms in sobolev spaces.  In Constructive theory of functions of several variables, pages 85–100. Springer, 1977.
  * [38] Marvin Eisenberger, Aysim Toker, Laura Leal-Taixé, and Daniel Cremers.  Deep shells: Unsupervised shape correspondence with optimal transport.  Advances in Neural information processing systems, 33:10491–10502, 2020.
  * [39] Kilian Fatras, Thibault Séjourné, Rémi Flamary, and Nicolas Courty.  Unbalanced minibatch optimal transport; applications to domain adaptation.  In International Conference on Machine Learning, pages 3186–3197. PMLR, 2021.
  * [40] Jean Feydy, Benjamin Charlier, François-Xavier Vialard, and Gabriel Peyré.  Optimal transport for diffeomorphic registration.  In Medical Image Computing and Computer Assisted Intervention- MICCAI 2017: 20th International Conference, Quebec City, QC, Canada, September 11-13, 2017, Proceedings, Part I 20, pages 291–299. Springer, 2017\.
  * [41] Alessio Figalli and Nicola Gigli.  A new transportation distance between non-negative measures, with applications to gradients flows with dirichlet boundary conditions.  Journal de mathématiques pures et appliquées, 94(2):107–130, 2010.
  * [42] Alessio Figalli and Federico Glaudo.  An invitation to optimal transport, Wasserstein distances, and gradient flows.  European Mathematical Society Publishing House, 2021.
  * [43] Martin A Fischler and Robert C Bolles.  Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography.  Communications of the ACM, 24(6):381–395, 1981.
  * [44] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer.  Pot: Python optimal transport.  Journal of Machine Learning Research, 22(78):1–8, 2021.
  * [45] Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis Bach.  Stochastic optimization for large-scale optimal transport.  Advances in neural information processing systems, 29, 2016.
  * [46] Steven Gold, Anand Rangarajan, Chien-Ping Lu, Suguna Pappu, and Eric Mjolsness.  New algorithms for 2d and 3d point matching: pose estimation and correspondence.  Pattern recognition, 31(8):1019–1031, 1998.
  * [47] Edouard Grave, Armand Joulin, and Quentin Berthet.  Unsupervised alignment of embeddings with wasserstein procrustes.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 1880–1890. PMLR, 2019.
  * [48] Kevin Guittet.  Extended Kantorovich norms: a tool for optimization.  PhD thesis, INRIA, 2002.
  * [49] Dirk Haehnel, Sebastian Thrun, and Wolfram Burgard.  An extension of the icp algorithm for modeling nonrigid objects with mobile robots.  In IJCAI, volume 3, pages 915–920, 2003.
  * [50] Osamu Hirose.  A bayesian formulation of coherent point drift.  IEEE transactions on pattern analysis and machine intelligence, 43(7):2269–2286, 2020.
  * [51] Hugues Hoppe, Tony DeRose, Tom Duchamp, John McDonald, and Werner Stuetzle.  Surface reconstruction from unorganized points.  In Proceedings of the 19th annual conference on computer graphics and interactive techniques, pages 71–78, 1992.
  * [52] Berthold KP Horn.  Closed-form solution of absolute orientation using unit quaternions.  Josa a, 4(4):629–642, 1987.
  * [53] Berthold KP Horn, Hugh M Hilden, and Shahriar Negahdaripour.  Closed-form solution of absolute orientation using orthonormal matrices.  JOSA A, 5(7):1127–1135, 1988.
  * [54] Qi-Xing Huang, Bart Adams, Martin Wicke, and Leonidas J Guibas.  Non-rigid registration under isometric deformations.  In Computer Graphics Forum, volume 27, pages 1449–1457. Wiley Online Library, 2008.
  * [55] Ruikun Huang, Junli Zhao, Fuqing Duan, Xin Li, Celong Liu, Xiaodan Deng, Zhenkuan Pan, Zhongke Wu, and Mingquan Zhou.  Automatic craniofacial registration based on radial curves.  Computers & Graphics, 82:264–274, 2019.
  * [56] Xiaoshui Huang, Guofeng Mei, Jian Zhang, and Rana Abbas.  A comprehensive survey on point cloud registration.  arXiv preprint arXiv:2103.02690, 2021.
  * [57] Veli Ilci and Charles Toth.  High definition 3d map creation using gnss/imu/lidar sensor integration to support autonomous vehicle navigation.  Sensors, 20(3):899, 2020.
  * [58] Bing Jian and Baba C Vemuri.  A robust algorithm for point set registration using mixture of gaussians.  In Tenth IEEE International Conference on Computer Vision (ICCV’05) Volume 1, volume 2, pages 1246–1251. IEEE, 2005.
  * [59] Kun Jin, Chaoyue Liu, and Cathy Xia.  Two-sided wasserstein procrustes analysis.  In IJCAI, pages 3515–3521, 2021.
  * [60] Andrew Edie Johnson and Sing Bing Kang.  Registration and integration of textured 3d data.  Image and vision computing, 17(2):135–147, 1999.
  * [61] Anupam Joshi and C-H Lee.  On the problem of correspondence in range data and some inelastic uses for elastic nets.  IEEE transactions on neural networks, 6(3):716–723, 1995.
  * [62] Leonid V Kantorovich.  On the translocation of masses.  In Dokl. Akad. Nauk. USSR (NS), volume 37, pages 199–201, 1942\.
  * [63] Leonid V Kantorovich.  On a problem of monge.  In CR (Doklady) Acad. Sci. URSS (NS), volume 3, pages 225–226, 1948\.
  * [64] Narendra Karmarkar.  A new polynomial-time algorithm for linear programming.  In Proceedings of the sixteenth annual ACM symposium on Theory of computing, pages 302–311, 1984.
  * [65] Soheil Kolouri, Kimia Nadjahi, Umut Simsekli, Roland Badeau, and Gustavo Rohde.  Generalized sliced Wasserstein distances.  Advances in Neural Information Processing Systems, 32, 2019.
  * [66] Soheil Kolouri, Se Rim Park, and Gustavo K Rohde.  The radon cumulative distribution transform and its application to image classification.  IEEE transactions on image processing, 25(2):920–934, 2015.
  * [67] Soheil Kolouri, Yang Zou, and Gustavo K Rohde.  Sliced wasserstein kernels for probability distributions.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5258–5267, 2016.
  * [68] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [69] Wonjun Lee, Rongjie Lai, Wuchen Li, and Stanley Osher.  Generalized unnormalized optimal transport and its fast algorithms.  Journal of Computational Physics, 436:110041, 2021.
  * [70] Christian Léonard.  A survey of the schr$\backslash$” odinger problem and some of its connections with optimal transport.  arXiv preprint arXiv:1308.0215, 2013.
  * [71] Zohar Levi and Craig Gotsman.  Smooth rotation enhanced as-rigid-as-possible mesh animation.  IEEE transactions on visualization and computer graphics, 21(2):264–277, 2014.
  * [72] Hao Li, Robert W Sumner, and Mark Pauly.  Global correspondence optimization for non-rigid registration of depth scans.  In Computer graphics forum, volume 27, pages 1421–1430. Wiley Online Library, 2008.
  * [73] Shiying Li and Caroline Moosmueller.  Measure transfer via stochastic slicing and matching.  arXiv preprint arXiv:2307.05705, 2023.
  * [74] Xingyu Liu, Charles R Qi, and Leonidas J Guibas.  Flownet3d: Learning scene flow in 3d point clouds.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 529–537, 2019.
  * [75] Antoine Liutkus, Umut Simsekli, Szymon Majewski, Alain Durmus, and Fabian-Robert Stöter.  Sliced-wasserstein flows: Nonparametric generative modeling via optimal transport and diffusions.  In International Conference on Machine Learning, pages 4104–4113. PMLR, 2019.
  * [76] Weixin Lu, Guowei Wan, Yao Zhou, Xiangyu Fu, Pengfei Yuan, and Shiyu Song.  Deepvcp: An end-to-end deep neural network for point cloud registration.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 12–21, 2019.
  * [77] Bin Luo and Edwin R Hancock.  A unified framework for alignment and correspondence.  Computer Vision and Image Understanding, 92(1):26–55, 2003.
  * [78] Jiayi Ma, Jia Wu, Ji Zhao, Junjun Jiang, Huabing Zhou, and Quan Z Sheng.  Nonrigid point set registration with robust transformation learning under manifold regularization.  IEEE transactions on neural networks and learning systems, 30(12):3584–3597, 2018.
  * [79] Jiayi Ma, Ji Zhao, and Alan L Yuille.  Non-rigid point set registration by preserving global and local structures.  IEEE Transactions on image Processing, 25(1):53–64, 2015.
  * [80] Graham McNeill and Sethu Vijayakumar.  A probabilistic approach to robust shape matching.  In 2006 International Conference on Image Processing, pages 937–940. IEEE, 2006.
  * [81] Guofeng Mei, Xiaoshui Huang, Litao Yu, Jian Zhang, and Mohammed Bennamoun.  Cotreg: Coupled optimal transport based point cloud registration.  arXiv preprint arXiv:2112.14381, 2021.
  * [82] Jean Meinguet.  Multivariate interpolation at arbitrary points made simple.  Zeitschrift für angewandte Mathematik und Physik ZAMP, 30(2):292–304, 1979.
  * [83] Ellon Mendes, Pierrick Koch, and Simon Lacroix.  Icp-based pose-graph slam.  In 2016 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR), pages 195–200. IEEE, 2016.
  * [84] Gaspard Monge.  Mémoire sur la théorie des déblais et des remblais.  Mem. Math. Phys. Acad. Royale Sci., pages 666–704, 1781.
  * [85] Andriy Myronenko and Xubo Song.  Point set registration: Coherent point drift.  IEEE transactions on pattern analysis and machine intelligence, 32(12):2262–2275, 2010.
  * [86] Andriy Myronenko, Xubo Song, and Miguel Carreira-Perpinan.  Non-rigid point set registration: Coherent point drift.  Advances in neural information processing systems, 19, 2006.
  * [87] Khai Nguyen and Nhat Ho.  Revisiting sliced wasserstein on images: From vectorization to convolution.  Advances in Neural Information Processing Systems, 35:17788–17801, 2022.
  * [88] James Orlin.  A faster strongly polynomial minimum cost flow algorithm.  In Proceedings of the Twentieth annual ACM symposium on Theory of Computing, pages 377–387, 1988.
  * [89] Zhichen Pan, Haoyao Chen, Silin Li, and Yunhui Liu.  Clustermap building and relocalization in urban environments for unmanned vehicles.  Sensors, 19(19):4252, 2019.
  * [90] Mark Pauly, Richard Keiser, and Markus Gross.  Multi-scale feature extraction on point-sampled surfaces.  In Computer graphics forum, volume 22, pages 281–289. Wiley Online Library, 2003.
  * [91] Gabriel Peyré, Marco Cuturi, et al.  Computational optimal transport: With applications to data science.  Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019.
  * [92] Benedetto Piccoli and Francesco Rossi.  Generalized wasserstein distance and its application to transport equations with source.  Archive for Rational Mechanics and Analysis, 211(1):335–358, 2014\.
  * [93] François Pomerleau, Francis Colas, Roland Siegwart, et al.  A review of point cloud registration algorithms for mobile robotics.  Foundations and Trends® in Robotics, 4(1):1–104, 2015.
  * [94] Vivien Potó, József Árpád Somogyi, Tamás Lovas, and Árpád Barsi.  Laser scanned point clouds to support autonomous vehicles.  Transportation Research Procedia, 27:531–537, 2017.
  * [95] Gilles Puy, Alexandre Boulch, and Renaud Marlet.  Flot: Scene flow on point clouds guided by optimal transport.  In Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part XXVIII, pages 527–544. Springer, 2020.
  * [96] Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas.  Pointnet: Deep learning on point sets for 3d classification and segmentation.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 652–660, 2017.
  * [97] Charles Ruizhongtai Qi, Li Yi, Hao Su, and Leonidas J Guibas.  Pointnet++: Deep hierarchical feature learning on point sets in a metric space.  Advances in neural information processing systems, 30, 2017.
  * [98] Hongxing Qin, Yucheng Zhang, Zhentao Liu, and Baoquan Chen.  Rigid registration of point clouds based on partial optimal transport.  In Computer Graphics Forum, volume 41, pages 365–378. Wiley Online Library, 2022.
  * [99] Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot.  Wasserstein barycenter and its application to texture mixing.  In International Conference on Scale Space and Variational Methods in Computer Vision, pages 435–446. Springer, 2011.
  * [100] Guillem Ramírez, Rumen Dangovski, Preslav Nakov, and Marin Soljačić.  On a novel application of wasserstein-procrustes for unsupervised cross-lingual learning.  arXiv preprint arXiv:2007.09456, 2020.
  * [101] Anand Rangarajan, Haili Chui, Eric Mjolsness, Suguna Pappu, Lila Davachi, Patricia Goldman-Rakic, and James Duncan.  A robust point-matching algorithm for autoradiograph alignment.  Medical image analysis, 1(4):379–398, 1997.
  * [102] Abtin Rasoulian, Robert Rohling, and Purang Abolmaesumi.  Group-wise registration of point sets for statistical shape models.  IEEE transactions on medical imaging, 31(11):2025–2034, 2012.
  * [103] Szymon Rusinkiewicz and Marc Levoy.  Efficient variants of the icp algorithm.  In Proceedings third international conference on 3-D digital imaging and modeling, pages 145–152. IEEE, 2001.
  * [104] Radu Bogdan Rusu and Steve Cousins.  3d is here: Point cloud library (pcl).  In 2011 IEEE international conference on robotics and automation, pages 1–4. IEEE, 2011.
  * [105] Vinit Sarode, Xueqian Li, Hunter Goforth, Yasuhiro Aoki, Rangaprasad Arun Srivatsan, Simon Lucey, and Howie Choset.  Pcrnet: Point cloud registration network using pointnet encoding.  arXiv preprint arXiv:1908.07906, 2019.
  * [106] Ryoma Sato, Marco Cuturi, Makoto Yamada, and Hisashi Kashima.  Fast and robust comparison of probability measures in heterogeneous spaces.  arXiv preprint arXiv:2002.01615, 2020.
  * [107] Meyer Scetbon and Marco Cuturi.  Low-rank optimal transport: Approximation, statistics and debiasing.  Advances in Neural Information Processing Systems, 35:6802–6814, 2022.
  * [108] Meyer Scetbon, Marco Cuturi, and Gabriel Peyré.  Low-rank sinkhorn factorization.  In International Conference on Machine Learning, pages 9344–9354. PMLR, 2021.
  * [109] Peter H Schönemann.  A generalized solution of the orthogonal procrustes problem.  Psychometrika, 31(1):1–10, 1966.
  * [110] Erwin Schrödinger.  Über die umkehrung der naturgesetze.  Verlag der Akademie der Wissenschaften in Kommission bei Walter De Gruyter u …, 1931.
  * [111] Thibault Séjourné, Clément Bonet, Kilian Fatras, Kimia Nadjahi, and Nicolas Courty.  Unbalanced optimal transport meets sliced-wasserstein.  arXiv preprint arXiv:2306.07176, 2023.
  * [112] Avinash Sharma, Radu Horaud, Jan Cech, and Edmond Boyer.  Topologically-robust 3d shape matching based on diffusion geometry and seed growing.  In CVPR 2011, pages 2481–2488. IEEE, 2011.
  * [113] Yiru Shen, Chen Feng, Yaoqing Yang, and Dong Tian.  Mining point cloud local structures by kernel correlation and graph pooling.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4548–4557, 2018.
  * [114] Zhengyang Shen, Jean Feydy, Peirong Liu, Ariel H Curiale, Ruben San Jose Estepar, Raul San Jose Estepar, and Marc Niethammer.  Accurate point cloud registration with robust optimal transport.  Advances in Neural Information Processing Systems, 34:5373–5389, 2021.
  * [115] Richard Sinkhorn.  A relationship between arbitrary positive matrices and doubly stochastic matrices.  The annals of mathematical statistics, 35(2):876–879, 1964.
  * [116] Martin Sinko, Patrik Kamencay, Robert Hudec, and Miroslav Benco.  3d registration of the point cloud data using icp algorithm in medical image analysis.  In 2018 ELEKTRO, pages 1–6. IEEE, 2018.
  * [117] Fabien Tâche, Wolfgang Fischer, Gilles Caprari, Roland Siegwart, Roland Moser, and Francesco Mondada.  Magnebike: A magnetic wheeled robot with high mobility for inspecting complex-shaped structures.  Journal of Field Robotics, 26(5):453–476, 2009.
  * [118] Art Tevs, Alexander Berner, Michael Wand, Ivo Ihrke, Martin Bokeloh, Jens Kerber, and Hans-Peter Seidel.  Animation cartography—intrinsic reconstruction of shape and motion.  ACM Transactions on Graphics (TOG), 31(2):1–15, 2012.
  * [119] Philip HS Torr and Andrew Zisserman.  Mlesac: A new robust estimator with application to estimating image geometry.  Computer vision and image understanding, 78(1):138–156, 2000.
  * [120] Yanghai Tsin and Takeo Kanade.  A correlation-based approach to robust point set registration.  In Computer Vision-ECCV 2004: 8th European Conference on Computer Vision, Prague, Czech Republic, May 11-14, 2004. Proceedings, Part III 8, pages 558–569. Springer Berlin Heidelberg, 2004.
  * [121] Cedric Villani.  Optimal transport: old and new.  Springer, 2009.
  * [122] Cédric Villani.  Topics in optimal transportation, volume 58.  American Mathematical Soc., 2021.
  * [123] Grace Wahba.  Spline models for observational data.  SIAM, 1990.
  * [124] Lingjing Wang, Jianchun Chen, Xiang Li, and Yi Fang.  Non-rigid point set registration networks.  arXiv preprint arXiv:1904.01428, 2019.
  * [125] Lingjing Wang, Xiang Li, Jianchun Chen, and Yi Fang.  Coherent point drift networks: Unsupervised learning of non-rigid point set registration.  arXiv preprint arXiv:1906.03039, 2019.
  * [126] Yue Wang and Justin M Solomon.  Deep closest point: Learning representations for point cloud registration.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 3523–3532, 2019.
  * [127] Zi-Ming Wang, Nan Xue, Ling Lei, and Gui-Song Xia.  Partial wasserstein adversarial network for non-rigid point set registration.  arXiv preprint arXiv:2203.02227, 2022.
  * [128] Alan Geoffrey Wilson.  The use of entropy maximising models, in the theory of trip distribution, mode split and route split.  Journal of transport economics and policy, pages 108–126, 1969\.
  * [129] Jinzhong Yang.  The thin plate spline robust point matching (tps-rpm) algorithm: A revisit.  Pattern Recognition Letters, 32(7):910–918, 2011.
  * [130] Zixin Yang, Richard Simon, and Cristian A Linte.  Learning feature descriptors for pre-and intra-operative point cloud matching for laparoscopic liver registration.  arXiv preprint arXiv:2211.03688, 2022.
  * [131] Zi Jian Yew and Gim Hee Lee.  Rpm-net: Robust point matching using learned features.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 11824–11833, 2020.
  * [132] Chi Yuan, Xiaoqing Yu, and Ziyue Luo.  3d point cloud matching based on principal component analysis and iterative closest point algorithm.  In 2016 International Conference on Audio, Language and Image Processing (ICALIP), pages 404–408. IEEE, 2016.
  * [133] Wentao Yuan, Benjamin Eckart, Kihwan Kim, Varun Jampani, Dieter Fox, and Jan Kautz.  Deepgmr: Learning latent gaussian mixture models for registration.  In Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part V 16, pages 733–750. Springer, 2020.
  * [134] Alan L Yuille and Norberto M Grzywacz.  A computational theory for the perception of coherent visual motion.  Nature, 333(6168):71–74, 1988.
  * [135] Andy Zeng, Shuran Song, Matthias Nießner, Matthew Fisher, Jianxiong Xiao, and Thomas Funkhouser.  3dmatch: Learning local geometric descriptors from rgb-d reconstructions.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 1802–1811, 2017.
  * [136] Juyong Zhang, Yuxin Yao, and Bailin Deng.  Fast and robust iterative closest point.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021\.
  * [137] Zhengyou Zhang.  Iterative point matching for registration of free-form curves and surfaces.  International journal of computer vision, 13(2):119–152, 1994.
  * [138] Zhiyuan Zhang, Yuchao Dai, and Jiadai Sun.  Deep learning based point cloud registration: an overview.  Virtual Reality & Intelligent Hardware, 2(3):222–246, 2020.
  * [139] Timo Zinßer, Jochen Schmidt, and Heinrich Niemann.  Point set registration with integrated scale estimation.  In International conference on pattern recognition and image processing, pages 116–119, 2005.

##  Appendix A Notations

Below, we set the notations consistently used in the main text and appendix in
this paper.

  * •

$\mathbb{R}_{+}=[0,\infty)$, set of all non-negative real numbers.

  * •

$\mathbb{R}_{++}=(0,\infty)$, set of all positive real numbers.

  * •

$\mathbb{R}^{D}$, the space of the point clouds. In this paper, without
specification, we set $D=2$ or $D=3$.

  * •

Source point set:

$X=\\{x^{n}\in\mathbb{R}^{D}\\}_{n=1}^{N}=[x^{1},\ldots,x^{N}]^{T}\in\mathbb{R}^{N\times
D}$

  * •

Target point set:

$Y=\\{y^{m}\in\mathbb{R}^{D}\\}_{m=1}^{M}=[y^{1},\ldots,y^{M}]^{T}\in\mathbb{R}^{M\times
D}$

  * •

Control points:

$C=\\{c^{k}\in\mathbb{R}^{D}\\}_{k=1}^{K}=[c^{1},\ldots,c^{K}]^{T}\in\mathbb{R}^{K\times
D}$

  * •

$A[n,:]$, the $n^{th}$ row of matrix $A$.

  * •

$\mathcal{D}\subset[1:N]$, a set of indices.

  * •

$A[\mathcal{D},:]$, the sub-matrix of $A$, formed by taking the rows
corresponding to the indices in $\mathcal{D}$, i.e. if $n\in\mathcal{D}$, then
$A[n,:]$ is in $A[\mathcal{D}]$.

  * •

$I_{D}$, $D\times D$ identity matrix.

  * •

$1_{N}$, $N\times 1$ vector with all 1 entries.

  * •

$0_{N}$, $N\times 1$ vector with all 0 entries.

  * •

$x_{d}$, the $d^{th}$ component of $x\in\mathbb{R}^{D}$.

  * •

$f:\mathbb{R}^{D}\mapsto\mathbb{R}^{D}$, is a deformation model that maps
source point cloud to target point cloud. In general

$f(x)=h(x)+B^{T}x+\beta$

where $h(x)$ is the nonrigid part of the model, $B\in\mathbb{R}^{D\times D}$
form the affine part and $\beta\in\mathbb{R}^{D}$ is the translation
parameter.

On particular case of the above model is

$f(x)=h(x)+R^{T}Sx+\beta.$

In this formulation, $S=\text{diag}(s_{1},\ldots s_{D})\in\mathbb{R}^{D\times
D}$ is a diagonal scaling matrix where each $s_{i}>0$. And
$R\in\mathbb{R}^{D\times D}$ is a $D$-dimensional rotation matrix.
It is worth noting that the model presented in [36] uses the formulation
$RSx+\beta$. To maintain consistency between the RBF and TPS models, we opt
for the $R^{T}SX+\beta$ formulation. Thus, the ICP step (35) is slightly
different from the equation (17) in this reference.

  * •

$f[d]$, $d^{th}$ component of $f$, where $d\in[1:D]$.

  * •

$O_{D}$, set of all $D\times D$ orthornormal matrices.

  * •

$\hat{Y}:=[\hat{y}^{1},\ldots,\hat{y}^{N}]^{T}\in\mathbb{R}^{N\times D}$, the
estimated target points, i.e. for each $n$, $\hat{y}^{n}=f(x^{n})$ where $f$
is a selected model.

  * •

correspondence between source and target:

In this paper, the phrase “correspondence between source and target”
encompasses two distinct but interconnected concepts.

First, during the step of solving the OPT problem (13) or (57), we have the
deformed source point set denoted as $\\{\hat{y}^{n}=f(x^{n})\\}_{n=1}^{N}$.
In this context, we say the optimal transportation plan $\gamma$ is a (soft)
correspondence between the (deformed) source points and the target points.

Second, after updating a subset of $\hat{y}$ points using either (15) or (58),
the set $\\{\hat{y}^{n}\\}_{n\in\mathcal{D}}$ is transformed to align with (a
subset of) the target set $\\{y_{m}\\}_{m=1}^{M}$. In this transformed state,
the point set $\\{\hat{y}\\}_{n\in\mathcal{D}}$ becomes a “representation” of
the target points. As a result, we refer to the set of pairs
$\\{(x^{n},\hat{y}^{n})\\}_{n\in\mathcal{D}}$ as a (partial) correspondence
between the source and the (representation of) target.

  * •

$\tilde{f}:\mathbb{R}^{D}\to\mathbb{R}$, the one dimensional model defined as

$\tilde{f}(x)=\tilde{h}(x)+b^{T}x+\tilde{\beta}$

where $\tilde{h}$ is the nonlinear deformation, $b\in\mathbb{R}^{D}$ and
$\tilde{\beta}\in\mathbb{R}$. Thus point deformation model $f$ can be regarded
as $(\tilde{f}_{1},\ldots\tilde{f}_{d})$, i.e.,
$f(x)=[\tilde{f}_{1}(x),\ldots,\tilde{f}_{d}(x)]^{T}$.

  * •

$\|\cdot\|$, $L_{2}$ norm in $\mathbb{R}^{D}$.

  * •

$\tilde{Y}=\\{\tilde{y}^{n}\in\mathbb{R}\\}_{n=1}^{N}=[\tilde{y}^{1},\ldots,\tilde{y}_{n}]^{T}$,
target points for $1D$ interpolation problem

$\displaystyle\min_{\tilde{f}}\sum_{n=1}^{n}\|\tilde{f}(x^{n})-\tilde{y}^{n}\|^{2}+\text{reg}(\tilde{f})$
(60)

where $\text{reg}(\tilde{f})$ is a regularization term. In this paper, the
above problem particularly becomes (38) and (41) when we discuss RBF model and
TPS model respectively.

  * •

$\phi(\cdot,\cdot),\phi_{G}(\cdot,\cdot),\phi_{T}(\cdot,\cdot)$: $\phi$ is a
kernel function, in particular:

$\phi_{G}(x,c;\sigma^{2})=e^{-\frac{\|c-x\|^{2}}{\sigma^{2}}}$

is the Gaussian kernel function with $\sigma$ denoting the kernel width;

$\phi_{T}(x,c;D)=\begin{cases}(\|x-c\|^{2})\ln(\|x-c\|)&\text{if }D=2\\\
\|x-c\|&\text{if }D=3\\\ \|x-c\|^{4-D}&\text{otherwise}\end{cases}$

is the Thin-Plate-Spline (TPS) kernel function.

  * •

$\Phi\in\mathbb{R}^{N\times K}$, $\Phi_{sub}=\Phi[\mathcal{D},:]$, where
$\mathcal{D}\subset[1:N]$: Kernel matrices where
$\Phi_{k,n}=\phi(x_{n},c_{k})$ and sub-matrix whose row is selected from
$\Phi$ via indices in $\mathcal{D}$.

  * •

$\alpha^{1},\ldots,\alpha^{K}\in\mathbb{R}^{D}$: the parameters for the kernel
mapping

$x\mapsto h(x):=\sum_{k=1}^{K}\alpha^{k}\phi(x,c^{k})$

in the point registration problem.

  * •

$\tilde{\alpha}_{1},\ldots,\tilde{\alpha}_{K}\in\mathbb{R}$: parameters of the
kernel mapping

$\tilde{x}\mapsto\tilde{h}(x):=\sum_{n=1}^{n}\alpha^{k}\phi(x,\tilde{c}^{k})$

in the 1D corresponding kernel interpolation problem (60).

  * •

$\mathcal{M}_{+}(\mathbb{R}^{D})$: the set of positive finite Radon measures
defined on $\mathbb{R}^{D}$.

  * •

$\mathcal{P}(\mathbb{R}^{D})$: the set of all probability measures defined on
$\mathbb{R}^{D}$.

  * •

$\mu,\nu\in\mathcal{M}_{+}(\mathbb{R}^{D})$: two positive finite Radon
measures.

  * •

$\pi_{1},\pi_{2}:(\mathbb{R}^{D})^{2}\to\mathbb{R}^{D}$: the canonical
projection mapping,

$\pi_{1}(\hat{y},y)=\hat{y};\pi_{2}(\hat{y},y)=y.$

  * •

$A$: A Borel set in space $\mathbb{R}^{D}$ or $(\mathbb{R}^{D})^{2}$.

  * •

$f_{\\#}\gamma$ where $\gamma\in\mathcal{M}_{+}((\mathbb{R}^{D})^{2})$,
$f:(\mathbb{R}^{D})^{2}\to\mathbb{R}^{D}$. The pushforward measure of $\gamma$
under $f$. In particular, for each Borel set $A\subseteq\mathbb{R}^{D}$,

$f_{\\#}\gamma(A)=\gamma(f^{-1}(A)).$

  * •

$|\mu|$: total mass of measure $\mu$.

  * •

$\Pi(\mu,\nu):=\\{\gamma\in\mathcal{P}((\mathbb{R}^{D})^{2}):\pi_{1\\#}\gamma=\mu,\pi_{2\\#}\gamma=\nu\\}$.

  * •

$\Pi(\frac{1}{N}1_{N},\frac{1}{M}1_{M}):=\\{\gamma\in\mathbb{R}_{+}^{N\times
M}:\gamma\frac{1}{M}1_{M}=\frac{1}{N}1_{N},\gamma^{T}\frac{1}{N}1_{N}=\frac{1}{M}1_{M}\\}$,
where $1_{N}$ ($1_{M}$) is $N$-size ($M$-size) vector whose entries are $1$.

  * •

$\mu\leq\nu$: denotes that $\mu$ is dominated by $\nu$, i.e., for each Borel
set $A\subseteq\mathbb{R}^{D}$, $\mu(A)\leq\nu(A)$.

  * •

$\Pi_{\leq}(\mu,\nu):=\\{\gamma\in\mathcal{M}_{+}((\mathbb{R}^{D})^{2}):\pi_{1\\#}\gamma\leq\mu,\pi_{2\\#}\gamma\leq\nu\\}.$

  * •

Optimal transport problem between $\mu,\nu\in\mathcal{P}(\mathbb{R}^{D})$:

$OT(\mu,\nu):=\inf_{\gamma\in\Pi(\mu,\nu)}\int
c(\hat{y},y)d\gamma(\hat{y},y).$

  * •

Empirical Optimal transport problem:

$\mu=\frac{1}{N}\sum_{n=1}^{N}\delta_{\hat{y}^{n}},\nu=\frac{1}{M}\sum_{m=1}^{M}\delta_{y^{m}}$,
where $\delta_{\hat{y}^{n}}$ denotes the Dirac measure centered at $x^{n}$
(and similarly for $\delta_{y^{m}}$).

$\text{OT}(\mu,\nu):=\min_{\gamma\in\Pi(\frac{1}{N}1_{N},\frac{1}{M}1_{M})}\sum_{n=1}^{N}\sum_{m=1}^{M}c(\hat{y}^{n},y^{m})\gamma_{n,m}$

  * •

OT Barycentric projection of $\nu$ with respect to $\mu$: Given $\gamma$, an
optimal solution for above empirical OT problem, the following empirical
measure is called baricentric projection of $\nu$ (with respect to $\mu$):

$\hat{\nu}:=\sum_{n=1}^{N}\frac{1}{N}\delta_{(M\gamma Y)[n:]^{T}}.$

In this formula, $M\gamma Y\in\mathbb{R}^{N\times D}$ and the $N$-points set
$\\{M\gamma Y[1]^{T},\ldots M\gamma Y[N]^{T}\\}$ can be regarded as a
“representation” of target point set $Y=\\{y^{1},\ldots y^{m}\\}$.

  * •

Optimal partial transport problem between
$\mu,\nu\in\mathcal{M}_{+}(\mathbb{R}^{D})$:

$OPT_{\lambda}(\mu,\nu):=\inf_{\gamma\in\Pi_{\leq}(\mu,\nu)}\int
c(\hat{y},y)d\gamma(\hat{y},y)+\lambda(|\mu|+|\nu|-2|\gamma|).$

where $\lambda\in\mathbb{R}_{+}$ is the penalty term for the creation and the
destruction of mass.

  * •

Empirical optimal partial transport:
$\mu=\sum_{n=1}^{N}\delta_{\hat{y}^{n}},\nu=\sum_{m=1}^{M}\delta_{y^{m}}$.

$\displaystyle\text{OPT}_{\lambda}(\mu,\nu):$
$\displaystyle=\min_{\gamma\in\Pi_{\leq}(1_{N},1_{M})}\sum_{n=1}^{N}\sum_{m=1}^{M}c(\hat{y}^{n},y^{m})\gamma_{n,m}+\lambda(N+M-2|\gamma|).$

  * •

Primal optimal partial transport:

$\displaystyle\text{primal-
OPT}(\mu,\nu;\zeta):=\inf_{\gamma\in\Pi_{\leq}(\mu,\nu)}$ $\displaystyle\int
c(\hat{y},y)d\gamma(\hat{y},y)$ $\displaystyle\text{s.t. }|\gamma|=\zeta.$

where $\zeta\in[0,\min\\{|\mu|,|\nu|\\}]$ is the mass constraint constant.

  * •

Empirical primal optimal partial transport:

$\displaystyle\text{primal-
OPT}(\mu,\nu;\zeta):=\inf_{\gamma\in\Pi_{\leq}(1_{N},1_{M})}$
$\displaystyle\sum c(\hat{y}^{n},y^{m})\gamma_{n,m}(\hat{y},y)$
$\displaystyle\text{s.t. }|\gamma|=\zeta.$

  * •

$L$: The Monge mapping in the OT problem or the optimal partial bijection in
empirical OPT problem. In particular, it means

$L:[1:N]\supset\text{Dom}(L)\to\text{Range}(L)\subset[1:M],$

or denoted as

$L:[1:N]\hookrightarrow[1:M],$

where $\text{Dom}(L),\text{Range}(L)$ is the domain and range of mapping $L$
respectively. In addition, $L$ is a bi-jection. We refer section III for the
detailed introduction and explanation of Monge mapping and corresponding
concepts.

  * •

$\mathbb{S}^{D-1}:=\\{\theta\in\mathbb{R}^{D}:\|\theta\|^{2}=1\\}$: the
uniform sphere in $\mathbb{R}^{D}$.

  * •

$\theta_{1},\ldots,\theta_{t^{\prime}}\in\mathbb{S}^{D-1}$: the projection
parameters we use for sliced-OT or sliced-OPT.

  * •

$T\in\mathbb{N}$: maximum number of iterations for a point cloud registration
method.

  * •

$\bar{X}=[Q_{1},Q_{2}]\begin{bmatrix}\mathcal{R}\\\ 0_{(N-(D+1))\times
N}\end{bmatrix}$, QR decomposition of $\bar{X}$.

On the left hand side, $\bar{X}:=[1_{N},X]$.

On the right hand side,
$Q_{1}\in\mathbb{R}^{N\times(D+1)},Q_{2}\in\mathbb{R}^{N\times(N-(D+1))}$,
$\mathcal{R}\in\mathbb{R}^{(D+1)\times(D+1)}$. $[Q_{1},Q_{2}]$ is an $N\times
N$ orthonormal matrix, and $\mathcal{R}$ is upper triangular matrix.

  * •

$\delta>0$: a constant used in the computation of primal-OPT problem. We refer
Appendix B for more details.

  * •

$\xi>0$: the weight parameter of entropic regularization in the Sinkhorn OT
(entropic OT) problem.

##  Appendix B Primal form of optimal partial transport problem and its
computation

Similar to the OPT problem (4) as we introduced in section III, there is a
one-by-one relation between the Primal form OPT (8) and the classical OT
problem (see Proposition 1 [17]). In particular, let
$\tilde{\mu}:=\mu+(\|\nu\|-\zeta)\delta_{\tilde{\infty}},\tilde{\nu}:=\nu+(\|\mu\|-\zeta)\delta_{\tilde{\infty}}$,
where $\tilde{\infty}$ is an auxiliary point. Define cost function $\tilde{c}$
with

$\displaystyle\tilde{c}(\hat{y},y)=\begin{cases}c(\hat{y}^{n},y^{m})&\text{if
}\hat{y}=\hat{y}^{n},y=y^{m}\\\ \lambda&\text{if
}\hat{y}=\hat{y}^{n},y=\tilde{\infty}\text{ or
}\hat{y}=\tilde{\infty},y=y^{m}\\\ 2\lambda+\delta&\text{if
}\hat{y}=y=\tilde{\infty}\end{cases}$ (61)

where $\lambda,\delta>0$. We define OT problem

$\text{OT}(\tilde{\mu},\tilde{\nu})=\min_{\tilde{\gamma}\in\Pi(\tilde{\mu},\tilde{\nu})}\sum_{n,m}\tilde{c}(\hat{y}^{n},y^{m})\tilde{\gamma}_{ij}$

where $x_{N+1},y_{M+1}=\tilde{\infty}$. We have

###### Proposition 1

The following function is a bijection between $\Pi(\tilde{\mu},\tilde{\nu})$
and $\Pi_{\leq}(\mu,\nu)$

$\Pi(\tilde{\mu},\tilde{\nu})\ni\tilde{\gamma}\mapsto\gamma[1:N,1:M]\in\Pi_{\leq}(\mu,\nu).$

One transportation plan $\tilde{\gamma}$ is optimal for above OT problem iff
$\tilde{\gamma}[1:N,1:M]$ is optimal for (8)

By the above relation between primal-OPT and OT, one can use the linear
programming solver [13] to solve the Primal form OPT problem, whose complexity
is $\mathcal{O}((N+M)NM)$.

##  Appendix C Relation between TPS-RPM method and Entropic OT problem

In this section, we discuss the relations between Entropic Optimal Transport
(EOT) and TPS-RPM. Based on this relation, we propose an improved version of
the TPS-RPM method.

Entropic optimal transport. The EOT problem, also known as the Sinkhorn
distance, was initially studied in the Schrodinger problem [110, 128, 70]. It
gained popularity in the machine learning community through the influential
paper[26], which demonstrates the significant computational speed-up
achievable using the Sinkhorn-Knopp algorithm. In classical EOT, the problem
is formulated between two probability measures, $\mu$ and $\nu$, both of which
belong to the space $\mathcal{P}(\mathbb{R}^{D})$, as follows:

$\displaystyle\text{EOT}_{\xi}(\mu,\nu):=\inf_{\gamma\in\Pi(\mu,\nu)}\int_{(\mathbb{R}^{D})^{2}}c(\hat{y},y)d\gamma+\xi
I(\gamma),$ (62)

where $\xi\in(0,\infty)$ is a constant and

$\displaystyle
I(\gamma):=\text{KL}(\gamma\parallel\mu\otimes\nu):=\int_{(\mathbb{R}^{D})^{2}}\ln\left(\frac{d\gamma}{d\mu\otimes
d\nu}\right)d\gamma.$ (63)

The term $I(\gamma)$ represents the “mutual information”, which is strictly
convex with respect to $\gamma$. When we consider the (normalized) $\gamma$ as
a joint distribution over the pair of random variables
$(\hat{Y}_{\mu},Y_{\nu})$, $I(\gamma)$ quantifies the amount of information
about $Y_{\nu}$ that is revealed when $\hat{Y}_{\mu}$ is observed, or vice
versa.

In the context of empirical distributions, where
$\mu=\sum_{n=1}^{N}\frac{1}{N}\delta_{\hat{y}^{n}}$ and
$\nu=\sum_{m=1}^{M}\frac{1}{M}\delta_{y^{m}}$, the EOT problem under
consideration becomes:

$\displaystyle EOT_{\xi}(\mu,\nu)$
$\displaystyle:=\inf_{\gamma\in\Pi(1_{N},1_{M})}\sum_{n,m}c(\hat{y}^{n},y^{m})\gamma_{n,m}+\xi\sum_{n,m}\ln(\gamma_{n,m})\gamma_{n,m}+\mathcal{C}_{1}$
$\displaystyle\simeq\inf_{\gamma\in\Pi(1_{N},1_{M})}\sum_{n,m}c(\hat{y}^{n},y^{m})\gamma_{n,m}+\xi\sum_{n,m}\ln(\gamma_{n,m})\gamma_{n,m}$
(64)

where $\mathcal{C}_{1}^{\prime}=\xi\sum_{n,m}\ln(NM)$ is a constant term. The
EOT problem (62) can be solved via the Sinkhorn-Knopp algorithm [91].

Entropic optimal partial transport. The EOT problem can be naturally extended
into the optimal partial transport setting [20]. In particular, given
$\mu,\nu\in\mathcal{M}_{+}(\Omega)$, the entropic OPT problem (EOPT) is
defined as

$\displaystyle\text{EOPT}_{\lambda,\xi}(\mu,\nu):=\inf_{\gamma\in\Pi_{\leq}(\mu,\nu)}\int_{(\mathbb{R}^{D})^{2}}c(\hat{y},y)d\gamma(\hat{y},y)+\lambda(|\mu-\pi_{1}\gamma|+|\nu-\pi_{2}\gamma|)+\xi
I(\gamma).$ (65)

and in the empirical setting, it becomes

$\displaystyle\text{EOPT}_{\lambda,\xi}(\mu,\nu)$
$\displaystyle:=\min_{\gamma\in\Pi_{\leq}(1_{N},1_{M})}\sum_{n,m}c(\hat{y}^{n},y^{m})\gamma_{n,m}+\lambda(N+M-2|\gamma|)+\xi\sum_{n,m}\gamma_{n,m}\ln(\gamma_{n,m})+\mathcal{C}_{1}$
$\displaystyle\simeq\min_{\gamma\in\Pi_{\leq}(1_{N},1_{M})}\sum_{n,m}c(\hat{y}^{n},y^{m})\gamma_{n,m}+\lambda(N+M-2|\gamma|)+\xi\sum_{n,m}\gamma_{n,m}\ln(\gamma_{n,m})$
$\displaystyle=\min_{\tilde{\gamma}\in\tilde{\Pi}(1_{N},1_{M})}\sum_{n,m}c(\hat{y}^{n},y^{m})\tilde{\gamma}_{n,m}-2\lambda\sum_{n=1}^{N}\sum_{m=1}^{M}\tilde{\gamma}_{n,m}+\xi\sum_{n=1}^{N}\sum_{m=1}^{M}\tilde{\gamma}_{n,m}\ln(\tilde{\gamma}_{n,m})+\mathcal{C}_{2}$
(66)

where in the last line, $\tilde{\Pi}\subset\mathbb{R}_{+}^{(N+1)\times(M+1)}$
is defined as set of $\tilde{\gamma}$ such that:

$\displaystyle\gamma 1_{M+1}[1:N]=1_{N},\gamma^{T}1_{N+1}[1:M]=1_{M}$
$\displaystyle\gamma 1_{M+1}[N+1]=M,\gamma^{T}1_{N+1}[M+1]=N.$ (67)

$\mathcal{C}_{2}=\lambda(N+M)$ is a constant. This equation holds from the
relation between the proposition 1, we also refer [16, 8] for more details.
Note, as an optimization problem, the constraint (67) can be released since
this constraint only involves the last row and column of $\tilde{\gamma}$,
i.e., $\tilde{\gamma}[N+1,:]$ and $\tilde{\gamma}[:,M+1]$, which is not
contained in the objective function (66). Remove the constant $\mathcal{C}$,
the problem (66) becomes

$\displaystyle\min_{\tilde{\gamma}}$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}c(\hat{y}^{n},y^{m})\tilde{\gamma}_{ij}-2\lambda\sum_{n=1}^{N}\sum_{m=1}^{M}\tilde{\gamma}_{ij}\ln(\tilde{\gamma}_{ij})+\xi\sum_{n=1}^{N}\sum_{m=1}^{M}\tilde{\gamma}_{ij}\ln(\tilde{\gamma}_{ij})$
(68) $\displaystyle\text{s.t.
}\tilde{\gamma}1_{M}[1:N]=1_{N},\tilde{\gamma}^{T}1_{N}[1:M]=1_{M}.$

Understanding TPS-RPM from the point of view of EOPT. In TPS-RPM [23, 129]
method, as we discussed before, the TPS model is defined as (40):

$\displaystyle f(x)=\sum_{n=1}^{N}\alpha^{n}\phi(x,x^{n})+B^{T}x+\beta$ (69)

and we aim to solve the following problem:

$\displaystyle\min_{f,\hat{\gamma}}$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}\|y^{m}-f(x^{n})\|^{2}\hat{\gamma}_{ij}+\lambda\|\nabla^{2}f\|^{2}-2\lambda\sum_{n=1}^{N}\sum_{m=1}^{M}\hat{\gamma}_{ij}+\xi\sum_{n=1}^{N}\sum_{m=1}^{M}\hat{\gamma}_{ij}\ln(\hat{\gamma}_{ij})$
(70)
$\displaystyle\text{s.t.}\hat{\gamma}\in\mathbb{R}_{+}^{(N+1)\times(M+1)},\hat{\gamma}1_{M}[1:N]=1_{N},\gamma^{T}1_{N}[1:M]=1_{M}.$

where $x^{N+1},y^{M+1}$ are auxiliary variables, called “outlier center”.

The method solves the problem described above in an iterative manner. Each
iteration consists of two steps: Step 1 involves updating the correspondence
$\hat{\gamma}$, and Step 2 involves updating $f$. Step 2 has been discussed in
section IV B and IV D. Step 1 is essentially the entropic OPT problem (68).
Indeed, when model $f$ is fixed,

$\displaystyle\min_{\hat{\gamma}}\sum_{n=1}^{N}\sum_{m=1}^{M}\|y^{m}-f(x^{n})\|^{2}\hat{\gamma}_{ij}+\lambda\|\nabla^{2}f\|^{2}-2\lambda\sum_{n=1}^{N}\sum_{m=1}^{M}\hat{\gamma}_{ij}+\xi\sum_{n=1}^{N}\sum_{m=1}^{M}\hat{\gamma}_{ij}\ln(\hat{\gamma}_{ij})$
$\displaystyle\simeq$
$\displaystyle\min_{\hat{\gamma}}\sum_{n=1}^{N}\sum_{m=1}^{M}\|y^{m}-f(x^{n})\|^{2}\hat{\gamma}_{ij}-2\lambda\sum_{n=1}^{N}\sum_{m=1}^{M}\hat{\gamma}_{ij}+\xi\sum_{n=1}^{N}\sum_{m=1}^{M}\hat{\gamma}_{ij}\ln(\hat{\gamma}_{ij})$

One can either solve it via equations (3)–(7) in [23], or the Sinkhorn
algorithm for entropic OPT problem [20, Algorithm 1].

An improved version of TPS-RPM. By the above understanding, we consider the
following entropic OPT problem:

$\displaystyle\min_{\gamma\in\Pi_{\leq}(\mu,\nu)}\int_{(\mathbb{R}^{D})^{2}}c(\hat{y},y)d\gamma+\xi
I(\gamma)$ $\displaystyle\quad\quad\quad\quad\text{s.t. }|\gamma|\geq\zeta.$
(71)

In the empirical setting, the above problem becomes (up to a constant)

$\displaystyle\min_{\gamma\in\Pi_{\leq}(1_{N},1_{M})}\sum_{n,m}c(\hat{y}^{n},y^{m})\gamma_{n,m}+\xi\sum_{ij}\gamma_{n,m}\ln(\gamma_{n,m})$
(72) $\displaystyle\quad\quad\quad\quad\quad\quad\text{s.t.
}\sum_{n=1}\sum_{m=1}^{M}\gamma_{n,m}=\zeta.$

Recall the relation between the primal-OPT and OT, we further simplify (72) as

$\displaystyle\min_{\hat{\gamma}}$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}c(\hat{y}^{n},y^{m})\hat{\gamma}_{ij}+\lambda^{\prime}(\sum_{n=1}^{N}\hat{\gamma}_{i,M+1}+\sum_{n=1}^{M}\hat{\gamma}_{N+1,j})+\delta\hat{\gamma}_{N+1,M+1}+\xi\sum_{n=1}^{N}\sum_{m=1}^{M}\hat{\gamma}_{ij}\ln(\hat{\gamma}_{ij})$
(73) $\displaystyle\text{s.t.
}\hat{\gamma}\in\mathbb{R}_{+}^{(N+1)\times(M+1)},\hat{\gamma}1_{M}[1:N]=1_{N},\hat{\gamma}^{T}1_{N}[1:M]=1_{M}$
$\displaystyle\hskip
15.00002pt\hat{\gamma}1_{M}[N+1]=M-\zeta,\hat{\gamma}^{T}1_{N}[M+1]=N-\zeta.$

where $\lambda^{\prime},\delta>0$ are constants. Note, un-rigorously speaking,
the above problem is a balanced entropic OT problem and it can be solved by
the Sinkhorn algorithm.

Thus, we propose a new version of TPS-RPM method. Consider the following
problem:

$\displaystyle\min_{f,\gamma}$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}\|y^{m}-f(x^{n})\|^{2}\gamma_{n,m}+\epsilon\|\nabla^{2}f\|^{2}+\xi\sum_{n=1}^{N}\sum_{m=1}^{M}\gamma_{n,m}\ln(\gamma_{n,m})$
(74) $\displaystyle\text{s.t. }\gamma\in\mathbb{R}_{+}^{N\times
M},\hat{\gamma}1_{M}\leq
1_{N},\gamma^{T}1_{N}=1_{M},\sum_{n=1}^{N}\sum_{m=1}^{M}\gamma_{n,m}=\zeta.$

In step of updating correspondence $\gamma$, it is equivalent to the entropic
primal-OPT problem:

$\displaystyle\min_{\gamma}$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}\|y^{m}-f(x^{n})\|^{2}\gamma_{n,m}+\xi\sum_{n=1}^{N}\sum_{m=1}^{M}\ln(\gamma_{n,m})\gamma_{n,m}$
$\displaystyle\text{s.t. }\gamma\in\mathbb{R}^{N\times M},\gamma 1_{M}\leq
1_{N},\gamma^{T}1_{N}\leq 1_{M},|\gamma|=\zeta,$

One can either solve it via Sinkhorn algorithm for primal OPT problem [9] or
convert it into the classical entropic OT problem (73)

##  Appendix D 2D Experiment

In this section, we present an evaluation of various baseline methods,
including CPD[85], the four non-rigid adaptations of Wasserstein Procrustes
methods [47]: OT-BBF, OT-TPS, SOT-RBF, and SOT-TPS, and the classical TPS-RPM
[23], along with its modified variant TPS-RPM(new) discussed in the previous
section, and the four novel methods proposed in Section IV in 2D setting. The
dataset for this experiment is the synthetic tropical fish, available at
https://github.com/siavashk/pycpd. In this section, we present an evaluation
of various baseline methods, including CPD[85], the four non-rigid adaptations
of Wasserstein Procrustes methods [47]: OT-BBF, OT-TPS, SOT-RBF, and SOT-TPS.
Additionally, we include TPS-RPM [23], its modified variant TPS-RPM(new)
discussed in the previous section, and the four novel methods proposed in
Section IV. Our dataset for this experiment is the synthetic tropical fish,
available at https://github.com/siavashk/pycpd.

As in the 3D experiment, we introduce random noise to the data uniformly
distributed over the dataset’s support, $[-2,2]^{2}$. The noise level, denoted
by $\eta$, represents the percentage of noisy data points relative to the
clean dataset and varies within the set $\\{0,10\%,20\%,30\%\\}$.

###### Remark 4

It’s worth noting that the data used in this experiment do not fulfill the
assumption outlined in equation (11) in Section IV. To accentuate the
differences between balanced and unbalanced optimal transport (OT) methods,
noise is only added to point cloud $Y$.

Accuracy.

.

Figure 3: We demonstrate the accuracy comparison of the 6 baseline methods:
CPD[85], 4 nonrigid Wasserstein procrustes[47, 100] methods, TPS-RPM [23] and
an variant introduced in last section, and our four partial OT based methods.
$\eta$ is the noise level. The source and target point clouds are in red and
green colors respectively.

See the figure 3 for the accuracy test. The parameter setting for each method
is similar to the 3D experiments. We find that all methods are capable of
successfully recovering the target point cloud when $\eta=0\%$. When $\eta$
increases to $10\%$, $20\%$, and $30\%$, the four balanced optimal transport
(OT)-based methods—namely OT-RBF, OT-TPS, SOT-RBF, and SOT-TPS—result in
relatively large errors.

Under conditions where the model constraint is not met, TPS-based methods
(including both the old and new variants of TPS-RPM, as well as OPT-TPS and
SOPT-TPS) generally exhibit superior performance. Furthermore, we notice that
Sinkhorn-based methods (old and new TPS-RPM) and the sliced-OPT approach also
perform well.

One plausible explanation for this could be that when the model constraint is
not satisfied, the correspondence between the deformed source points
$\\{f(x^{n})\\}_{n=1}^{N}$ and the target points $\\{y^{m}\\}_{m=1}^{M}$, as
determined by OPT, becomes limited reliable. In such cases, approaches
incorporating soft correspondence, like the Sinkhorn method, or introducing
randomness in the correspondence estimation step, like Sliced-OPT methods,
result in more stable performance.

Since the dataset is small (about $90$ to $100$), each method requires less
than one second, we do not report the wall-clock time in this experiment.

Registration error of all methods in 2D fish experiment
$\eta=0\%$ $\eta=10\%$ $\eta=20\%$ $\eta=30\%$ baselines CPD [85] 0.171 0.171
0.158 0.158 OT-RBF 0.032 0.172 0.285 0.363 OT-TPS 0.031 0.146 0.218 0.278 SOT-
RBF 0.032 0.206 0.383 0.621 SOT-TPS 0.031 0.213 0.337 0.456 TPS-RPM[23] 0.031
0.092 0.110 0.112 TPS-RPM(new) 0.032 0.035 0.037 0.038 ours OPT-RBF $0.032$
$0.194$ $0.173$ $0.192$ OPT-TPS $0.032$ $0.175$ $0.163$ $0.145$ SOPT-RBF
$0.031$ $0.102$ $0.138$ $0.114$ SOPT-TPS $0.031$ $0.032$ $0.031$ $0.033$

TABLE II: This table displays the accuracy performance metrics for six
baseline methods as well as our proposed methods. The variable $\eta$ is the
noise level. For each pair of method and $\eta$, the corresponding cell
contain the the error (59).

##  Appendix E pseudo-code of OPT-TPS, SOPT-RBF, SOPT-TPS

In this section, we present the pseudo-code for the algorithms described in
Section IV, Subsections B, C, and D. The input variables
$X,Y,\zeta,T,\epsilon,C$ are elaborated upon in Subsection E of Section IV. In
Algorithms 2 and 4, $\phi_{T}$ represents the TPS (Thin Plate Spline) kernel
function. The corresponding parameter $D$ is denoted as ‘param’, which can be
manually specified. Algorithms 3 and 4 make use of the variable $t$,
representing the number of projections. In the first line of both methods, the
variable $\lambda$ is initialized based on intuitive understanding, as
outlined in [8, Lemma 3.2]. Generally, it is set to
$K\|\text{mean}(Y)-\text{mean}(X)\|^{2}$, where $K\geq 2$.

Input: $X,Y,\zeta,T,\phi_{T},\text{param},\epsilon$

Output: $B,\beta,\Phi,\alpha$

1 Initilize $R,S,\beta,\alpha,B\leftarrow SR$

2 Initilize $\Phi\leftarrow\phi_{T}(X,X^{T},\text{param})$

3 for  _$T^{\prime}=1,2,\ldots T$_ do

/* Step 1 */

4 $\hat{Y}\leftarrow\Phi\alpha+XB+\beta^{T}1_{N}$

5 run lines 5-8 in Algorithm 1

/* Step 2 */

6 if  _condition for non-rigid is False_ then

7 $\hat{Y}^{\prime}\leftarrow\hat{Y}_{sub}-\Phi_{sub}\alpha$

8 By $(X_{sub},\hat{Y}^{\prime})$, update $R,S,\beta$ via (35),(36),(37),
$B\leftarrow SR$

9 else

10 By $(X,\Phi,\hat{Y})$, update $B,\beta,\alpha$ via (56)

11

Algorithm 2 OPT-TPS

Input: $X,Y,\zeta,T,t,\phi,\text{param},\epsilon,C\leftarrow X$

Output: $R,S,\beta,\Phi,\alpha$

1 Initilize $R,S,\beta,\alpha,\lambda$

2 Initilize $\Phi\leftarrow\phi(X,C^{T},\text{param})$

3 for  _$T^{\prime}=1,2,\ldots T$_ do

/* Step 1 */

4 $\hat{Y}\leftarrow\Phi\alpha+XSR+\beta^{T}1_{N}$

5 Sample $\\{\theta_{t^{\prime}}\\}_{t^{\prime}=1}^{t}\subset\mathbb{S}^{D-1}$

6 $\mathcal{D}\leftarrow\emptyset$

7 for  _t’=1,2,…t_ do

8
$\hat{Y}_{\theta_{t^{\prime}}}=\theta^{T}_{t^{\prime}}Y,Y_{\theta_{t^{\prime}}}=\theta^{T}_{t^{\prime}}Y$

9 Compute optimal tansportation plan $L_{t^{\prime}}$ for
$\text{OPT}_{\lambda}(\hat{Y}_{\theta_{t^{\prime}}},Y_{\theta_{t^{\prime}}})$.

10 $\mathcal{D}_{t^{\prime}}\leftarrow\text{Dom}(L_{t^{\prime}})$

11 update $\hat{y}^{n},\forall n\in\mathcal{D}_{t^{\prime}}$ via (58)

12 If $\\#(\mathcal{D}_{t^{\prime}})<\zeta$, increase $\lambda$; otherwise,
decrease $\lambda$.

13 $\mathcal{D}\leftarrow\mathcal{D}\cup\mathcal{D}_{t^{\prime}}$

14
$(X_{sub},\hat{Y}_{sub},\Phi_{sub})=(X[\mathcal{D},:],\hat{Y}[\mathcal{D},:],\Phi[\mathcal{D},:])$

/* Step 2 */

15 run lines 9-13 in Algorithm 1.

Algorithm 3 SOPT-RBF

Input: $X,Y,\zeta,T,t,\phi_{T},\text{param},\epsilon$

Output: $B,\beta,\Phi,\alpha$

1 Initilize $R,S,\beta,\alpha,\lambda,B\leftarrow SR$

2 Initilize $\Phi\leftarrow\phi_{T}(X,X^{T},\text{param})$

3 for  _$T^{\prime}=1,2,\ldots T$_ do

/* Step 1 */

4 $\hat{Y}\leftarrow\Phi\alpha+XB+\beta^{T}1_{N}$

5 run lines 5-14 in Algorithm 3

/* Step 2 */

6 run lines 6-10 in Algorithm 2.

Algorithm 4 SOPT-TPS

###### Remark 5