Year: 2021

Abstract:
A large class of problems in machine learning involve data sets in which each data point is a point cloud in R. Many applications involving point cloud data sets require a means of measuring not only distances, but also angles, volumes, derivatives, and other more advanced geometric notions. Existing approaches for the most part tend to be ad hoc, and lack coordinate-invariance. In this paper we develop a Riemannian geometric structure for point cloud data. By interpreting a point cloud as samples from some underlying probability distribution, the space of point cloud data can be given the structure of a statistical manifold, with the Fisher information metric acting as a natural Riemannian metric. The only requirement on the part of the user is the choice of a meaningful underlying probability distribution, which is more intuitive and natural to make than what is required in existing ad hoc formulations. Two autoencoder case studies involving point cloud data are presented to demonstrate the advantages of our statistical manifold framework: (i) interpolating between two 3D point cloud data sets to smoothly deform one object into another; (ii) finding an optimal set of latent space coordinates that minimizes distortion.

heading: 1 INTRODUCTION
text: A growing number of machine learning problems involve data sets in which each data point is a point cloud in RD, e.g., 3D point cloud obtained by depth cameras. Typical applications include measuring the degree of similarity betweeen two point clouds – the point clouds may be measurements obtained from a depth camera, for example – for which a distance metric on the space of point clouds is needed. Some widely used distance metrics used in this context include the Hausdorff distance (both the original and averaged versions), the chamfer distance (Hausdorff, 1914), and the earth mover distance (Rubner et al., 2000).
However, the distance metric measures just one aspect of point cloud data; other applications may require more advanced concepts and tools. For example, in the case of a moving point cloud, one may wish to measure some of its more dynamic aspects like its velocity or other quantities that require a notion of higher-order derivatives. Applications such as Monte Carlo sampling may require the construction of an isotropic Gaussian distribution on the underlying space, in which case the notion of “angle” is needed in addition to distance (more technically, an inner product on the tangent space is required (Girolami & Calderhead, 2011; Gemici et al., 2016; Mallasto & Feragen, 2018).
Related to the above, the idea of interpreting a point cloud as samples from some underlying probability distribution is well-known, and has been applied to problems ranging from point set registration (Jian & Vemuri, 2005; Wang et al., 2006; Myronenko & Song, 2010; Hasanbelliu et al., 2014; Zhou et al., 2014; Min et al., 2018; Li et al., 2021) to point cloud de-noising (Zaman et al., 2017; Luo & Hu, 2021). In these approaches, divergence measures from probability and informa-
tion theory have been utilized to compute the similarity between point clouds. While effective for certain applications, these divergence measures still only capture just one aspect of point cloud data, and cannot be used to measure other quantities of a more geometric nature.
For more advanced applications, a rigorous mathematical characterization of the space of point cloud data is an essential ingredient to a more comprehensive, robust, and correct (in the sense of being coordinate-invariant and geometrically well-defined) analysis of the types described above, particularly one based on Riemannian geometry. As such, the first contribution of this paper is a Riemannian geometric structure for the space of point cloud data.
The key idea behind our approach draws upon the information geometry framework of (Amari & Nagaoka, 2000; Amari, 2016), in which the space of probability density functions is given the structure of a Riemannian manifold (the statistical manifold), with the Fisher information acting as a natural Riemannian metric (the Information Riemannian metric, or info-Riemannian metric for short). The connection between point cloud data and information geometry is established by constructing a 1-1 mapping from the space of point cloud data to the space of probability density functions, i.e., a point cloud X = {x1, ..., xn | xi ∈ RD} is mapped to a density function p(x;X) on RD in a 1-1 fashion as illustrated in Figure 1.
Two case studies involving autoencoders are presented to demonstrate the benefits of our approach. In the first case study, a pre-trained autoencoder is used to encode two 3D point clouds – one representing a cylinder, one a cone – and the minimal geodesic (or path of shortest length) with respect to the info-Riemannian metric is then constructed between these two objects. The shape evolution obtained for the info-Riemannian metric is seen to be far more natural and intuitive than that obtained for the straight line interpolant in latent space.
In the second case study, we use the info-Riemannian statistical manifold framework to find a set of distortion minimizing latent space coordinates, in the sense that (Euclidean) straight lines in the latent space closely approximate minimal geodesics on the statistical manifold. Such a set of coordinates offers a more discriminative representation for the data manifold (Chen et al., 2020) that results in, e.g., higher linear SVM classification accuracy vis-á-vis existing state-of-the art methods. Experiments are carried out with both synthetic and standard benchmark datasets (ShapeNet (Chang et al., 2015), ModelNet (Wu et al., 2015)).

heading: 2 STATISTICAL MANIFOLDS AND THE FISHER INFORMATION METRIC
text: We begin by extending the original definition of a statistical manifold as follows:
Definition 1. Given an m-dimensional topological manifold1 Θ and a 1-1 map from Θ to the space of probability density functions θ 7→ p(x; θ), the image of this mapping, denoted S := {p(x; θ)|θ ∈ Θ}, is an m-dimensional statistical manifold.
In the original definition Θ is taken to be an open subset of Rm (Amari & Nagaoka, 2000; Amari, 2016)).
By endowing S with a Riemannian metric, S can be given the structure of a Riemannian manifold, allowing for lengths, angles, and volumes to be defined on S in a coordinate-invariant manner. The Fisher information metric serves as a natural Riemannian metric on S: the elements (gij) of the Fisher information metric G(θ) ∈ Rm×m can be expressed as
gij(θ) :=
∫ p(x; θ) ∂ log p(x; θ)
∂θi
∂ log p(x; θ)
∂θj dx, i, j = 1, . . . ,m, (1)
where θ = (θ1, ..., θm) are local coordinates on S. Defining infinitesimal length on S by ds2 = dθTG(θ)dθ, the length of a curve θ(t) on S can then be computed as the integral ∫ T 0 ds. Further details on statistical manifolds and the Fisher information metric can be found in, e.g., (Amari & Nagaoka, 2000; Efron & Hinkley, 1978; Rissanen, 1996; Han & Park, 2014).
1A topological manifold is a locally Euclidean Hausdorff topological space.

heading: 3 STATISTICAL MANIFOLD FRAMEWORK FOR POINT CLOUD DATA
text: With the above statistical manifold preliminaries, we now construct a Riemannian geometric structure for the space of point cloud data. Section 3.1 defines a statistical manifold from the point cloud data, while Section 3.2 uses the Fisher information metric to construct a Riemannian metric for point cloud data. To keep the definitions and results simple, we shall assume throughout all point cloud data consists of exactly n distinct points in RD, i.e., each point cloud is of the form X = {x1, ..., xn |xi ∈ RD, xi ̸= xj if i ̸= j}. The set of all point clouds is denoted X . Later we discuss methods for dealing with point clouds that violate our assumptions. The proofs of propositions in this section are in Appendix B.

