Linear optimal transport subspaces for point set classification

The fundamental principle of optimal transport theory relies on quantifying the amount of effort (measured as the product of mass and distance) required to rearrange one distribution to another, which gives rise to the Wasserstein metric between distributions. In the present study, we utilize a linearized version of this metric, as outlined in [23], which is constructed formally through a tangent space approximation of the underlying manifold.

Following the construction in [23], we define the linear optimal transport transform for probability measures in $\mathcal{P}_{2}(\mathbb{R}^{L})$, which is the set of absolutely continuous measures with bounded finite second moments and densities ¹¹1Any $\mu\in\mathcal{P}_{2}(\mathbb{R}^{L})$ has the following two properties (i) bounded second moment, i.e. $\int\|x\|^{2}d\mu(x)<\infty$; (ii) absolute continuity with respect to the Lebesgue measure on $\mathbb{R}^{L}$ with bounded density, i.e., $\mu$ has a density function $f_{\mu}$ defined on $\mathbb{R}^{L}$ with $\|f_{\mu}\|_{\infty}<\infty$. . For simplicity, let us fix a reference measure $\sigma$ as the Lebesgue measure on a convex compact set of $\mathbb{R}^{L}$. Thanks to Brenier’s theorem [24], there is a unique minimizer $T_{\sigma}^{\mu}$ to the following optimal transportation problem

where the push-forward (transport) relation $T_{\sharp}\sigma=\mu$ is defined via $\mu(B)=\sigma(T^{-1}(B))$ for any measurable set $B\subseteq\mathbb{R}^{L}$. The linear optimal transport (LOT) transform is given by the following correspondence

where each probability measure $\mu$ is identified with the optimal transport map $T_{\sigma}^{\mu}:\mathbb{R}^{L}\rightarrow\mathbb{R}^{L}$ from a fixed reference $\sigma$ to $\mu$, which lies in a linear space. This square-root of the minimum is called the Wasserstein-2 distance between $\sigma$ and $\mu$ [25]. The LOT metric between two probability distributions $\mu,\nu\in\mathcal{P}_{2}(\mathbb{R}^{L})$ is ²²2Note that $\|T\|_{\sigma}:=\Big{(}\int_{\mathbb{R}^{L}}\|T(x)\|^{2}d\sigma(x)\Big{)}^{1/2}$.

For simplicity, we denote $\widehat{\mu}$ as the LOT transform of $\mu$, i.e., $\widehat{\mu}=T_{\sigma}^{\mu}$ where $\sigma$ is fixed.

It turned out the linearization ability of LOT is closely related to the scope of the following so-called composition property [26, 27]

where $g\in\mathcal{T}_{L}$, and $\mathcal{T}_{L}$ is the set of all diffeomorphisms from $\mathbb{R}^{L}$ to $\mathbb{R}^{L}$. In particular, given a convex $\mathcal{G}\subseteq\mathcal{T}_{L}$, the LOT embedding of deformed measures via maps in $\mathcal{G}$ become convex ³³3 Note in general $\mathcal{G}_{\sharp}\mu$ is not convex as $(\lambda_{1}g_{1}+\lambda_{2}g_{2})_{\sharp}\mu\neq\lambda_{1}{g_{1}}_{\sharp}\mu+\lambda_{2}{g_{2}}_{\sharp}\mu$. if all $g\in\mathcal{G}$ satisfies the above composition property (4), which is shown more formally below.

When the dimension $L\geq 2$, it is shown in [26] that $g$ can only be “basic” transformations (more specifically, translations or isotropic scalings or their compositions) for the composition property (4) to hold for arbitrary $\mu$’s. Luckily, [27] proposes an approximate composition property for perturbations of the aforementioned basic transformations, the set of which we denote as $\mathcal{A}=\{h(x)=ax+b:a>0,b\in\mathbb{R}^{L}\}$.

With the above approximate composition property, one can show the following approximate convexity analog of Proposition II.1 using Lemma A.3, A.4 of [27]:

For the analysis of discrete point set data, a discrete version of the Linear Optimal Transport (LOT) embedding is required. In this particular case, both the reference $\sigma$ and target $\mu$ are chosen as discrete probability measures, represented by point sets in $\mathbb{R}^{L}$. A point set in a $L$-dimensional space is a finite set of points in $\mathbb{R}^{L}$. A point set $\Omega_{s}$ with $N$ points can be thought as the image of an injective map $s:\{1,\cdots,N\}\rightarrow\mathbb{R}^{L}$ ⁵⁵5Note that a point set may be associated with many injective maps, e.g. the image sets of $s\circ\gamma$ and $s$ are the same for any permutation $\gamma$.. Given a point set $\Omega_{s}$ with $N$ points, we define a discrete probability distribution associated with the point set as

Given a diffeomorphism $g\in\mathcal{T}_{L}$, the push-forward distribution of $P_{s}$ under $g$ is given as

Let $\mathcal{F_{N,L}}$ denote the collection of injective maps from $\{1,\cdots,N\}$ to $\mathbb{R}^{L}$. Given $s,r\in\mathcal{F_{N,L}}$, the optimal transportation (Wasserstein-2) distance between associated distributions $P_{s}$ and $P_{r}$ can be obtained by solving the linear programming problem given below:

where $\pi_{ij}\geq 0$, and $\sum_{i=1}^{N}\pi_{ij}=\sum_{j=1}^{N}\pi_{ij}=1/N$ for all $i,j=1,\cdots,N$. Let us fix some $r\in\mathcal{F_{N,L}}$ and use $P_{r}$ as a reference. It turned out that any minimizer matrix $\pi^{*}$ to the optimal transport problem in (10) is a permutation matrix[25]. In other words, there is a permutation $\sigma_{s}^{*}:\{1,\cdots,N\}\rightarrow\{1,\cdots,N\}$ such that

Hence with $r$ being fixed, an optimal transport map between $P_{r}$ and $P_{s}$ can be determined by $\sigma^{*}_{s}$ and $s$. The LOT transform for $P_{s}$ is defined as [23] ⁶⁶6Note one can write $s\circ\sigma^{*}_{s}=\begin{bmatrix}s(\sigma_{s}^{*}(1)),\cdots,s(\sigma_{s}^{*}(N))\end{bmatrix}^{T}.$ Note also that $\sigma^{*}_{s}$ may not be unique in general, we follow the implementation in [23] to estimate one of them.

and the LOT distance between two point set measures is

where $s,q\in\mathcal{F_{N,L}}$.

Based on the aforementioned theoretical discussions, we put forward a straightforward non-iterative training approach for the classification method. This involves computing a projection matrix that maps each sample in the LOT space onto the subspace $\widehat{\mathbb{V}}^{(k)}$ (as outlined in [21]), generated by the 2$\delta$-convex set $\widehat{\mathbb{S}}^{(k)}$. Specifically, we estimate the projection matrix by applying the following procedure:

Given a set of sample training data, denoted as $\{P_{s_{1}^{(k)}},P_{s_{2}^{(k)}},\cdots\}$, the first step in our proposed method is to apply the LOT transform on them using a reference distribution $P_{r^{(k)}}$. This results in the generation of transformed samples, denoted as $\{\widehat{P}_{s_{1}^{(k)}},\widehat{P}_{s_{2}^{(k)}},\cdots\}$. The reference distribution $P_{r^{(k)}}$ is obtained by selecting a point set at random from the training set, followed by the introduction of random perturbations. Subsequently, we estimate $\widehat{\mathbb{V}}^{(k)}$ as follows:

The proposed method also provides a structure to mathematically encode invariances with respect to deformations that are known to be present in the data [21, 28]. In this paper, we prescribe methods to encode invariances with respect to a set of affine transformations: translation, isotropic and anisotropic scaling, and shear. Detailed descriptions of the deformation types used for encoding invariances and the corresponding methodologies are explained as follows:

1. 1.

   Translation: Let $g(\mathbf{x})=\mathbf{x}+\mathbf{x_{0}}$ be the translation by $\mathbf{x}_{0}=((\mathbf{x}_{0})_{1},(\mathbf{x}_{0})_{2},\cdots,(\mathbf{x}_{0})_{L})\in\mathbb{R}^{L}~{}\mbox{and}~{}P_{s_{g}}=g_{\#}P_{s}$. Using equation (6), we have that $\widehat{P}_{s_{g}}=\widehat{g_{\#}P_{s}}\approx g\circ\widehat{P}_{s}=\widehat{P}_{s}+\mathbf{x}_{0},~{}\mbox{where}~{}\widehat{P}_{s}=((\widehat{P}_{s})_{1},(\widehat{P}_{s})_{2},\cdots,(\widehat{P}_{s})_{L})$. Consequently,

   	$\displaystyle\widehat{P}_{s_{g}}\approx\widehat{P}_{s}+\mathbf{x}_{0}=\widehat{P}_{s}+\left((\mathbf{x}_{0})_{1},(\mathbf{x}_{0})_{2},\cdots,(\mathbf{x}_{0})_{L}\right)$
   	$\displaystyle=\widehat{P}_{s}+(\mathbf{x}_{0})_{1}\left(1,0,0,\cdots\right)+(\mathbf{x}_{0})_{2}\left(0,1,0,\cdots\right)+\cdots$
   	$\displaystyle+(\mathbf{x}_{0})_{L}\left(0,0,\cdots,1\right).$

   Therefore, as in [21, 28], we define the spanning set for translation as

   	$\displaystyle\mathbb{U}_{T}=\{u_{t}(1),u_{t}(2),\cdots,u_{t}(L)\},~{}\mbox{where}$
   	$\displaystyle u_{t}(1)=(1,0,0,\cdots),u_{t}(2)=(0,1,0,\cdots),\cdots,$
   	$\displaystyle u_{t}(L)=(0,0,\cdots,1).$

2. 2.

   Isotropic scaling: Let $g(\mathbf{x})=a\mathbf{x}$ be the normalized isotropic scaling of $P_{s}$ by $a$, where $a\in\mathbb{R}_{+}$ and $P_{s_{g}}=g_{\#}P_{s}$. Using equation (6), we have that $\widehat{P}_{s_{g}}\approx g\circ\widehat{P}_{s}=a\widehat{P}_{s}$. As in [21, 28], an additional spanning set for isotropic scaling is not required as the subspace containing $\widehat{P}_{s}$ naturally contains its scalar multiplication $a\widehat{P}_{s}$. Therefore, the spanning set for isotropic is defined as $\mathbb{U}_{D_{0}}=\varnothing$.

3. 3.

   Anisotropic scaling: Let $g(\mathbf{x})=\breve{\mathcal{D}}\mathbf{x}$ be the normalized anisotropic scaling of $P_{s}$, where $\breve{\mathcal{D}}=\begin{bmatrix}a_{1},&0,&\cdots\\ 0,&a_{2},&\cdots\\ \vdots&\vdots&\ddots\end{bmatrix}$, $a_{i}\neq a_{j}$, $a_{i}\in\mathbb{R}_{+}$, and $P_{s_{g}}=g_{\#}P_{s}$. Using equation (6), we have that $\widehat{P}_{s_{g}}\approx g\circ\widehat{P}_{s}=\breve{\mathcal{D}}\widehat{P}_{s}=\left(a_{1}(\widehat{P}_{s})_{1},a_{2}(\widehat{P}_{s})_{2},\cdots,a_{L}(\widehat{P}_{s})_{L}\right)$. Consequently,

   	$\displaystyle\widehat{P}_{s_{g}}\approx\breve{\mathcal{D}}\widehat{P}_{s}=a_{1}((\widehat{P}_{s})_{1},0,0,\cdots)+a_{2}(0,(\widehat{P}_{s})_{2},0,\cdots)+$
   	$\displaystyle\cdots+a_{L}(0,0,\cdots,(\widehat{P}_{s})_{L}).$

   Therefore, the spanning set for anisotropic scaling is defined as

   	$\displaystyle\mathbb{U}_{D}=\{u_{d}(1),u_{d}(2),\cdots,u_{d}(L)\},~{}\mbox{where}$
   	$\displaystyle u_{d}(1)=((\widehat{P}_{s})_{1},0,0,\cdots),u_{d}(2)=(0,(\widehat{P}_{s})_{2},0,\cdots),$
   	$\displaystyle\cdots,u_{d}(L)=(0,0,\cdots,(\widehat{P}_{s})_{L}).$

4. 4.

   Shear: Let $g(\mathbf{x})=\mathcal{H}\mathbf{x}$ be the normalized shear of $P_{s}$, where $\mathcal{H}=\begin{bmatrix}1,&k_{12},&\cdots\\ k_{21},&1,&\cdots\\ \vdots&\vdots&\ddots\end{bmatrix}$ and $P_{s_{g}}=g_{\#}P_{s}$. Here, the shear matrix $\mathcal{H}$ is constructed using the shear factors $k_{ij}\in\mathbb{R}$, which are located at the non-diagonal positions of $\mathcal{H}$. Using equation (6), we have that $\widehat{P}_{s_{g}}\approx g\circ\widehat{P}_{s}=\mathcal{H}\widehat{P}_{s}=((\widehat{P}_{s})_{1}+k_{12}(\widehat{P}_{s})_{2}+k_{13}(\widehat{P}_{s})_{3}+\cdots,(\widehat{P}_{s})_{2}+k_{21}(\widehat{P}_{s})_{1}+k_{23}(\widehat{P}_{s})_{3}+\cdots,\cdots,(\widehat{P}_{s})_{L}+k_{L1}(\widehat{P}_{s})_{1}+k_{L2}(\widehat{P}_{s})_{2}+\cdots)$. Consequently,

   	$\displaystyle\widehat{P}_{s_{g}}\approx\mathcal{H}\widehat{P}_{s}=\widehat{P}_{s}+k_{12}((\widehat{P}_{s})_{2},0,0,\cdots)+$
   	$\displaystyle k_{13}((\widehat{P}_{s})_{3},0,0,\cdots)+k_{14}((\widehat{P}_{s})_{4},0,0,\cdots)+\cdots+$
   	$\displaystyle k_{21}(0,(\widehat{P}_{s})_{1},0,\cdots)+k_{23}(0,(\widehat{P}_{s})_{3},0,\cdots)+\cdots+$
   	$\displaystyle k_{L1}(0,0,\cdots,(\widehat{P}_{s})_{1})+k_{L2}(0,0,\cdots,(\widehat{P}_{s})_{2})+\cdots.$

   Therefore, the spanning set for shear is defined as

   	$\displaystyle\mathbb{U}_{S}=\{u_{s}(1,2),u_{s}(1,3),\cdots,u_{s}(L,L-1)\},~{}\mbox{where}$
   	$\displaystyle u_{s}(1,2)=((\widehat{P}_{s})_{2},0,0,\cdots),u_{s}(1,3)=((\widehat{P}_{s})_{3},0,0,\cdots),$
   	$\displaystyle\cdots,u_{s}(L,L-1)=(0,0,\cdots,(\widehat{P}_{s})_{L-1}).$

Finally, in light of the preceding discussion, we can approximate the enriched subspace $\widehat{\mathbb{V}}^{(k)}_{E}$ as

where $\mathbb{U}_{A}=\mathbb{U}_{T}\cup\mathbb{U}_{D_{0}}\cup\mathbb{U}_{D}\cup\mathbb{U}_{S}$.

To classify a given test sample $P_{s}$, we first apply the LOT transform to $P_{s}$ to obtain its corresponding LOT space representation $\widehat{P}_{s,r^{(k)}}$ with respect to the reference $P_{r^{(k)}}$ (which was pre-selected duing the training phase). Assuming that the test samples originate from the generative model presented in equation (13) (or equation (14)), we can determine the class of an unknown test sample $P_{s}$ using the following expression:

where $d(\cdot,\cdot)$ is the distance between the test sample and the trained subspaces in the LOT transform space. We can estimate the distance between $\widehat{P}_{s,r^{(k)}}$ and the trained subspaces using $d^{2}\left(\widehat{P}_{s,r^{(k)}},\widehat{\mathbb{V}}^{(k)}_{E}\right)\sim||\widehat{P}_{s,r^{(k)}}-B^{(k)}B^{(k)T}\widehat{P}_{s,r^{(k)}}||^{2}_{L_{2}}$, where the matrix $B^{(k)}$ contains the basis vectors of the subspace $\widehat{\mathbb{V}}^{(k)}_{E}$ arranged in its columns.

Our objective is to analyze how the proposed method performs compared to state-of-the-art approaches in terms of classification accuracy, required training data, and robustness in out-of-distribution scenarios in limited training data setting. To achieve this, we created training sets of varying sizes from the original training set for each dataset under examination. We then trained the models using these training sets and assessed their performance on the original test set. Each train split was generated by randomly selecting (without replacement) samples from the original training set, and we repeated the experiments for each split size ten times. The same train-test data samples were used for all algorithms in each split.

In order to assess the effectiveness of the proposed approach, we utilized several comparison methods. These included PointNet [1], DGCNN [17], and multilayer perceptron (MLP) [32] in FSpool feature embedding space [16]. We also conducted a comparative analysis with various conventional machine learning techniques across different set feature embedding spaces. These included logistic regression (LR), kernel support vector machine (k-SVM), multilayer perceptron (MLP), and nearest subspace (NS) classifier models [32] in GeM1, GeM2, GeM4 [13], COVpool [14, 15], and FSpool [16] embedding spaces. The performance of the proposed method was evaluated in relation to these baselines. We conducted these evaluations in addition to performing out-of-distribution experiments. In the proposed method, we selected the number of basis vectors for the subspaces $\widehat{\mathbb{V}}^{(k)}_{E}$ such that the total variance explained by the chosen basis vectors in the $k$-th class captured up to 99% of the total variance explained by that class.

To assess the relative performance of the methods, we evaluated them on several datasets, including Point cloud MNIST [33, 34], ModelNet [35], and ShapeNet [36] datasets. We additionally applied random translations, anisotropic scaling, and shear transformations to both the training and test sets of the datasets. For the ShapeNet dataset, we tested the methods under two experimental setups: the regular setup, where both the training and test sets contained point sets at the same deformation magnitude level, and the out-of-distribution setup, where the training and test sets contained point sets at different deformation magnitude levels.

We first evaluated the effectiveness of the proposed method by comparing it with other state-of-the-art techniques on two synthetic datasets. The synthetic datasets were generated by selecting one sample per class from the point cloud MNIST and ShapeNet datasets, followed by introducing random translations, anisotropic scaling, and shear transformations to each selected sample to generate training and test sets. Specifically, the training set consisted of two samples per class, while the test set comprised 25 samples per class. The obtained comparative results are displayed in Fig. 1. As observed, the proposed method substantially outperformed the other methods in this synthetic scenario.

We conducted the performance evaluation of the proposed method by comparing it with several state-of-the-art techniques, including PointNet, DGCNN, and MLP in FSpool feature embedding space, on the MNIST, ShapeNet, and ModelNet datasets. Fig. 2 presents the average test accuracy values obtained for different numbers of training samples per class. The results demonstrate that our proposed method outperformed the other methods across the range of training sample sizes used to train the models. Notably, the proposed method’s accuracy vs. training size curves exhibited a smoother trend in most cases compared to the other methods.

To assess the effectiveness of the proposed method under the out-of-distribution setting, we introduced a gap between the magnitudes of deformations in the training and test sets. Specifically, we used $\mathcal{G}_{out}$ as the deformation set for the ‘out-distribution’ test set, while $\mathcal{G}_{in}$ was the deformation set for the ‘in-distribution’ training set. We trained the models using the ‘in-distribution’ data and tested using the ‘out-distribution’ data. For our out-of-distribution experiment, we used the ShapeNet dataset with small deformations as the ‘in-distribution’ training set and the ShapeNet dataset with larger deformations as the ‘out-distribution’ test set (see Fig. 3). The results show that the proposed method outperformed the other methods by an even more significant margin under the challenging out-of-distribution setup, as shown in Fig. 3. Under this setup, the proposed method obtained accuracy figures closer to that in the standard experimental setup (i.e., ShapeNet in Fig. 2). On the other hand, the accuracy of the other methods declined significantly under the out-of-distribution setup compared to the standard experimental setup (see ShapeNet results in Figs. 2 and 3).

We further evaluated the proposed method against various set embedding-based approaches in combination with classical machine learning methods. The study involved comparing the proposed method with different classifier techniques, including LR, k-SVM, MLP, and NS [32], that were employed with various set-to-vector embedding methods, such as GeM (1,2,4) [13], COVpool [14, 15], and FSpool [16]. Fig. 4 illustrates the percentage test accuracy results obtained from these modified experiments, along with the results of the proposed method for comparison. As shown in Fig. 4, the proposed method outperformed all these models in terms of test accuracy.