heading: 3.1 STATISTICAL MANIFOLD OF POINT CLOUD DATA
text: The core idea for constructing the statistical manifold is to interpret a point cloud X as a set of n samples drawn from some underlying probability density function. Using a kernel density estimator (Parzen, 1962; Davis et al., 2011), a parametric probability density function p(x;X) can be defined in which X itself is the parameter:
Definition 2. Given a positive kernel function K : RD → R such that ∫ RD K(u) du = 1 and a D ×D symmetric positive-definite matrix Σ (the bandwidth matrix), the kernel density estimate
p(x;X) := 1 n √ |Σ| n∑ i=1 K(Σ− 1 2 (x− xi)) (2)
is said to be a statistical representation of the point cloud X ∈ X . The set of statistical representations is denoted S := {p(x;X) |X ∈ X}.
To ensure that S is a statistical manifold, recall from Definition 1 that the following two conditions need to be satisfied: (i) X is a topological manifold; (ii) A 1-1 mapping h : X → S, X 7→ p(x;X) must be defined. The first condition can be satisfied with the “distinct points” assumption:
Proposition 1 (Corollary 2.2.11. in (Knudsen, 2018)). The set of point clouds in which each point cloud is a set of n distinct points of dimension D, is an nD-dimensional topological manifold.
To satisfy the second condition, additional assumptions are needed. The following proposition provides a sufficient condition for h to be 1-1:
Proposition 2. If the set of functions {K(Σ− 12 (x− xi))|xi ∈ F} are linearly independent2 for any arbitrary finite subset F ⊂ RD with |F| ≤ 2n, the mapping h : X → S is 1-1.
With the above proposition, any kernel function that satisfies the linear independence condition is sufficient to ensure the existence of a 1-1 mapping h. For our purposes the standard and widely used normal kernel function satisfies the linear independence condition.
Proposition 3. Under the distinct points assumption, the standard (multivariate) normal kernel function
K(u) = 1√
(2π)D exp(−u
Tu 2 ) (3)
with the scaled identity bandwidth matrix, i.e., Σ = σ2I , satisfies the linear independence condition in Proposition 2.
From Propositions 2 and 3 we have established that, under the distinct points assumption and using the standard normal kernel function, the mapping h : X → S is 1-1; S can therefore be given the structure of statistical manifold. While other choices of kernel function are possible, throughout the remainder of the paper we use the standard normal kernel function. Figure 2 illustrates statistical manifold representations of some example point clouds.
2The linear independence of a set of functions implies that only a trivial linear combination of the functions equals the zero function.

heading: 3.2 INFORMATION RIEMANNIAN METRIC FOR POINT CLOUD DATA SPACE
text: We now equip the point cloud statistical manifold S with the Fisher information metric, which we refer to as the info-Riemannian metric and denote by H. The first task is to define a local coordinate system on the space of point clouds X . Toward this end, we use the matrix representation X ∈ Rn×D of a point cloud X. Observe that the matrix representation is not unique; given an n×n permutation matrix P ∈ Rn×n, then X and PX represent the same point cloud X. Fortunately, this does not cause problems since p(x;X) is defined in a permutation-invariant way, i.e., p(x;X) = p(x;PX) for any n×n permutation matrix P . Throughout we use italics to denote local coordinate representations, e.g., X has local coordinates X ∈ Rn×D, the tangent vector V ∈ TXX has local coordinates V ∈ Rn×D. The info-Riemannian metric H can be expressed in local coordinates coordinates X as follows:
Hijkl(X) :=
∫ p(x;X) ∂ log p(x;X)
∂Xij ∂ log p(x;X) ∂Xkl dx, (4)
for i, k = 1, ..., n and j, l = 1, ..., D. Given two tangent vectors V,W ∈ TXX with respective matrix representations V,W ∈ Rn×D, their inner product is then computed as follows:
⟨V,W⟩X := n∑
i,k=1 D∑ j,l=1 Hijkl(X)V ijW kl. (5)
We note that the coordinate expression of the info-Riemannian metric Hijkl(X) results in a permutation-invariant inner product, i.e., ∑ Hijkl(X)V ijW kl = ∑ Hijkl(PX)(PV ) ij(PW )kl for any n× n permutation matrix P , showing that the metric is geometrically well-defined. Using the standard normal kernel function, the coordinate expression of the info-Riemannian metric Hijkl has a simple analytic expression as follows:
Proposition 4. With the standard (multivariate) normal kernel function and the bandwidth parameter σ, the information Riemannian metric is given by
Hijkl(X) =
∫ p(x;X) K(x−xiσ )K( x−xk σ )
( ∑n
m=1 K( x−xm σ )) 2
[ (x− xi)(x− xk)T σ4 ] jl dx, (6)
Figure 3 shows that, given two moving point cloud data whose velocity matrices have equal Euclidean norm (i.e., ||V||2 = ∑n i=1 ∑D j=1 V
ijV ij), the velocity norms under the info-Riemannian metric are significantly different: the velocity A has a value of 0.2626, while the velocity B has a value of 2.2 × 10−8. In particular, observe that the tangential velocity in the case B, which does not change the overall distribution of the point cloud, has a very small velocity norm under the info-Riemannian metric as it should.
As another illustrative example that highlights the difference between the info-Riemannian metric and Euclidean metric (i.e., ⟨V,W⟩ := ∑n i=1 ∑D j=1 V
ijW ij), Figure 4 shows random walks for point cloud data under these metrics. Consider a randomly sampled and normalized velocity matrix under Euclidean metric. 3D velocity vectors in the sampled velocity matrix are equally likely to point any direction. On the other hand, 3D velocity vectors in the sampled velocity matrix under info-Riemannian metric are more likely to point tangential directions (e.g., more likely the case B than the case A in Figure 3). As a result, unlike the standard Euclidean metric (top), the infoRiemannian metric (bottom) produces random walks that stay close to the initial sphere without significant changes in its overall distribution pattern.

heading: 4 APPLICATIONS TO POINT CLOUD AUTOENCODERS
text: Riemannian geometric formulations of autoencoders for representation learning have recently been introduced and extensively studied in (Shao et al., 2018; Arvanitidis et al., 2018; Yang et al., 2018a; Chen et al., 2018; Kalatzis et al., 2020; Arvanitidis et al., 2020; Chen et al., 2020). In these works, the image of the decoder function is viewed as a low-dimensional manifold embedded in the highdimensional data space – we refer to this manifold as the decoded manifold – and a Riemannian metric for the decoded manifold is obtained by projecting the data space Riemannian metric to this manifold. In contrast, this perspective cannot be reasonably extended to existing point cloud autoencoders (e.g., FoldingNet (Yang et al., 2018b), AtlasNet (Groueix et al., 2018), AtlasNetV2 (Deprelle et al., 2019), and TearingNet (Pang et al., 2021)), due to the absence of a geometrically well-formulated Riemannian manifold structure.
In this section, using the info-Riemannian metric, we extend this perspective by defining a Riemannian metric for the decoded manifold of the point cloud autoencoder. With this info-Riemannian metric, we examine two case studies: (i) interpolation between two points of latent space via the minimal geodesic; (ii) learning an optimal set of latent space coordinates that best preserves distances and angles (or intuitively, minimizes distortion).
Consider a point cloud decoder function with the m-dimensional latent space f : Rm → Rn×D, where the output is expressed in terms of the matrix representation. The projection of the infoRiemannian metric on the point cloud statistical manifold to the decoded manifold is then expressed in latent space coordinates z ∈ Rm as follows:
Gab(z) := n∑ i,k=1 D∑ j,l=1 Hijkl(f(z))(Jf ) ij a (z)(Jf ) kl b (z), (7)
where Jf denotes the Jacobian of f , and the indices a, b both range from 1 to m. The latent space is then assigned a Riemannian metric Gab(z); the following two applications rely on G(z) ∈ Rm×m. Geodesic interpolation: The latent space metric G(z) can be used to find the minimal geodesic curve connecting two point clouds (i.e., the shortest length curve in the decoded manifold). Let z1, z2 be the encoded values of the two point clouds in the latent space Rm. In terms of the metric G(z), the geodesic curve connecting these two points can be determined as a solution to the following optimization problem (Do Carmo, 2016):
min z(t) ∫ 1 0 ż(t)TG(z(t))ż(t) dt, (8)
subject to z(0) = z1 and z(1) = z2. Parametrizing z(t) by a cubic spline with fixed boundary points z1, z2 then leads to an unconstrained optimization problem. To avoid excessive memory consumption when computing the objective function and its gradient, instead of the usual Riemann sum approximation of the integral, we interpret the integral as an expectation over the uniform distribution t ∼ U(0, 1) and accordingly use the mini-batch sampling technique. Learning optimal latent space coordinates: The latent space metric G(z) can be used to formulate a regularization term when training an autoencoder to learn an optimal set of latent space coordinates; by “optimal” we mean G(z) = cI for some positive scalar c, so that the decoder preserves distances and angles as much as possible. Recently, a regularization technique for this purpose has been introduced in (Chen et al., 2020). Specifically, the following regularization term is added to the reconstruction loss function: Ez∼P [∥G(z)− cI∥2F ], (9) where ∥ · ∥F is the Frobenius norm, c = Ez∼P [ 1mTr(G(z))], and P is defined via the modified mix-up augmentation, i.e., z ∼ P ⇐⇒ z = αz1 + (1 − α)z2 where z1, z2 are sampled from the set of encoded training data and α ∼ U(−η, 1 + η) for η > 0. For latent spaces whose dimension m is large, in order to avoid the expensive and memory-consuming computation of G(z) ∈ Rm×m, we use the following regularization term in the subsequent experiments:
Ez∼P [Ev∼N (0,I)[∥vTG(z)v − cvT v∥2]], (10) where we use mini-batch sampling to estimate the expectations. This can be done much more efficiently since we need only compute the Jacobian-vector product, i.e., ∑ a(Jf ) ij a v a.

heading: 5 EXPERIMENTAL RESULTS
text: We now verify the effectiveness of the info-Riemannian metric for the two point cloud autoencoder applications described above using both a synthetic 3D basic shape dataset and standard benchmark point cloud datasets.
The synthetic 3D basic shape dataset consists of cylinders, cubes, cones, and ellipsoids with various aspect ratios of the shape parameters (e.g., radius versus height for the cylinder). We sample 512 points from the surface mesh of the shapes using a greedy sample elimination algorithm, so that each sampled point is approximately the same distance from its neighborhood points (Yuksel, 2015). Each point cloud is then normalized so that distance between the two farthest points is one. Further details about the synthetic 3D basic shape dataset generation are provided in Appendix D.
The standard benchmark point cloud dataset consists of ModelNet (Wu et al., 2015) and ShapeNet (Chang et al., 2015), where ModelNet consists of ModelNet10 and ModelNet40, each of which consist of 10 and 40 shape classes, respectively. Each point cloud data has 2048 points; we normalize these into a unit sphere as done in (Yang et al., 2018b).

heading: 5.1 SYNTHETIC 3D BASIC SHAPE DATASET
text: In Section 5.1.1, we use a dataset consisting of cones, cylinders, and ellipsoids, which are split into training/validation/test sets of size 3196/800/804. We then confirm the validity of the proposed metric by comparing the results of several shape interpolation methods in the latent space.
In Section 5.1.2, to study the effects of the regularization term when learning the optimal latent coordinates (with respect to the info-Riemannian metric), we use a dataset consisting of boxes, cones, and ellipsoids divided into training/validation/test sets of size 720/240/240.
We use DGCNN as the encoder (Wang et al., 2019) and a fully-connected neural network as the decoder. The latent space is assumed to be two-dimensional. For the reconstruction loss term, we use the Chamfer distance. The regularization term in Equation (9) is multiplied by a coefficient λ > 0 and added to the reconstruction loss term. Further details about the network architectures and training are provided in Appendix D.

heading: 5.1.1 EXAMPLE 1: CONE, CYLINDER, AND ELLIPSOID
text: Figure 5 shows the test data encoded in the latent space together with the interpolation trajectories and generated point clouds from those interpolants. In the case of intra-class interpolations (i.e.,
interpolants between cylinders), the linear interpolant clearly passes through the red ellipsoid region in the latent space, with some of the generated point clouds clearly ellipsoids. The geodesic interpolations under the Euclidean and info-Riemannian metrics both avoid ellipsoid regions in the latent space. In particular, the geodesic interpolants between two cylinders are also cylinders; this is wellaligned with human intuition. However, if we look at the generated point clouds in detail, while the info-Riemannian metric produces clearly blue cylinders, some of the generated point clouds with the Euclidean metric are non-blue cylinders (i.e., relatively closer to the ellipsoid region) with noisy side surfaces. For the inter-class interpolations (i.e., interpolants between a cylinder and a cone), the linear interpolant also clearly passes through the red ellipsoid region. The geodesic interpolation under the Euclidean metric produces many non-blue and non-green color shapes during the transition from cylinders to cones, while geodesic interpolation under the info-Riemannian metric produces such cases far less. Overall, it can be observed that the geodesic interpolants under the info-Riemannian metric have minimal shape class changes.
Figure 6 shows the latent spaces produced by the regularized autoencoders using the Euclidean and Info-Riemannian metrics. Compared to the latent space in Figure 5 (trained without regularization), the following observations can be made: (i) the encoded latent spaces of cones and cylinders are flattened, and (ii) the ellipsoids and cones become more discriminative. Furthermore, we note that the encoded latent space curves of cones and cylinders in the right figure (where the infoRiemannian metric is used) are clearly flatter than those in the left figure (where the Euclidean metric is used). In other words, if we linearly interpolate between
two cylinders or two cones, the interpolants in the right case will most likely remain in the same class, unlike the left case (the generated point clouds from the linear interpolants are provided in Appendix E).

heading: 5.1.2 EXAMPLE 2: BOX, CONE, AND ELLIPSOID
text: Figure 7 shows the test data encoded in the latent space together with the visualization of the Riemannian metric, the fitted Gaussian Mixture Model and its samples, and pairwise Euclidean dis-
tances. From the first column of Figure 7, by comparing the results with and without regularization, the following two key observations can be made: (i) in the vanilla autoencoder case (upper case), the major axes of the gray ellipses are aligned with the decision boundary (i.e., a hypersurface that partitions the different class regions), which implies that shapes of different classes are actually more distant in the learned manifold (under the info-Riemmanian metric) than shown in the latent space, and (ii) in the regularized autoencoder case (lower case), by encouraging the metric to be isotropic (i.e., turning ellipses into circles), the gaps between different class regions are widened. The second column confirms that the components of the GMM are better separated after regularization; each component of the GMM on the regularized autoencoder generates high-quality, even samples from the same class shape as shown in the third column. The heat maps of the pairwise distances in the last column also indicate that shapes in different classes are more distant, and therefore more easily separable in the latent space of the regularized autoencoder.

heading: 5.2 STANDARD BENCHMARK DATA
text: To show that the regularization technique with the info-Riemannian metric can benefit unsupervised representation learning from the perspective of discriminative representation learning, we compare the transfer classification accuracy of ShapeNetCore.v2 to ModelNet following the same experimental procedure outlined in (Yang et al., 2018b). When training autoencoders with ShapeNet, random rotations about an axis parallel to the direction of gravity are applied to each point cloud. We use four different point cloud autoencoders: FcNet and FoldingNet adopted from (Yang et al., 2018b), PointCapsNet adopted from (Zhao et al., 2019), and DGCNN-FcNet using DGCNN (Wang et al., 2019); the latent space is 512-dimensional for all. The four autoencoders are trained with and without regularization. In the former case, the regularization terms of Equation (10) under both the Euclidean and info-Riemannian metrics are used while varying the regularization coefficient λ. We distinguish between regularized autoencoders using the Euclidean and info-Riemannian metrics by an “+E” or “+I” after the network name. After network training is finished, we train linear SVM classifiers with the encoded data for ModelNet10 and ModelNet40. These are split into training/test sets of sizes 3991/909 and 9843/2468, respectively. Further experimental detail are provided in Appendix D.
Table 1 shows a comparison of transfer classification accuracy from ShapeNet to ModelNet10 (MN10) and ModelNet40 (MN40) for various recent state-of-the-art methods. In the left column (Adopted from References), the numbers are adopted from previous papers (the experimental procedures may differ slightly from ours). In the right column (Implemented by Authors), we report the best numbers obtained (adopted from Appendix E). As shown in the right column, regularization using the info-Riemannian metric improves classification accuracy over vanilla autoencoders, with
higher accuracy compared to the Euclidean metric case. Although their performance is not directly comparable due to differences in the experimental procedures, it can be seen that our regularized autoencoders are comparable to the state-of-the-art methods and, at least for our implementation, shows higher classification accuracy than vanilla autoencoders.
We also conduct additional experiments to determine how much more robust the representation obtained with our regularization approach is for noisy point cloud data. We add noise with different levels of standard deviation (1%, 5%, 10%, and 20% of the diagonal length of the point cloud bounding box) to point cloud data (see Appendix D for details). Then FcNet is trained with and without regularization in the same way as above.
Table 2 shows a comparison of transfer classification accuracy in the presence of noise. As the noise level increases, the classification accuracy obviously decreases in both cases, but the reduction is more dramatic for the vanilla autoencoder case. Additional experimental results including semisupervised classification task are included in Appendix E. Overall, it is indeed somewhat surprising that unsupervised classification accuracy can be improved with a simple regularization technique in lieu of a complex neural network architecture or loss function.

heading: 6 DISCUSSION AND CONCLUSION
text: This paper has proposed a new Riemannian geometric structure for the space of point cloud data. We have defined a statistical representation of point cloud data and constructed a statistical manifold in a mathematically rigorous way. Then a natural Riemannian metric – Fisher information metric – is assigned to the point cloud statistical manifold, which provides geometrically well-defined measures needed for applications. We demonstrate its advantages through two applications involving point cloud autoencoders: (i) minimal geodesic interpolants under info-Riemannian metric have minimal shape changes compared to the standard linear interpolants, and (ii) the optimal latent coordinates learned using our method produce more discriminative representation spaces than existing methods.
As a potential issue, the “fixed number of points” assumption used in our construction of the statistical manifold may be violated in real world problems. In such cases, we can easily mitigate this issue by matching the number of points in each point cloud through a simple upsampling/downsampling algorithm. Further, the kernel function used in our current implementations, the standard normal kernel function, may not be an optimal choice. Other choices can be explored to enhance our algorithms as long as the conditions of Proposition 2 are satisfied.

heading: REPRODUCIBILITY STATEMENT
text: We have provided clear explanations of the assumptions and statements for the theoretical results, and included complete proofs of the propositions in Appendix B. Also, we have included experiment settings in detail as much as possible in Appendix C and D such as the number of training/validation/test splits of datasets, data pre-processing methods, neural network architectures, and hyperparameters used in training (e.g., batch size, number of epochs, learning rate, etc). In particular, we have included the details of the synthetic dataset generation in Appendix D.

heading: APPENDIX

heading: A EXISTING GEOMETRIC/STATISTICAL METHODS FOR POINT CLOUD DATA

heading: A.1 GEOMETRIC METHODS
text: The Hausdorff distance measures the distance between two non-empty subsets of a metric space (Hausdorff, 1914; 2008). Given two point clouds X = {x1, ..., xn | xi ∈ RD} and Y = {y1, ..., yn | yi ∈ RD} and metric ∥x− y∥2 in RD, the Hausdorff distance can be computed as follows:
max(max x∈X (min y∈Y ∥x− y∥2),max y∈Y (min x∈X ∥x− y∥2)).
The Hausdorff distance is susceptible to outliers; hence, in practice, the average Hausdorff distance is used more often:
1 |X| ∑ x∈X min y∈Y ∥x− y∥2 + 1 |Y| ∑ y∈Y min x∈X ∥x− y∥2,
where |X| denotes the number of elements in the set X. A slightly modified version of this is often referred to as the Chamfer distance (Yang et al., 2018b). The popular point cloud registration algorithm ICP relies on these classes of metrics (Besl & McKay, 1992).
Another popular similarity measure between two point cloud data is the Earth Mover’s Distance (EMD) (Rubner et al., 2000): ∑
x∈X
min ϕ:X→Y
∥x− ϕ(x)∥2,
where ϕ is a bijective mapping. Although the EMD is computationally more expensive than the above Hausdorff distances, comparing point clouds with optimal matching in EMD provides a more robust and well-behaved similarity measure.
The Chamfer distance and EMD are often used to measure the distances between two point clouds, but they are typically computationally expensive. Recently, the sliced Wasserstein distance and its variants have been proposed to more efficiently measure distances (Nguyen et al., 2021). In another study, each point cloud data is represented as a matrix of pairwise Euclidean distances between all points, and the Frobenius norm of the difference between the two matrices is used as the distance between two point clouds Cosmo et al. (2020).
However, all these distance metrics measure just one aspect of point cloud data; other applications may require more advanced concepts and tools. For example, in the case of a moving point cloud, one may wish to measure some of its more dynamic aspects like its velocity or other quantities that require a notion of higher-order derivatives. Applications such as Monte Carlo sampling may require the construction of an isotropic Gaussian distribution on the underlying space, in which case the notion of “angle” is needed in addition to distance (more technically, an inner product on the tangent space is required (Girolami & Calderhead, 2011; Gemici et al., 2016; Mallasto & Feragen, 2018).
We believe a rigorous mathematical characterization of the space of point cloud data is an essential ingredient to a more comprehensive, correct, and robust analysis of the types described above, particularly one based on Riemannian geometry. In this paper, we rigorously define a Riemannian geometric structure for the space of point cloud data.

heading: A.2 STATISTICAL METHODS
text: Interpreting the point cloud data as a set of samples from some underlying probability distribution is very intuitive and natural, and has been adopted in many previous works (Jian & Vemuri, 2005; Wang et al., 2006; Myronenko & Song, 2010; Hasanbelliu et al., 2014; Zhou et al., 2014; Min et al., 2018; Li et al., 2021). Commonly, a mixture model is used to describe the point cloud, written as follows:
p(x;w, θ) := k∑ i=1 wiϕ(x|θi),
where the wi are weights and ϕ(x|θi) are primitive density functions with parameter θi (e.g., ϕ(x|θi) can be a standard Gaussian with mean θi). Given a point cloud X := {x1, . . . , xn|xi ∈ RD}, the parameters w, θ are either fit with data or specified by the user, then the mixture model is used as a statistical representation of the point cloud data. Besides the most popular choice for ϕ, the Gaussian (Jian & Vemuri, 2005), other choices such as the t-distribution (Zhou et al., 2014) or hybrid model (Min et al., 2018) have been explored.
The main purpose behind using statistical representations in existing works are to use the wellknown information-theoretic divergence measures such as the KL-divergence to compute the similarity between point clouds. However, these divergence measures (which are approximate distance metrics) only measure just one aspect of point cloud just like the distance measures discussed in A.1. Effective mathematical concepts and tools for defining and measuring other important geometric aspects of point cloud are still lacking.

heading: B PROOF OF THE PROPOSITIONS

heading: B.1 PROOF OF PROPOSITION 2
text: Proof. Let’s consider two point clouds {yi}ni=1 and {zi}ni=1. To show that the mapping h : X → S is 1-1 (especially injective), we have to prove the following statement:
p(x; {yi}ni=1) = p(x; {zi}ni=1) =⇒ {yi}ni=1 = {zi}ni=1. (11) The conditional statement can be rewritten as follows:
n∑ i=1 K(Σ− 1 2 (x− yi)) = n∑ i=1 K(Σ− 1 2 (x− zi)). (12)
Let denote B = {yi}ni=1 ∩ {zi}ni=1 and |B| = m, and assume that m < n. Then the above equation is reduced to ∑
y∈{yi}ni=1−B
K(Σ− 1 2 (x− y))− ∑ z∈{zi}ni=1−B K(Σ− 1 2 (x− z)) = 0. (13)
Since the sets {yi}ni=1 −B and {zi}ni=1 −B are disjoint and each set has n−m elements, the LHS has 2(n − m) ≤ 2n terms and the terms are different to each other. Then, by the assumption, the above 2(n−m) terms are linearly independent, so the above equation cannot hold by the definition of the linear independence. In the other words, there is a contradiction and we can find that the assumption m < n is wrong. Therefore, m = n, so {yi}ni=1 = {zi}ni=1.

heading: B.2 PROOF OF PROPOSITION 3
text: Proof. First, let’s think about the 1-dimensional case, i.e., x ∈ R. Consider an arbitrary point set {xi}Ni=1 and the corresponding set of functions {K(σ− 1 2 (x−xi))|i = 1, 2, ..., N}, where the kernel function K is
K(u) = 1√ 2π exp(−u 2 2 ). (14)
Suppose there exist non-zero constants a1, a2, ..., aN such that N∑ i=1 aiK(σ − 12 (x− xi)) = N∑ i=1 1√ 2πσ ai exp(− (x− xi)2 2 ) = 0 ∀x ∈ R. (15)
The above expression can be reduced as follows:
exp(−x 2
2 ) N∑ i=1 ai exp(xix) exp(− x2i 2 ) = 0 ∀x ∈ R. (16)
Let λi = ai exp(−x 2 i
2 ) for i = 1, 2, ..., N , then the above equation is simplified as:
N∑ i=1 λi exp(xix) = 0 ∀x ∈ R. (17)
In the other words, it is enough to show that the functions in the set {exp(xix)|i = 1, 2, ..., N} are linearly independent.
There is a useful tool for proving the linear independence of the functions:
Theorem 1 (Wronskian’s Theorem). Let f1, f2, ..., fn : I → R be n − 1 times differentiable functions on an interval I and the Wronskian W (f1, f2, ..., fn) be a function on I defined by:
W (f1, f2, ..., fn)(x) = ∣∣∣∣∣∣∣∣∣∣∣ f1(x) f2(x) . . . fn(x) f ′1(x) f ′ 2(x) . . . f ′ n(x) f ′′1 (x) f ′′ 2 (x) . . . f ′′ n (x) ... ... . . . ...
f (n−1) 1 (x) . . . . . . f (n−1) n (x)
∣∣∣∣∣∣∣∣∣∣∣ ∀x ∈ R, (18)
where this is the determinant of the square matrix constructed by placing the functions and their derivatives in an appropriate way. Then, if the Wronskian of this set of functions is not identically zero, then the set of functions is linearly independent.
Using the above theorem, we can construct the Wronskian W (x) with the set of functions {exp(xix)|i = 1, 2, ..., N}. Substituting zero for x, then we can obtain:
W (0) = ∣∣∣∣∣∣∣∣∣∣ 1 1 . . . 1 x1 x2 . . . xN x21 x 2 2 . . . x 2 N ... ... . . . ...
xN−11 . . . . . . x N−1 N
∣∣∣∣∣∣∣∣∣∣ (19)
The above can be simplified by the mathematical induction:∣∣∣∣∣∣∣∣∣∣ 1 1 . . . 1 x1 x2 . . . xN x21 x 2 2 . . . x 2 N ... ... . . . ...
xN−11 . . . . . . x N−1 N
∣∣∣∣∣∣∣∣∣∣ = ∣∣∣∣∣∣∣∣∣∣ 1 0 . . . 0 x1 x2 − x1 . . . xN − x1 x21 x 2 2 − x21 . . . x2N − x21 ... ... . . . ...
xN−11 . . . . . . x N−1 N − x N−1 1
∣∣∣∣∣∣∣∣∣∣ (20)
= ∣∣∣∣∣∣∣∣∣∣ 1 0 . . . 0 x1 x2 − x1 . . . xN − x1 x21 x 2 2 − x21 . . . x2N − x21 ... ... . . . ...
xN−11 . . . . . . x N−1 N − x N−1 1
∣∣∣∣∣∣∣∣∣∣ (21)
= (x2 − x1) · · · (xn − x1) ∣∣∣∣∣∣∣∣∣∣ 1 . . . 1 x2 . . . xN x22 . . . x 2 N ... . . .
... xN−22 . . . x N−2 N
∣∣∣∣∣∣∣∣∣∣ (22)
= · · · (23)
= n∏ i,j=1 i>j (xi − xj) (24)
Since all the points {xi}Ni=1 are distinct to each other, the above determinant value must be non-zero. Then by Wronskian’s Theorem, the functions {exp(xix)|i = 1, 2, ..., N} are linearly independent, so are the functions {K(σ− 12 (x− xi))|i = 1, 2, ..., N}. Next, let’s think about the general D-dimensional case with D ≥ 2; x ∈ RD. The Wronskian’s Theorem introduced above cannot be used, instead, we adopt a generalized Wronskian’s Theorem applicable to multivariable functions.
Theorem 2 (Generalized Wronskian’s Theorem). Let f1, f2, ..., fn : RD → R be multivariable functions and the generalized Wronskian gW (f1, f2, ..., fn) as a function defined by:
gW (f1, f2, ..., fn)(x) = ∣∣∣∣∣∣∣∣∣∣ f(x) D1f(x) D2f(x)
... Dn−1f(x)
∣∣∣∣∣∣∣∣∣∣ ∀x ∈ RD, (25)
where f(x) = (f1(x), ..., fn(x)), Djf(x) are row vectors, Dj is any partial derivative of order not greater than j, and all Dj are distinct for j = 1, 2, ..., n − 1. We note that if at least one of the generalized Wronskians of this set of functions is not identically zero, then the set of functions is linearly independent.
Similar to the 1-dimensional case xi ∈ R, consider an arbitrary point set {xi}Ni=1 and the corresponding set of functions {K(σ− 12 (x− xi))|i = 1, 2, ..., N}, where the kernel function K is now:
K(u) = 1√
(2π)D exp(−u
Tu 2 ). (26)
Suppose that there exist non-zero constants b1, b2, ..., bN such that
N∑ i=1 biK(σ −D2 (x− xi)) = N∑ i=1 1√ (2πσ)D bi exp(− (x− xi)T (x− xi) 2 ) = 0 ∀x ∈ RD. (27)
The above expression can be reduced as follows:
exp(−x Tx
2 ) N∑ i=1 bi exp(x T i x) exp(− xTi xi 2 ) = 0 ∀x ∈ RD (28)
Let ωi = bi exp(−x T i xi 2 ) for i = 1, 2, ..., N . Then the final expression of the equation is as follows:
N∑ i=1 ωi exp(x T i x) = 0 ∀x ∈ RD. (29)
In other words, it is enough to show that the functions {exp(xTi x)|i = 1, 2, ..., N} are linearly independent.
To apply the generalized Wronskian theorem, we need to find a proper set of partial derivatives. To find such set, our strategy is to use a coordinate transformation with a rotation matrix R ∈ SO(D), i.e., rotating vectors. To be more specific, we will find R that makes the first components of yi = Rxi for i = 1, 2, ..., N all different and use the transformed coordinates in the following proof. For this purpose, we first show that such rotation matrix R exists.
Intuitively, we can find a vector r1 so that x1, ..., xn are projected to r1 as all distinct points. If R has this r1 vector as the first column vector, then R can make the first components of yi = Rxi for i = 1, 2, ..., N all different. More precisely, the first components of the two points xi, xj are same if r1 is in the (D − 1)-dimensional hyperplane {r1 ∈ RD|rT1 (xi − xj) = 0}. Since there are only N(N − 1)/2 combinations of the pairs of (xi, xj) for i, j = 1, 2, ..., n, we can find a vector r1 such that rT1 (xi − xj) ̸= 0 ∀i ̸= j and i, j = 1, 2, ..., n. (30) We pick a coordinate transform matrix R ∈ SO(D) with this r1 as the first column vector. Then this R can make the first components of yi = Rxi for i = 1, 2, ..., N all different.
Using the R, rewrite the set of functions {exp(xTi x)|i = 1, 2, ..., N} as follows:
{exp(xTi x)|i = 1, 2, ..., N} = {exp(xTi RTRx)|i = 1, 2, ..., N} (31) = {exp(yTi z)|i = 1, 2, ..., N}, (32)
where z = Rx and yi = Rxi for i = 1, 2, ..., N . We again note that the first elements of yi for i = 1, ..., N are all distinct to each other.
To use the Theorem 2 in the transformed coordinates z, we can pick a generalized Wronskian gW (z) with the corresponding partial derivatives:
D1 = ∂
∂z1 , D2 =
∂2
∂z21 , ..., Dn−1 =
∂n−1
∂zn−11 , (33)
where z = (z1, ..., zD). With this generalized Wronskian and the set of functions {exp(yTi z)|i = 1, 2, ..., N}, substituting zero vector for z ∈ RD, we can obtain:
gW (0) = ∣∣∣∣∣∣∣∣∣∣ 1 1 . . . 1 y11 y12 . . . y1N y211 y 2 12 . . . y 2 1N ... ... . . . ...
yN−111 . . . . . . y N−1 1N
∣∣∣∣∣∣∣∣∣∣ , (34)
where y1i is the first element of yi for i = 1, 2, ..., N . Since y1i’s are all distinct to each other, we can obtain, just like the 1-dimensional case:
gW (0) = n∏ i,j=1 i>j (y1i − y1j) ̸= 0. (35)
By generalized Wronskian’s Theorem, the functions {exp(xTi x)|i = 1, 2, ..., N} are linearly independent, so are the functions {K(σ− 12 (x − xi))|i = 1, 2, ..., N}, also in the D-dimensional case.

heading: B.3 PROOF OF PROPOSITION 4
text: Proof. Since ∂log p(x;X)∂Xij = 1 p(x;X) ∂p(x;X) ∂Xij , the Riemannian metric Hijkl in equation (4) is∫
p(x;X) 1
p2(x;X)
∂p(x;X) ∂Xij ∂p(x;X) ∂Xkl dx.
By plugging p(x;X) in equation (3) in ∂p(x;X))∂Xij , we get the following expression:
n √ |Σ| ∂p(x;X)
∂Xij =
∂
∂Xij n∑ a=1 K(Σ−1/2(x− xa))
= ∂
∂Xij K(Σ−1/2(x− xi))
= JK(Σ −1/2(x− xi))Σ−1/2
∂
∂Xij (x− xi)
= JK(Σ −1/2(x− xi))Σ−1/2
∂
∂Xij (x− xi)
= −[JK(Σ−1/2(x− xi))Σ−1/2]j ,
JK : RD → RD is the Jacobian of the kernel function K. This consequently leads to
Hijkl(X) :=
∫ p(x;X) (JK ∣∣ h(x,xi) Σ− 1 2 )j(JK ∣∣ h(x,xk) Σ− 1 2 )l
( ∑n
m=1 K(h(x, xm))) 2
dx,
where h(x, xi) = Σ− 1 2 (x− xi).
If K is the standard normal kernel function, then we get JK(x) = −K(x)xT . By plugging this in the above equation with Σ = σ2I , we get the following:
Hijkl(X) =
∫ p(x;X) (JK ∣∣ h(x,xi) Σ− 1 2 )j(JK ∣∣ h(x,xk) Σ− 1 2 )l
( ∑n
m=1 K(h(x, xm))) 2
dx
= ∫ p(x;X) K(h(x, xi))K(h(x, xk))
( ∑n
m=1 K(h(x, xm))) 2
[ (x− xi)(x− xk)T σ4 ] jl dx.

heading: C RANDOM WALKS OVER THE POINT CLOUD RIEMANNIAN MANIFOLD
text: This section describes details of the random walk experiments over the point cloud statistical manifold equipped with a Riemannian metric Hijkl, i, k = 1, ..., n, j, l = 1, 2, 3. First, an initial shape is defined by slightly deforming the reference sphere with the addition of a small velocity matrix V . Then, sequences of point clouds are generated by repeating the following process: i) given a point cloud matrix X , reorder it as a 3n-dimensional vector and the Riemannian metric at X as 3n × 3n matrix, ii) sample a velocity vector from the Gaussian distribution whose covariance is the inverse of the Riemannian metric, iii) reorder the sampled velocity vector to n×3 matrix representation and divide it by its norm induced by the metric, and iv) add it to X .
D IMPLEMENTATION DETAILS FOR THE EXPERIMENTS

heading: D.1 SYNTHETIC 3D BASIC SHAPE DATASET GENERATION
text: For the synthetic 3D basic shape dataset, we define 5 shape classes consist of cylinder, cone, elliptic cone, ellipsoid, and box. Figure 8 shows the representative shape of the shape corresponding to each class and the shape parameters required to define the shape class. We again note that we sample 512 points from the surface mesh of the shapes using a greedy sample elimination algorithm, and each point cloud is then normalized so that distance between the two farthest points is one. In the process of normalizing the point cloud, the dimension of the shape parameter is reduced by one dimension. As an example, in the case of cylinder, the cylinders are contained in a two-dimensional space defined by a radius dimension r and a height dimension h, but it is reduced to one dimension through normalization. In summary, cylinders and cones with two shape parameters are placed in one-dimensional space, and elliptic cones, ellipsoids, and boxes with three parameters are placed in two-dimensional space.
In Section 5.1.1, we use a dataset consisting of cones, cylinders, and ellipsoids, which are split into training/validation/test sets of size 3196/800/804. The detail ranges of the exact value of the shape parameters are shown in Table 3.
In Section 5.1.2, we use a dataset consisting of boxes, cones, and ellipsoids divided into training/validation/test sets of size 720/240/240. The detail ranges of the aspect ratios of the shape parameters are shown in Table 4.

heading: D.2 DETAILS FOR EXPERIMENTS ON SYNTHETIC 3D BASIC SHAPE DATASET
text: We used an encoder with a structure similar to the classification network used in DGCNN (Wang et al., 2019). The input point cloud with dimension 3×512 passes through five EdgeConv layers with point-wise latent space dimensions (64, 64, 128, 256) and a max pooling layer (we do not use a batch normalization layer unlike the original DGCNN classification network), then we can obtain a 1024-dimensional feature vector. Other settings are the same (e.g., k = 20, leaky relu activation) Then this feature vector again passes through three fully-connected neural networks with dimension (512, 256, 2) with leaky relu activation function and linear output activation function; the latent space is two-dimensional. For the decoder model, we simply use a fully-connected neural network as the decoder. The two-dimensional vector on the latent space passes through three fully-connected neural networks with dimension (256, 512, 3×512) with relu activation function and linear output activation function; the output is a 3D point cloud with the number of points 512.
Section 5.1.1: To train the networks, we use ADAM with a learning rate of 0.001 and batch size of 16; the total number of the epochs is 500. The mean value of MEDs of the dataset is 0.0339, and we use the bandwidth value k to 0.5. We use Chamfer distance as the reconstruction loss; for regularization figures, the regularization term with the version of Equation (9) is multiplied by a coefficient λ = 107 with the info-Riemannian metric and by a coefficient λ = 1 with the Euclidean metric and added to the reconstruction loss term for each metric case. The value of η is set to be 0.0. For geodesic computation, we parametrize the curve z(t) by a cubic spline with fixed boundary points z1, z2 and 10 control points. The control points are first initialized with equally spaced linear interpolants between z1 and z2. Then, for each iteration of optimization, we randomly sample 40 points on ti ∼ U(0, 1), i = 1, ..., 40 and calculate an expectation 140 ∑40 i=1 ż(ti)
TG(z(ti))ż(ti) over the sampled points as the approximation of the objective function. We use ADAM with a learning rate of 0.001 and the total number of the iterations is 5000.
Section 5.1.2: To train the networks, we use ADAM with a learning rate of 0.001 and batch size of 16; the total number of the epochs is 3000. The mean value of MEDs of the dataset is 0.0341, and we use the bandwidth value k to 0.5. We use Chamfer distance as the reconstruction loss, and the regularization term with the version of Equation (9) is multiplied by a coefficient λ = 107 and added to the reconstruction loss term. The value of η is set to be 0.0.
Color assigning method: To visually indicate which class generated point cloud belong to, we color these according to the ratio of the Chamfer distances to the nearest point cloud for each class. In detail, the smallest value (distance to nearest point cloud) is found by comparing the distance between the given point cloud and all point clouds of each class in the dataset. Since we are using 3 shape classes in both examples, we call the nearest distance to each class d1, d2, and d3. After
that, the vector d = (d1, d2, d3) is normalized with 2-norm so that the 2-norm of the vector to be 1. Finally, the value 0.2 × Softmax(1/d1, 1/d2, 1/d3) is regarded as the ratio of the distances and a color is assigned to a given point cloud according to this ratio (i.e., linear weighted sum in the RGB coordinate).

heading: D.3 DETAILS FOR EXPERIMENTS ON STANDARD BENCHMARK DATASET
text: We use four different point cloud autoencoders: FcNet, FoldingNet, PointCapsNet, and DGCNNFcNet; the latent space is 512-dimensional. For FcNet and FoldingNet, we use the exactly same point cloud autoencoder structures both adopted from (Yang et al., 2018b). For PointCapsNet, we also use the exactly same point cloud autoencoder structure adopted from (Zhao et al., 2019); we use 16× 32 capsules to restrict the latent space to a reasonable size of 512. For DGCNN-FcNet, we use DGCNN classification network as encoder (i.e., the same encoder architecture used in experiments on synthetic 3D basic shape dataset, see Appendix C.2), and the same decoder structure from FcNet as decoder (i.e., three fully-connected neural networks with dimension (1024, 2048, 3×2048) with relu activation function and linear output activation function). To train the networks, we use ADAM with a learning rate of 0.0001, betas of [0.9, 0.999], and weight decay of 0.000001 and batch size of 16; the total number of the epochs is 500. The mean value of MEDs of the dataset is 0.0356, and we use the bandwidth value k to 0.8. We use Chamfer distance as the reconstruction loss and regularization term with the version of Equation (10) with the value of η to be 0.2. The regularization term is multiplied by various coefficients, where the values of the regularization coefficients are summarized in Appendix D.

heading: D.4 DETAILS FOR EXPERIMENTS ON STANDARD BENCHMARK DATASET WITH NOISE
text: We use the exactly same point cloud autoencoder structures adopted from (Yang et al., 2018b), FcNet, where the latent space is 512-dimensional. We add noise to each point x in point cloud of the dataset (ShapeNet, ModelNet10, and ModelNet40) according to x 7→ x +mv, where v is uniformly sampled on the unit sphere and m is sampled from the Gaussian distribution with zero mean and different levels of standard deviation (1%, 5%, 10%, and 20% of the diagonal length of the point cloud bounding box). The training configuration is the same with the case of Appendix D.3 except the followings. The regularization term is multiplied by λ = 8000. The mean values of MEDs of the dataset are 0.0320, 0.0364, 0.0442, and 0.579 for the cases of the noise levels 1%, 5%, 10%, and 20%, respectively, and we use the bandwidth value k to 0.8.

heading: E ADDITIONAL EXPERIMENTAL RESULTS

heading: E.1 SYNTHETIC DATASET

heading: E.1.1 QUANTITATIVE RESULTS
text: In this section, we quantitatively evaluate how much the regularization term improves class separability. In addition to the example shown in Figure 7, more diverse synthetic datasets are made and experiments are conducted. We use datasets consisting of boxes, elliptic cones, and ellipsoids as the case in Figure 7. But in these additional experiments, for each shape class, we use more diverse parameter configurations. In details, we generate short, normal, and tall shapes for each shape class, and the aspect ratios of the shape parameters are shown in Table 5. We conduct a total of 27 experiments with 33 combinations.
Similar to the experiment of Section 5.1.2, each dataset is divided into training/validation/test sets of size 720/240/240. Training configurations are the same with the experiment of Section 5.1.2, except that the mean values of MEDs are different to each other (but we consistently use the bandwidth value k to 0.5) and the total number of epochs is 500.
In the obtained representation spaces, after fitting the Gaussian mixture model by using the training and validation data, the clustering scores are measured with the test data. The Normalized Mutual Information (NMI), Adjusted Rand Index (ADI), and Silhouette Coefficient (SC) are used as metrics. The results are shown in Table 6. For all three clustering measures, regularized autoencoder shows
higher values, indicating that our regularization technique makes the autoencoder learn more optimal latent spaces.
Some representative examples are shown in Figure 9. The trend of the experimental results is similar to the experimental results in Section 5.1.2. For all five results, the gray ellipses are aligned well with the decision boundary in the vanilla autoencoder (we show the decision boundary in Figure 9 while it is not included in the main manuscript due to lack of space). In the regularized autoencoder, these gray ellipses (or Riemannian metrics) try to become isotropic, so the gaps on the decision boundaries get widened. As a result, different class clusters become farther away from each other. The interpretations about results on GMM and pairwise distances are same as the ones in Section 5.1.2.

heading: E.1.2 LINEAR INTERPOLATIONS USING REGULARIZED AUTOENCODERS
text: The generated point clouds from the representative intra-class linear interpolants between two cylinders and two cones with the regularized autoencoders with the Euclidean metric and infoRiemannian metric are drawn in Figure 10.

heading: E.2 STANDARD BENCHMARK DATASET

heading: E.2.1 PERFORMANCE ANALYSIS WITH VARYING REGULARIZATION COEFFICIENTS
text: Figure 11 shows graphs of the classification accuracy versus reconstruction error for the trained AEs measured on ModelNet datasets, for a range of regularization coefficients. The reconstruction error is measured by the modified Chamfer distance as in (Yang et al., 2018b). Compared to vanilla autoencoders (red), regularized autoencoders under the info-Riemannian metric (blue) show overall higher classification accuracy regardless of the regularization coefficients. At the same time, they do not significantly increase the reconstruction error. On the other hand, when comparing the performance of regularized autoencoders under the Euclidean metric (green), most of these are clearly inferior to the vanilla autoencoder; the others perform even worse. Overall, regularization under the info-Riemannian metric is much more robust to the choice of regularization coefficients compared to using the Euclidean metric.
The linear SVM classification accuracy and reconstruction error (modified Chamfer distance) according to regularization coefficient are shown in Table 9 and Table 10. The tables are also arranged
according to autoencoder models (FcNet vs. Foldingnet vs. PointCapsNet vs. DGCNN-FcNet) and regularization types (Vanilla vs. Euclidean vs. info-Riemannian).
Fi gu
re 10
: T
he ge
ne ra
te d
po in
tc lo
ud s
fr om
th e
lin ea
ri nt
er po
la nt
s of
th e
re gu
la ri
ze d
au to
en co
de rs
w ith
th e
E uc
lid ea
n m
et ri
c (U
pp er
)a nd
in fo
-R ie
m an
ni an
m et
ri c
(L ow
er ).

heading: E.2.2 LEARNING CURVES FOR NOISY POINT CLOUD DATA
text: Figure 12 shows how the linear SVM classification accuracy and reconstruction error (modified Chamfer distance) evolve as the training proceeds, where datasets are ModelNet10 and ModelNet40 and noise levels are 1%, 5%, 10%, and 20% (details about noise are in Appendix D.4). Compared to vanilla autoencoders (light colored lines), regularized autoencoders under Info-Riemannian metrc (dark colored lines) show overall higher classification accuracy, while they show similar levels of reconstruction errors. The increase in classification accuracy becomes more pronounced as the noise level increases. Especially, when the noise level is 10% and 20%, the classification accuracy of the vanilla autoencoder and regularized autoencoder under Euclidean metric gradually decreases as the learning progresses (as the epoch increases). However, such phenomenons do not appear in the
regularized autoencoders under the Info-Rimennian metric. This result implies that our method is very advantageous in situations where there is noise in the data.

heading: E.2.3 EXPERIMENTAL RESULTS ON SEMI-SUPERVISED CLASSIFICATION
text: We train FcNet whose the latent space is 512-dimensional with and without regularization. The training configuration is the same with the case of Appendix D.3 except the following: the regularization term of the regularized autoencoder (i.e., FcNet + I) is multiplied by λ = 8000. In this case, when we training linear SVM classifier, we only use the different numbers of training data (1%, 5%, 10%, and 50% of the overall training data).
Table 7 shows a comparison of transfer classification accuracy according to the percentage of training label used. The overall trend is similar to the results of the experiment with noise in Section 5.2 and Appendix E.2.2. Regularized autoencoders (i.e., FcNet + I) show overall higher classification accuracy compared to vanilla autoencoders (i.e., FcNet), while their reconstruction errors are not significantly different to vanilla autoencoders’. Moreover, the increase in classification accuracy becomes more pronounced as the label rate decreases. In other words, our performance is more effective as the number of labels decreases. Figure 13 shows how the linear SVM classification accuracy and reconstruction error evolve as the training proceeds, where label rate levels are 1%, 5%, 10%, and 50%. Also, similarly, when the label rate is 1%, the classification accuracy of the vanilla autoencoder gradually decreases as the learning progresses (as the epoch increases), but such phenomenons do not appear in the regularized autoencoders. This result implies that our method is also very advantageous in semi-supervised settings.

heading: F COMPUTATIONAL ASPECT
text: First of all, the proposed info-Riemannian metric
Hijkl(X) = ∑ r K(xr−xiσ )K( xr−xk σ ) ( ∑n m=1 K( xr−xm σ )) 2 [ (xr − xi)(xr − xk)T σ4 ] jl
requires to compute all pairwise distances between points in X . This can be computed in parallel using GPU, and it requires a similar level of computation to that of calculating the widely-used average Hausdorff distance metric between point cloud data.
On the other hand, since Hijkl is a huge 4D tensor, there may be a memory issue, but fortunately, there is no need to store this tensor for both autoencoder application tasks studied in this paper: (i) geodesic interpolation and (ii) learning optimal latent space coordinates. Because the calculations commonly required to perform both tasks can be done without storing the tensor Hijkl entirely. Specifically, the commonly-required calculations are as follows: given a neural network decoder f : Rm → Rn×D and a latent value z ∈ Rm with tangent vector v = (v1, . . . , vm) ∈ Rm, we need to compute ∑
i,j,k,l,a,b
Hijkl(f(z))(Jf ) ij a (z)(Jf ) kl b (z)v avb,
where Jf is the Jacobian of f . This can be memory-efficiently done as described in the pytorch-style pseudocode:
1 2 def Info_Riemannian_velocity_norm(f, z, v, sigma): 3 """ 4 input: decoder function ’f’ # torch.nn.Module, output size: 5 (batch_size, num_pnts, x_dim) 6 latent value ’z’ # (batch_size, latent_dim) 7 tangent vector ’v’ # (batch_size, latent_dim) 8 bandwith ’sigma’ # scalar 9 output: velocity square ’v_sq’ # (batch_size)
10 """ 11 # Jacobian vector product 12 X, Jv = torch.autograd.functional.jvp(f, z, v) 13 14 # calculate Info-Riemannian metric 15 batch_size, num_pts, x_dim = X.size() 16 delta_X = (X.unsqueeze(2) - X.unsqueeze(1)) / sigma 17 K_delta_X = Kernel(delta_X.view(-1, x_dim)) 18 K_delta_X = K_delta_X.view(batch_size, num_pts, num_pts) 19 K_delta_X = K_delta_X / K_delta_X.sum(dim=2).unsqueeze(2) 20 21 # summation 22 v_sq = torch.einsum(’nxi, nxij, nij -> nx’, K_delta_X, delta_X, Jv) 23 v_sq = torch.einsum(’nx, nx -> n’, v_sq, v_sq) 24 return v_sq / (sigma ** 2) 25 26 def Kernel(u): 27 """ 28 input: vector ’u’ # (batch_size, vec_dim) 29 output: kernel values ’k_u’ # (batch_size) 30 """ 31 vec_dim = u.size(1) 32 k_u = 1/(2*pi)**(vec_dim/2) * torch.exp(-torch.norm(u, dim=1)**2/2) 33 return k_u
Listing 1: Pytorch-style pseudocode for the calculation of velocity square under Info-Riemannian metric
For training the regularized autoencoders, it is not only required to compute the above quantity but also back-propagate through it. In the following, we report how costly adding the regularization term is, compared to the vanilla autoencoder case. We measure the ellapsed time per one training epoch under GeForce RTX 3090. For representative examples, we use FcNet, FcNet + E, and FcNet + I, and the ellapsed tims for each method is shown in Table 8.

heading: G DIFFERENT KERNEL CHOICES FOR DENSITY ESTIMATIONS
text: In this section, we qualitatively analyze how the statistical representation can change when a kernel function other than the Gaussian kernel function is used. As representative ones, we compare the Gaussian, Uniform, Triangular, Epanechnikov, and Sigmoid kernel functions. All but Gaussian and Sigmoid kernel functions have the finite support. For visualization, comparison is made using 2d point clouds, and data in the shape of circles, ellipses, squares, stars, and hearts are used.
Figure 14 visualizes the probability distributions. Basically, the estimated probability distributions are not that different. However, when using kernel functions with finite support, estimated distributions are not very smooth compared to using a Gaussian or Sigmoid kernel function that have infinite support. Overall, it is difficult to say definitively which choice is better, but we choose the one that produces smooth estimations, the Gaussian kernel.

Review summaries:
This paper proposes to represent point cloud data as samples drawn from statistical distribution, and project the data to an underlying statistical manifold on which various analyses can be conducted in a manner that is more favorable than doing so in a regular Euclidean space.   The authors thoroughly describe how points are converted into samples from a probability density function, and explain the metric that determines the coordinate system on the statistical manifold, which describe the changes to the distribution of point clouds.  To introduce this metric to represent point clouds, the authors propose to use this to find minimal geodesic curve between two data in a latent space. They also propose to add a novel constraint to the error metric of autoencoders, so that the encoded latent space possesses the Fisher information metric.  The effect of this proposal is firstly analyzed by two synthetic experiments. First, the authors show that using the proposal, they are able to find the geodesic curve along the manifold that consists of data with the same label. Then, they demonstrate that the latent space of the constrained autoencoder successfully captures the Riemann metric, as samples of many classes are placed along a single line, or are cleanly separated.  Finally, the authors transfer the learnt representation between dataset to demonstrate that the proposed metric space is suitable for representation of various point cloud data.

The paper concerns data sets of point clouds in Euclidean space and proposes to analyze those as samples from underlying probability distributions. Each point cloud being represented by a probability distribution can now be seen as a point in a statistical manifold on which the Fisher information metric provides a natural Riemannian structure. The authors presents two experiments in this setup where they train autoencoders to represent point clouds, interpolate between point clouds using geodesics, and where the latent representation is trained so that straight line are close to geodesics.

The paper studies 3D point cloud data. The core message is that one can regard a point cloud as a collection  of samples from an underlying distribution and, hence, one can imagine a manifold of point clouds and construct an associated distance metric via the Fisher information. That distance can then be used for certain manipulation and analysis tasks, e.g., to interpolate ("morph") between point clouds along the manifold, or to regularise representation learning. A central assumption of the paper is that it is easier to choose a meaningful probability distribution than to make other a-priori assumptions "ad hoc" about the point cloud. However it is not discussed how to actually choose a suitable distribution, the paper limits itself to Gaussian kernel density (seemingly also "ad hoc" for mathematical convenience).