Exploring the Common Principal Subspace of Deep Features in Neural Networks
Haoran LiuHaoyi XiongYaqing WangHaozhe AnDongrui WuDejing Dou
The first author is affiliated to Texas A&M University, College Station, TX at present. This work was done when the first author was an intern under supervision from the second author at Baidu Research. Corresponding Author: Haoyi Xiong<EMAIL_ADDRESS>
H. Liu, H. Xiong, Y. Wang, H. An, D. Wu, and D. Dou

H. LiuH. Xiong Y. Wang D. Dou Big Data Lab, Baidu Research

H. AnDepartment of Computer Science,

University of Maryland, College Park, MD
D. Wu School of Artificial Intelligence and Automation

Huazhong University of Science and Technology

We find that different Deep Neural Networks (DNNs) trained with the same dataset share a common principal subspace in latent spaces, no matter in which architectures (e.g., Convolutional Neural Networks (CNNs), Multi-Layer Preceptors (MLPs) and Autoencoders (AEs)) the DNNs were built or even whether labels have been used in training (e.g., supervised, unsupervised, and self-supervised learning). Specifically, we design a new metric $\mathcal{P}$-vector to represent the principal subspace of deep features learned in a DNN, and propose to measure angles between the principal subspaces using $\mathcal{P}$-vectors. Small angles (with cosine close to $1.0$) have been found in the comparisons between any two DNNs trained with different algorithms/architectures. Furthermore, during the training procedure from random scratch, the angle decrease from a larger one ($70^\circ-80^\circ$ usually) to the small one, which coincides the progress of feature space learning from scratch to convergence. Then, we carry out case studies to measure the angle between the $\mathcal{P}$-vector and the principal subspace of training dataset, and connect such angle with generalization performance. Extensive experiments with practically-used Multi-Layer Perceptron (MLPs), AEs and CNNs for classification, image reconstruction, and self-supervised learning tasks on MNIST, CIFAR-10 and CIFAR-100 datasets have been done to support our claims with solid evidences.

Keywords: Interpretability of Deep Learning, Feature Learning, and Subspaces of Deep Features

§ INTRODUCTION

Blessed by the capacities of feature learning, deep neural networks [1] have been widely used to perform learning in various learning settings (e.g., supervised, unsupervised, and self-supervised learning), ranging from classification, to generation [2, 3]. To better euclid the features learned by deep models, numerous works have studied on interpreting the features spaces of the well-trained models as follows.

Backgrounds and Related Works.
As early as 2013, [4] proposed to visualize the features learned by deep convolutional neural networks (CNN) and made sense of discriminative learning via deep feature extraction.
For generative models, [5, 6] studied the interpolation of latent spaces while [6] discovered an user-controlled way to manipulate the images generated through the surrogation of latent spaces via manifolds.

Later, [7] presented the visual concepts learned in the feature spaces of discriminative models through network dissection on specific datasets while the same group of researchers also proposed GAN dissection [8] – an interactive way to manipulate the semantics and style of image synthesis.

[9] compared GAN and Gaussian Mixture Models (GMMs) to understand the capacity of distribution learning in GAN.  [10] proposed to improve understanding and interpolation of Autoencoders using adversarial regularizer while [11] compared AEs with its variational derivatives to interpret the latent spaces.
More recently, the authors in [12] studied the “steerability” of GAN, where authors discovered point-to-point editing paths for content/style manipulation. [13] uncovered the phenomena that DNN classifiers with piecewise linear activation tend to map the input data to linear subregions. Apparently, many impressive studies are not well discussed here [14, 15, 16].

While existing studies primarily focus on the interpolation of a given model to discover mappings from the feature space to outputs of the model (e.g., classification [7] and generation [12]), the work is so few that compares the feature spaces learned by deep models of varying architectures (e.g., MLP/CNN classifiers versus Autoencoders) for different learning paradigms (e.g., supervised/self-supervised classification [17, 18], unsupervised data reconstruction or de-noising [11], etc.). Furthermore, the dynamics of feature space learning over the training procedure still has not been known yet. In our research, we aim at compare feature spaces of DNN models based on various architectures/paradigms and try to understand how the space evolves during the training process and how feature space learning connects to data distributions and performance of models.

Intuitions and Hypotheses.  Based on the same training dataset, it is not difficult to imagine well-trained DNN classifiers of various architectures in supervised learning settings would map these samples into feature vectors that share certain linear subspaces. To yield an appropriate analysis, let simply push back from the classification results – in any well-trained model, the feature vectors should be capable of being projected to the same set of ground-truth labels (as well-trained models fit training datasets well) through a Fully-Connected (FC) Layer (i.e., a linear transform) and/or a Softmax operator, even though the parameters of networks are different. We further doubt that such subspace might be not only shared by supervised learners but also with Autoencoders (AEs) which are trained to reconstruct input data without any label information in an unsupervised manner, or even shared with self-supervised DNN classifiers (e.g., SimCLR [17]) which train CNN feature extractor and the classifiers separately in an ad-hoc manner.

We specifically hypothesize that (H.I:) there exists certain common subspace shared by the feature spaces of well-trained deep models using the same training datasets, even though the architectures (MLPs, CNNs, and AEs) and the learning paradigms (supervised, unsupervised, and self-supervised) are significantly different. As the training procedure usually initializes the DNN models from random weights, we further hypothesize that (H.II:) there might exist a process of convergence in feature space learning over training iterations, where feature spaces in the early stage of training procedure would vary from each other for different random initializations however they would gradually evolve to shape the common subspace for convergence. Finally, we hypothesize that (H.III:) the convergence to the shared common subspace would connect to the data distribution and performance of models, as such behavior indicates how well the features are learned from data.

Results and Contributions.

To test above three hypotheses, this work makes contributions in proposing new measures to the DNN features and conducting extensive experiments for empirical studies. Given a set of feature vectors extracted from a DNN using either the training or testing datasets, we form a #samples $\times$ #features matrix[#samples and #features refer to the numbers of samples and features respectively.], then we use the top left singular vector of such matrix as the vector representing the the principal subspace [19]  (for more about feature extraction from different architectures and computation, please refer to Section 2).
For convenient, we name the principal subspace as the $\mathcal{P}$-vector of the model based on training/testing sets accordingly.

We train deep models using various DNN architectures, multiple learning paradigms, and datasets, with the checkpoint restored per epoch. Then, we use these models to discover some interesting phenomenons. Results and evidences to support the three hypotheses are summarized as follows.

(1) Small Angles between Principal Subspaces of Deep Features Learned by Various DNNs (H.I).

Given multiple DNNs trained using the same dataset, we estimate their $\mathcal{P}$-vectors and calculate the angles between their $\mathcal{P}$-vectors to measure the similarity in their deep feature spaces. Our experiments find that the angles between $\mathcal{P}$-vectors of these DNNs are small with a high cosine similarity close to $1$, no matter how DNNs are trained or modeled. For example, Figure 1 (a)–(b). show the high cosine measure (close to $1$) for angles between $\mathcal{P}$-vectors of any two well-trained DNNs (CIFAR-10) for eight different architectures/tasks. The experiment results confirm the existence of common subspaces in deep features learned from the same training set, which is invariant under various architectures/tasks.

[Cosine (Train)]

[Cosine (Test)]

Common Principal Subspaces. We present cosine (in the range of [0,1]) of angles between principal subspaces of well-trained models in various architectures under different learning paradigms, where angles are measured as the angle between $\mathcal{P}$-vectors based on training and testing datasets respectively. A well-trained model here is the one trained under the suggested settings after 200 epochs for supervised/self-supervised CNN classifiers and 100 epochs for unsupervised AEs. Experiments are carried out in 5 independent trials (5 random seeds) with averaged results reported.

(2) Converging Trends of Angles over the Number of Training Epochs (H.II). For any model in the training procedure, we checkout their checkpoints restored after every epoch and extract a $\mathcal{P}$-vector for every checkpoint. Our experiments find that, the principal subspace of deep features at beginning of a training procedure is very different with the well-trained one, but it would slowly converge to the well-trained one and get closer and closer to the common subspace over epochs. For example, Figure <ref> (a)–(f) demonstrate the angles between the $\mathcal{P}$-vectors of the well-trained model and the checkpoint per epoch, where consistent converging trends from large angles to small ones (e.g., $\approx10^\circ$ for supervised CNN classifiers and SimCLR; $\leq 10^\circ$ for ConvAEs/DenoiseAEs; and $\approx30^\circ$ for SupCon [18]) could be found in the comparisons among eight different architectures/tasks.

[Supervised (Train)]
[Supervised (Test)]

[Unsupervised (Train)]
[Unsupervised (Test)]

[Self-Sup. (Train)]
[Self-Sup. (Test)]

Convergence to the Common Subspace measured using Angles between $\mathcal{P}$-vectors based on CIFAR-10. We present angles between principal subspaces of the well-trained model and the checkpoint per training epoch (from 0$^{th}$–199$^{th}$ epochs) of three learning paradigms.

(3) The Angle between Principal Subspaces of Deep Features and Data Distribution and its Connections to Generalization (H.III). We compare principal subspaces between deep features (i.e., $\mathcal{P}$-vector) and raw data (i.e., the top left singular vector of the #samples $\times$#dimensions data matrix). For convenience, we name them as model $\mathcal{P}$-vector and data $\mathcal{P}$-vector respectively, and measure the angles between data and model $\mathcal{P}$-vectors, as the proxy measurement of angles between the principal subspaces of deep features and data. We also find the converging trends of such angles over training epochs, where we can again see the well-trained models would incorporate smaller angles ($\approx50^\circ$) than ones in the early stage of training processes ($\approx80^\circ$). We correlate such angles with training/testing accuracy of DNNs, where we observe significant negative correlations in most cases of experiments. We further demonstrate the feasibility of using such angles to predict the generalization performance of a DNN model with their training data only.

The most relevant studies to our work are [7, 13, 20, 21]. For discriminative models, [7] recovered visual features learned by CNN classifiers with a priorly labeled dataset, and quantified then compared the feature learning capacities (namely “interpretability” in the work) of different DNN models through patterns matching with the ground truth. Compared to [7], we carry out the empirical studies on a wide range of datasets without any prior information on their features and observe consistent phenomena in the distribution of samples in the feature spaces.

Furthermore, while [13] studied the properties of regions where a supervised DNN classifier with piecewise linear activation behaves linearly, our work observes the common linear subspaces shared by the features learned by the networks that are trained with different architectures (e.g., MLP/CNN classifiers and AEs with ReLU activation) and paradigms (e.g., supervised, unsupervised, and self-supervised learning). Furthermore, both our work and [20] compare feature/representation through SVD, while we perform SVD to investigate the distribution of samples in principal subspace of deep features and [20] uncovered the latent structures in input-and-output. To the best of our knowledge, we make unique contributions compared to the above work.

§ MEASURING PRINCIPAL SUBSPACES OF DEEP FEATURES WITH $\MATHCAL{P}$-VECTORS

§.§ Methodology

As shown in Figure <ref>, we extract deep feature vectors from DNNs and estimate the $\mathcal{P}$-vector in two steps.

Feature Vector Extraction. Given a model, either the well-trained one or a checkpoint obtained during the training process, we extract the feature vector for every sample, with respect to the architectures. For DNN classifiers (either under supervised/self-supervised learning), we use the output of CNN feature extractor (i.e., the input to the Fully-Connected Layer) as the feature vector of the given sample, while we vectorize the output bottleneck layer as the feature vector for AEs. Note that, in our research, we consider AEs with symmetric architectures of encoders and decoders only.

Deep Feature Vectors Extraction and $\mathcal{P}$-vectors Estimation

Singular Value Decomposition for $\mathcal{P}$-Vector estimation. Given the feature vector for every sample, we form a #samples$\times$#features matrix and perform SVD to obtain the top left singular vector as the $\mathcal{P}$-vector. Furthermore, there is no need to solve the singular vectors of the complete spectrum, as only the top singular vector is requested for $\mathcal{P}$-vector estimation. In this way, we propose to use Randomized SVD [22] that compresses the feature domain and approximate the low-rank structure of SVD for acceleration purpose. Actually, we compare the numerical solution of Randomized SVD and Common SVD for ResNet-50 on CIFAR-10 dataset (#features=256 and #samples= 50,000), where we need to perform SVD on a #samples$\times$#features matrix and the $\mathcal{P}$-vector/top left singular vector should be with  dimensions. Compared to the vanilla SVD, around 109x (from 44.88 seconds to 0.41 seconds) speedup has been achieved by Randomized SVD while no significant numerical errors having been found.

§.§ Measuring Angles between $\mathcal{P}$-Vectors

Given the same set of samples, $\mathcal{P}$-vectors of any two DNNs should be in the equal length, as they are both the top singular vectors in the sample side. Thus, we can measure the angle between two $\mathcal{P}$-vectors as a proxy of the angle between two principal subspaces. A larger cosine of the angle between two vectors (e.g., close to 1.0) usually refers to the evidence that the two DNNs share common principal subspace in their features.

Note that, in high-dimensional spaces, the chance of orthogonality between two random vectors appears more frequently, due to the curse of dimensionality [23]. Thus, given a sample set such as CIFAR-10 with 50,000 samples, when the cosine measure close to 1.0 or the angle between the two $\mathcal{P}$-vectors (with 50,000 dimensions) is small, we can conclude that the two networks would share a subspace in the feature spaces in high confidence.

Of-course, there exists other ways estimating angles between principal subspaces [24].

Singular Value Distribution and Explained Variance Ratios of the Matrix of Feature Vectors

§.§ Statistical Properties of $\mathcal{P}$-vectors

A possible threat to validity of using $\mathcal{P}$-vectors to analyze the feature space of DNN is that $\mathcal{P}$-vectors might fail to capture the necessary information to represent the features learned. To validate the relevance of $\mathcal{P}$-vectors, we analyze the statistical properties of $\mathcal{P}$-vectors as follows.

Significance of $\mathcal{P}$-vectors. We carried out Singular Value Decomposition on the matrix of feature vectors, using CIFAR-10 and CIFAR-100 datasets both based on ResNet-50 models, and obtain the distribution of singular values over indices. More specifically, we compute the distributions of singular values for the feature matrices obtained in the 1$^{st}$, 60$^{th}$, 120$^{th}$, 160$^{th}$ and 200$^{th}$ epochs to monitor the change of singular value distributions throughout the training procedure. It has been observed in Figures <ref>(a) and (b) that a “cliff” pattern in the distribution becomes more and more significant after epochs of training for both CIFAR-10 and CIFAR-100 datasets – a very small number (less than 10) of top singular values might dominate the whole distribution. In Figure <ref>(c), we plot the curve of explained variance ratio $\sigma_k^2/{\sum_{j=1}^d\sigma_j^2}$ for every pair of singular vectors, using well-trained models of 200 epochs based on CIFAR-10 and CIFAR-100, where $\sigma_k$ refers to the $k^{th}$ singular value and $d$ is the rank of matrix. The explained variance ratio of the top-1 singular vectors (i.e., the $\mathcal{P}$-vector and the top-1 right singular vector) is more than 50% while the second top singular vectors are less than 10%. Results show the use of $\mathcal{P}$-vectors could represent the principal subspace of deep features learned. In addition to top-1 singular vectors (the $\mathcal{P}$-vector), Appendix (A.7) also presents the results of using top-$2,3,\dots$ singular vectors, no consistent observations have been obtained (not representative).

Distribution of Values in a $\mathcal{P}$-vector. A potential threat to validity in this work is that when $\mathcal{P}$-vectors were trivial with identical values in every dimension, the angles between these vectors would close to zero. To make our analysis more rigorous, we take a closer look at the $\mathcal{P}$-vector. In A.6 of Appendix, we include a study to observe the distribution of values in every dimension of the  $\mathcal{P}$-vector. With $\mathcal{P}$-vectors extracted from ResNet50, We first retrieve the value on every dimension of the vector and count the frequency that each value appears in the vector. In Figures <ref> (a)–(e) and (f)–(j), we plot the histograms characterizing the frequency of values appeared in every dimension of the P-vector with smoothed curves  for CIFAR-10 and CIFAR-100 training datasets respectively, where results of models trained after the 0$^{th}$ (prior to the training procedure), 60$^{th}$, 120$^{th}$, 160$^{th}$, and 200$^{th}$ epoch have been presented. It is obvious that the values in every dimension of the $\mathcal{P}$-vectors are not identical. Though the shapes of histograms are similar, the frequency distributions are quite different from the ranges and values' perspectives. Thus, we can conclude that $\mathcal{P}$-vectors are different for different datasets. In this way, we can conclude that the $\mathcal{P}$-vectors for the same model in various training epochs are different; the $\mathcal{P}$-vectors for various datasets are different; and $\mathcal{P}$-vectors are not trivial while values are not identical.

§ UNCOVERING COMMON PRINCIPAL SUBSPACE USING ANGLES BETWEEN $\MATHCAL{P}$-VECTORS

[MNIST (Train)]
[MNIST (Test)]

[CIFAR-100 (Train)]
[CIFAR-100 (Test)]
Cosine of angles between principal subspaces of deep features, measured using $\mathcal{P}$-vectors, for models trained under default settings.

In this section, we elaborate the phenomena that we mentioned in Figures <ref> and <ref> of Section 1 using rich datasets, and provide details of evidences to support the hypotheses H.I and H.II.

Common Principal Subspaces. To demonstrate the existence of common subspaces, we propose to measure the angles between principal subspaces of deep features learned using the same datasets with various architectures under different tasks. Specifically, as was discussed in Section 3, we use the angle between $\mathcal{P}$-vectors as the proxy of measurements.
Figure <ref> (a)–(b) presents the cosine of angles between the $\mathcal{P}$-vectors of well-trained models (trained with various architectures and tasks using CIFAR-10 dataset).

We carry out the same experiments on CIFAR-10 and CIFAR-100 datasets, and present the cosine of angles between the $\mathcal{P}$-vectors of well-trained models in Figure <ref>. All experiments based on both training and testing datasets demonstrate that, with a relatively large cosine (close to $1.0$), the models well-trained using the same dataset, no matter what types of architectures or whether labels has been used in the feature learning of various tasks, share a common principal subspace. In addition to above results, we also include the experiment results based on ImageNet [25] dataset, which is sufficiently large and further confirm our observations.

Converging Trends. To understand the dynamics of feature learning that shapes the common principal subspace from scratch, we measure the change of angles over epochs between the $\mathcal{P}$-vectors of training models (i.e., checkpoint for every epoch) and well-trained ones in comparisons. In Figure <ref> (a)–(f), for every model of various tasks, we present angles between the $\mathcal{P}$-vectors of the well-trained model and itself's training checkpoint per epoch. A consistent convergence could be observed. As the $\mathcal{P}$-vectors of well-trained models are close to each other (see in Figures <ref> and <ref>), we can conclude the feature vectors extracted by these models would evolve on time and gradually converge to share the common subspace during the learning procedure.

[Sup. vs. Sup.]
[Sup. vs. Unsup.]
[Sup. vs. Self-Sup.]

[Unsup. vs. Sup.]
[Unsup. vs. Unsup.]
[Unsup vs Self-Sup]

[Self. vs Sup.]
[Self. vs Unsup.]
[Self. vs Self-Sup.]
Angles between the principal subspaces of the well-trained model and the checkpoint per training epoch using $\mathcal{P}$-vectors based on CIFAR-10.

We carry out a more comprehensive model-to-model comparison using CIFAR-10 dataset. Figures. <ref> (a)–(c) present the converging trends of $\mathcal{P}$-vector angles between the well-trained supervised models and the checkpoints per epoch of supervised, unsupervised, and self-supervised learning models, where we use the well-trained Wide-ResNet28 (trained with 200 epochs under suggest settings) as the reference of supervised models. Figures. <ref> (d)–(f) present $\mathcal{P}$-vector angles between the well-trained ConvAE (trained with 100 epochs under suggest settings as the reference of unsupervised learning) and the checkpoints per epoch of supervised, unsupervised, and self-supervised learning models. Figures. <ref> (g)–(i) present the $\mathcal{P}$-vector angles between the well-trained SimCLR representations (trained under suggest settings [17] as the reference of self-supervised learning) and the checkpoints per epoch for supervised classification, unsupervised image reconstruction, and self-supervised contrastive learning tasks. The converging trends in all comparison further validate our hypotheses.

Note that, even though a model has been selected as the reference of well-trained models for the comparison in every setting, the angle between $\mathcal{P}$-vectors of the model to itself would not converge to zero in Figure <ref>. As was stated in the caption of Figure <ref>, we carried out experiments with 5 independent trails with different random seeds. Among these models, only one would be selected as the reference. In Figure. <ref>, we plot the average angles for the 5 trails, where the angles might not be able to converge to zero even when the target and the reference are with the same architectures.

[Supervised CIFAR-10]
[Supervised CIFAR-100]

[Unsupervised CIFAR-10]
[Unsupervised CIFAR-100]

[Self-Supervised CIFAR-10]
[Self-Supervised CIFAR-100]
Angles between principal subspaces, measured using $\mathcal{P}$-vectors based on CIFAR-10 and CIFAR-100, between well-trained models and checkpoints per training iteration in the first epoch.

Discussion. To better visualize every comparison, we use the maximal angles achieved during the first epoch to represent the angles between $\mathcal{P}$-vectors corresponding to the first epoch in Figures <ref> and <ref>. Actually, in the first epoch, there would involve some non-monotonic trends for the angles varying over the number of iterations. Figure <ref> presents the angles between the $\mathcal{P}$-vectors of the training and well-trained models over the number of iterations in the first epoch in three settings of learning, where we use the training set of CIFAR-10 and CIFAR-100 for the experiments (Please refer to the results based on the testing set of CIFAR-10 and CIFAR-100 in Appendix A.2). While such “zigzag” curves could be found in the first few training epochs, the angles between the $\mathcal{P}$-vectors drop down and converge to smaller ones over the number of iterations in the rest of learning procedure.

Note that in Appendix (A.7), we include the comparison between the top singular vectors other than the $\mathcal{P}$-vectors, i.e., the angles of top-2, 3, 4, 5 and 6 left singular vectors between models, where we cannot observe the phenomena.

§ PREDICTING GENERALIZATION PERFORMANCE USING $\MATHCAL{P}$-VECTORS

In this section, we study the angles between the principal subspaces of (raw) data and deep features, and use such angles to predict the generalization performance of models (H.III).

§.§ Measuring Angles between Principal Subspace of Raw Data and Deep Features with $\mathcal{P}$-vectors

Given the raw data matrix, i.e., a #samples$\times$#data dimensions matrix, we obtain the Data $\mathcal{P}$-vector[We use the term “model $\mathcal{P}$-vector” to represent a $\mathcal{P}$-vector estimated using feature vectors of a deep model, while using “data $\mathcal{P}$-vector” as the top left singular vector of the raw data matrix.] of these samples using the top left singular vector of the raw data matrix, which represents the principal subspace of raw features (or the position of every sample projected by the principal component of the data). To understand the connection between models and data, We carry out case studies using CIFAR-10 dataset, including (1) measuring the angles between the principal subspace of deep features and the raw data and (2) understanding the changes of angles over training epochs, using model and data $\mathcal{P}$-vectors.

[Supervised CIFAR-10]
[Supervised CIFAR-100]

[Unsupervised CIFAR-10]
[Unsupervised CIFAR-100]

[Self-Supervised CIFAR-10]
[Self-Supervised CIFAR-100]

Angles between principal subspaces of raw data and the checkpoint per training epoch, measured using model and data $\mathcal{P}$-vectors with CIFAR-10 (C10) and CIFAR-100 (C100).

Figure. <ref> shows the converging trend of angles between data and model $\mathcal{P}$-vectors over number of training epochs for supervised, unsupervised, and self-supervised learning models. The angles between the data and model $\mathcal{P}$-vectors start from about orthogonal and generally decrease and converge to about 50$^\circ$ degree to 60$^\circ$ degree. Note that the angle points corresponding to the first epoch on the curves are the largest $\mathcal{P}$-vector angles during the training procedure of the first epoch from the random scratch. The results indicate that the principal subspace of the well-trained models are more close to the principal subspace of the raw data, no matter which architectures are used for what kind of learning tasks or whether labels are used for training.

[VggNet16 (Train)]
[VggNet16 (Test)]

[ResNet20 (Train)]
[ResNet20 (Test)]

[ResNet56 (Train)]
[ResNet56 (Test)]
Strong and consistent correlations between the model performance (training and testing accuracy) and the angels between principal subspaces of raw data and deep features  using $\mathcal{P}$-vectors and CIFAR-10 datasets.

[ResNet20 (Train)]
[ResNet20 (Test)]

[ResNet56 (Train)]
[ResNet56 (Test)]

[ResNet110 (Train)]
[ResNet110 (Test)]

Strong and consistent correlations between the model performance (training and testing accuracy) and the angels between principal subspaces of raw data and deep features  using $\mathcal{P}$-vectors and CIFAR-100 datasets.

§.§ Correlation between the Angles and Deep Learning Performance

We further explore the relationship between the performance of the various models on different datasets and the angles between model and data $\mathcal{P}$-vectors for trained ones. Experiments are delivered using ResNet-20/56/110 and VGG-16 on CIFAR-10/100 datasets. As shown in Fig.<ref>, there exists a strong and consistent correlation between the training/testing accuracy of models and the angles between the model and data $\mathcal{P}$-vectors. Note that in these experiments, we use the samples in the training dataset to estimate the model and data $\mathcal{P}$-vectors while avoiding the use of validation information, so as the demonstrate the connection of $\mathcal{P}$-vectors all based on training samples to the generalization performance. For CIFAR-100 dataset, we use ResNets with pre-activation enabled.

To avoid the dominance of some outliers, we conduct correlation analysis between the rank of model performance and the rank of angles by the use of Spearman's coefficients and p-values. With $0.05$ as the threshold for significance, we find significance in the correlations between the angles between model and data $\mathcal{P}$-vectors and the training/testing accuracy for all above cases, except the correlation between the angles and training accuracy on ResNet-100 using CIFAR-100 dataset (p-value=$0.0927$).

In Appendix (A.8), we also include an additional correlation analysis based on log-log plots, where we can further validate our observation. Further for all models on all datasets, the correlation between angles and the testing accuracy is slightly higher than the correlations between angles and training accuracy (and significantly lower p-value).

In this way, we could conclude that (1) both training and testing accuracy are correlated to the angles between the model and data $\mathcal{P}$-vectors, (2) the strong correlations between angles and the testing accuracy might not be caused by the correlations between the angles and training accuracy, as the earlier ones are even stronger, (3) the angles between the model and data $\mathcal{P}$-vectors would be a reasonable performance indicator, as they are strongly, consistently, and significantly correlated to the testing accuracy. This observation could be interpreted that a significant (locally) linear term exists in the well-trained model [13], which makes DNN feature principal subspace correlate to the data principal subspace.

§.§ Applications to Generalization Performance Prediction for Deep Models

To further demonstrate the feasibility of using the $\mathcal{P}$-vector as a “validation-free” measure of generalization performance, we use the experiment settings of “Predicting Generalization in Deep Learning Competition” at NeurIPS 2020 [26] to evaluate the angle between principal subspaces of deep DNN features and raw data on predicting the generalization performance of models using the training dataset. The competition offers a large number of deep models trained with various hyper-parameters and DNN architectures, while the official evaluator for the competition first predicts the generalization performance of every model using the proposed measure, then verify the prediction results through the mutual information (the higher the better) between the proposed measures and the (observable) ground truth of generalization gaps.

In the experiments, we propose to use the angles between the model and data $\mathcal{P}$-vectors using the training dataset as the metrics of generalization performance. In the comparisons with the proposed $\mathcal{P}$-vectors, we include several baseline measures in generalization performance predictors, including VC Dimension [27], Jacobian norm w.r.t intermediate layers [28], Distance from the convergence point to initialization [29], and the Sharpness of convergence point [30]. In addition to these methods, we also propose “Pseudo Validation Accuracy” as a measure for comparisons, where this measure first uses random data augmentation apply to the original set of training data to generate a set of “pseudo validation samples”, then tests the accuracy of the model using “pseudo validation samples”.

Table <ref> presents the comparisons between the proposed measures and baselines. It shows that when the proposed measure – angles between the model and data $\mathcal{P}$-vectors – stands alone, the measures significantly outperform the baseline methods including Jacobian norm w.r.t intermediate layers and the VC dimensions. However, through complementing with other metrics, the metrics based on $\mathcal{P}$-vector angles could be further improved in predicting the generalization performance and finally outperform all baseline methods when combining with “Pseudo validation accuracy”. Note that we combine the results of two metrics through weighted aggregation ([31], with a constant weight 0.05) of two ranking lists that are sorted according to the two metrics respectively.

Methods	   	Prediction Score

VC Dimension [27]	   	0.020

Jacobian norm w.r.t intermediate layers' features [28]	   	2.061

$\mathcal{P}$-vector	   	3.325

Distance to the Initialization Point [29]	   	4.921

Distance to Initialization Point + $\mathcal{P}$-vector	   	4.971

Sharpness of the Convergence Point [30]	   	10.667

Pseudo Validation Accuracy	   	13.531

Pseudo Validation Accuracy + $\mathcal{P}$-vector	   	15.618

Scores of different method to predict the generalization gap using CIFAR-10 and SHVN.

§ CONCLUSION

In this work, we study principal subspaces of deep features learned in DNNs, we first propose the $\mathcal{P}$-vector (the top left singular vector of the #samples$\times$#feature matrix extracted from DNN) as the proxy measurement of the principal subspace and use the angles between $\mathcal{P}$-vectors to compare principal subspaces.

We observe that, no matter which architectures or whether the labels have been used to train the models, the angle between principal subspaces of different DNNs (measured using $\mathcal{P}$-vectors)

would decrease to smaller ones (e.g., around $10^\circ$ for most models in this study with cosine similarity higher than 0.985), when the models are trained using the same dataset. We thus conclude that these DNNs would share a common principal subspace in deep feature space.

Furthermore, during the training procedure from the random scratch, the principal subspace of deep features would

slowly approach to the principal subspace of data vectors from an almost-orthogonal status (e.g., 80$^\circ$–90$^\circ$) to smaller angles (e.g., 50$^\circ$–60$^\circ$).

Finally, we find that angles between principal subspaces of deep features and data vectors are strongly correlated to the performance of models while they are capable of predicting generalization performance, even when $\mathcal{P}$-vectors are all estimated using training dataset only. As was discussed, we believe the empirical observations obtained here are partially due to the local linearity of DNN models [13], while we have validate the significance of $\mathcal{P}$-vectors to measure the principal subspace of deep features in Section  3.2. Our future work may focus on the theoretical understanding to these phenomena.

§ DECLARATIONS

Not applicable.

[1]
Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.
Deep learning.
Nature, 521(7553):436–444, 2015.

[2]
Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley,
Sherjil Ozair, Aaron Courville, and Yoshua Bengio.
Generative adversarial nets.
In Advances in Neural Information Processing Systems, pages
2672–2680, 2014.

[3]
Alec Radford, Luke Metz, and Soumith Chintala.
Unsupervised representation learning with deep convolutional
generative adversarial networks.
arXiv preprint arXiv:1511.06434, 2015.

[4]
Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman.
Deep inside convolutional networks: Visualising image classification
models and saliency maps.
In ICLR Workshop, 2013.

[5]
Tom White.
Sampling generative networks.
arXiv preprint arXiv:1609.04468, 2016.

[6]
Jun-Yan Zhu, Philipp Krähenbühl, Eli Shechtman, and Alexei A Efros.
Generative visual manipulation on the natural image manifold.
In European Conference on Computer Vision, pages 597–613,

[7]
David Bau, Bolei Zhou, Aditya Khosla, Aude Oliva, and Antonio Torralba.
Network dissection: Quantifying interpretability of deep visual
In Conference on Computer Vision and Pattern Recognition, 2017.

[8]
David Bau, Jun-Yan Zhu, Hendrik Strobelt, Bolei Zhou, Joshua B. Tenenbaum,
William T. Freeman, and Antonio Torralba.
Gan dissection: Visualizing and understanding generative adversarial
In International Conference on Learning Representations, 2019.

[9]
Eitan Richardson and Yair Weiss.
On GANs and GMMs.
In Advances in Neural Information Processing Systems, pages
5847–5858, 2018.

[10]
David Berthelot, Colin Raffel, Aurko Roy, and Ian Goodfellow.
Understanding and improving interpolation in autoencoders via an
adversarial regularizer.
In International Conference on Learning Representations, 2019.

[11]
Thilo Spinner, Jonas Körner, Jochen Görtler, and Oliver Deussen.
Towards an interpretable latent space: an intuitive comparison of
autoencoders with variational autoencoders.
In IEEE VIS, 2018.

[12]
Ali Jahanian, Lucy Chai, and Phillip Isola.
On the "steerability" of generative adversarial networks.
In International Conference on Learning Representations, 2020.

[13]
Xiao Zhang and Dongrui Wu.
Empirical studies on the properties of linear regions in deep neural
In International Conference on Learning Representations, 2020.

[14]
Anh Nguyen, Alexey Dosovitskiy, Jason Yosinski, Thomas Brox, and Jeff Clune.
Synthesizing the preferred inputs for neurons in neural networks via
deep generator networks.
In Advances in Neural Information Processing Systems, pages
3387–3395, 2016.

[15]
Georgios Arvanitidis, Lars Kai Hansen, and Søren Hauberg.
Latent space oddity: On the curvature of deep generative models.
In International Conference on Learning Representations, 2018.

[16]
Tom Sercu, Sebastian Gehrmann, Hendrik Strobelt, Payel Das, Inkit Padhi, Cicero
Dos Santos, Kahini Wadhawan, and Vijil Chenthamarakshan.
Interactive visual exploration of latent space (IVELS) for peptide
auto-encoder model selection.
In ICLR DeepGenStruct Workshop, 2019.

[17]
Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.
A simple framework for contrastive learning of visual
arXiv preprint arXiv:2002.05709, 2020.

[18]
Prannay Khosla, Piotr Teterwak, Chen Wang, Aaron Sarna, Yonglong Tian, Phillip
Isola, Aaron Maschinot, Ce Liu, and Dilip Krishnan.
Supervised contrastive learning.
arXiv preprint arXiv:2004.11362, 2020.

[19]
Hervé Abdi and Lynne J Williams.
Principal component analysis.
Wiley interdisciplinary reviews: computational statistics,
2(4):433–459, 2010.

[20]
Andrew M Saxe, James L McClelland, and Surya Ganguli.
A mathematical theory of semantic development in deep neural
Proceedings of the National Academy of Sciences,
116(23):11537–11546, 2019.

[21]
Dennis Lee, Christian Szegedy, Markus N Rabe, Sarah M Loos, and Kshitij Bansal.
Mathematical reasoning in latent space.
arXiv preprint arXiv:1909.11851, 2019.

[22]
Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp.
Finding structure with randomness: Probabilistic algorithms for
constructing approximate matrix decompositions.
SIAM Review, 53(2):217–288, 2011.

[23]
Vladimir Pestov.
On the geometry of similarity search: Dimensionality curse and
concentration of measure.
arXiv preprint cs/9901004, 1999.

[24]
Áke Björck and Gene H Golub.
Numerical methods for computing angles between linear subspaces.
Mathematics of computation, 27(123):579–594, 1973.

[25]
Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.
Imagenet classification with deep convolutional neural networks.
In Advances in Neural Information Processing Systems, pages
1097–1105, 2012.

[26]
Yiding Jiang, Pierre Foret, Scott Yak, Daniel M Roy, Hossein Mobahi,
Gintare Karolina Dziugaite, Samy Bengio, Suriya Gunasekar, Isabelle Guyon,
and Behnam Neyshabur.
Neurips 2020 competition: Predicting generalization in deep learning.
arXiv preprint arXiv:2012.07976, 2020.

[27]
Vladimir Vapnik.
The Nature of Statistical Learning Theory.
Springer Science & Business Media, 2013.

[28]
Yiding Jiang, Pierre Foret, Scott Yak, Behnam Neyshabur, Isabelle Guyon,
Hossein Mobahi, Gintare Karolina, Daniel Roy, Suriya Gunasekar, and Samy
Predicting generalization in deep learning, competition at NeurIPS

[29]
Vaishnavh Nagarajan and J Zico Kolter.
Generalization in deep networks: The role of distance from
arXiv preprint arXiv:1901.01672, 2019.

[30]
Yiding Jiang, Behnam Neyshabur, Hossein Mobahi, Dilip Krishnan, and Samy
Fantastic generalization measures and where to find them.
arXiv preprint arXiv:1912.02178, 2019.

[31]
Vasyl Pihur, Susmita Datta, and Somnath Datta.
RankAggreg, an R package for weighted rank aggregation.
BMC Bioinformatics, 10(1):62, 2009.

§ APPENDIX

This appendix includes results of additional experiments.

§.§ A.0 Comparison of angles between  $\mathcal{P}$-vector extracted from well-trained models Using Different Architectures on ImageNet

[ImageNet (Train)]

[ImageNet (Test)]
Cosine of angles between principal subspaces of deep features, measured using $\mathcal{P}$-vectors, for models trained under default settings.

§.§ Comparison of angles between checkpoints and well-trained $\mathcal{P}$-vector Using Different Architectures on CIFAR-100

In the main text, we presented the result on CIFAR-10 dataset. To generalize the observations, we repeated the experiments on CIFAR-100 dataset to validate our hypothesis of the convergence of the angles between model checkpoints and well-trained model $\mathcal{P}$-vectors. We investigate the change of angles over the $\mathcal{P}$-vectors of training model checkpoints per epoch with comparison to the $\mathcal{P}$-vectors of well-trained models (model of epoch 200 in our case).  As shown in Fig.<ref>, a gradually decreasing manner of the curves for the angle between $\mathcal{P}$-vectors and all angles between $\mathcal{P}$-vectors cross models with different supervisory manners generally converge to a value that smaller than 10$^\circ$ degree. We can conclude that the hypothesis of the existence of common subspace during the learning procedure also stands on the experiments with CIFAR-100 dataset.

[Supervised (Train)]
[Supervised (Test)]

[Unsupervised (Train)]
[Unsupervised (Test)]

[Self-Sup. (Test)]

Convergence to the Common Feature Subspace with CIFAR-100.  Curves of angles of $\mathcal{P}$-vectors between the well-trained model and its checkpoint per training epoch of three learning supervisory manners. The trends of convergence for the angles can be observed in all models.

§.§ The Non-monotonic Trend within the First Epoch

The variation of angles between $\mathcal{P}$-vectors for the well-trained model and its checkpoint per training epoch of each iteration in the first epoch. The non-monotonic trends within the first epoch also incorporate in the experiments using the testing sets of CIFAR-10 and CIFAR-100 datasets. Fig.<ref> shows the curves indicating the variation of angles between the training model $\mathcal{P}$-vectors and the well-trained model $\mathcal{P}$-vectors in the iterations in the first epoch. As we use 128 as the batch size in training procedure, the number of iterations for updates is 391 per epoch. We obtain the observation of a non-monotonic trend that the angle first rises with the random initialization and drop down. And in the rest of training process, the angles keeps the approximately monotonically decreasing and converging to small values. The experiments shows consistent result and conclusion on the testing set of CIFAR-10 and CIFAR-100 with the discussion in section 3.

[Supervised CIFAR-10]
[Supervised CIFAR-100]

[Unsupervised CIFAR-10]
[Unsupervised CIFAR-100]

[Self-Supervised CIFAR-10]
[Self-Supervised CIFAR-100]
Angles between principal subspaces, measured using $\mathcal{P}$-vectors based on the Testing Sets of CIFAR-10 and CIFAR-100, between well-trained models and checkpoints per training iteration in the first epoch.

§.§ Model-to-model Common Subspace

We also test and verify the model-to-model common subspace shared by models trained with different supervisory manners on CIFAR-100 dataset. Experiments carried out to evaluate the angles between $\mathcal{P}$-vectors for checkpoints of all models and $\mathcal{P}$-vectors for well-trained supervised, unsupervised and self-supervised models, where we use the well-trained Wide-ResNet28/Convolution Auto-encoder and SimCLR model (trained with 200 epochs under suggest settings) as the reference of supervised, unsupervised and self-supervised models, respectively. As shown in Fig.<ref>, a consistent convergence for the curves of the angles can be observed and support our hypothesis that the dynamics learning procedure construct the common subspace gradually.

[Sup. vs. Sup.]
[Sup. vs. Unsup.]
[Sup. vs. Self-Sup.]

[Unsup. vs. Sup.]
[Unsup. vs. Unsup.]
[Unsup vs Self-Sup]

[Self. vs Sup.]
[Self. vs Unsup.]
[Self. vs Self-Sup.]

Convergence of $\mathcal{P}$-vector angles between checkpoint per epoch and well-trained models using CIFAR-100.

§.§ The Non-monotonic Trend in the First Epoch of Comparison of Angles between Model and Raw Data $\mathcal{P}$-vectors

We also explore the construction procedure for the common subspace share between feature vectors and the raw data during the training process.  Experiments are carried out to compare the space of models and raw data $\mathcal{P}$-vectors on the training dataset. As shown in Fig.<ref>, we observe a non-monotonic trend that the angle first rises with the random initialization and drop down. The angles keeps the approximately monotonically decreasing and converging to small values in following training epochs. The experiments shows consistent result on both CIFAR-10 and 100 dataset. Note that we follow the default random data augmentation policy to pre-process the training dataset.

[Supervised CIFAR-10]
[Supervised CIFAR-100]

[Unsupervised CIFAR-10]
[Unsupervised CIFAR-100]

[Self-Sup CIFAR-10.]
[Self-Sup CIFAR-100.]

Angles between the $\mathcal{P}$-vectors of the training model and the raw datasets over the number of iterations in the first epoch using CIFAR-10 /CIFAR-100.

§.§ Case Studies on the Angle Variety Between Model and Raw Data for Each Layer

To explore the dynamic variation process by zooming in to every layer variation in each epoch, we perform the case study using Resnet-18 structure on CIFAR-10 dataset. As shown in Fig.<ref>, we give the angles according to layers of 5 different epochs, where the x-axis indicating the indices of residual blocks in the network structure and y-axis refers to the angles between the model checkpoint and raw data $\mathcal{P}$-vectors. We observed that in early training stage, the angles between the model checkpoint and raw data $\mathcal{P}$-vectors keeps an increasing manner when the features passing through layers and turn into a decrease trend towards the stacked layers in the late training stage. This set of experiments further support our hypothesis that the dynamics learning procedure construct the common subspace gradually through training process.

[Epoch 5.]
Angles changes through layer between the $\mathcal{P}$-vectors of the training and raw data over the number of iterations in the first epoch using CIFAR-10.

§.§ Distribution of Values in the $\mathcal{P}$-vector

Please refer to Figure <ref> for the results of experiments carried out on CIFAR-10 and CIFAR-100 datasets using ResNet-50. To have a better view of the distribution drift, we further provide the KDE-smoothed frequency map of the $\mathcal{P}$-vector for wide-resnet training using CIFAR-100 dataset, as shown in Figure <ref>.

[CIFAR 10 Epoch 0]
[CIFAR 10 Epoch 60]
[CIFAR 10 Epoch 120]

[CIFAR 10 Epoch 160]
[CIFAR 10 Epoch 200]

[CIFAR 100 Epoch 0]
[CIFAR 100 Epoch 60]
[CIFAR 100 Epoch 120]

[CIFAR 100 Epoch 160]
[CIFAR 100 Epoch 200]

Frequency of Values appeared in the P-vector versus training epochs.

[CIFAR 100 Epoch 0,1]

[CIFAR 100 Epoch 1,60,120,160,200]

KDE Smoothed Frequency of Values appeared in the P-vector versus training epochs.

§.§ No Convergence found in comparisons between the Top Singular Vectors Other Than the $\mathcal{P}$-vectors

Please refer to Figure <ref> for the results of experiments carried out on CIFAR-10 datasets using ResNet-50.

[Top 1 singular vector.]
[Top 2 singular vector.]

[Top 3 singular vector.]

[Top 4 singular vector.]

[Top 5 singular vector.]

[Top 6 singular vector.]

Angles between the top-k left singular vector ($k=1$ is the $\mathcal{P}$-vector) of the training and well-trained models over the number of epochs in the training process (Resnet-50, CIFAR-10). Note the first plot point refers to the feature matrix after trained for one epoch.

§.§ Log-Log Plots that correlate the $\mathcal{P}$-vector Angles and the Model Performance

Please refer to Figure <ref> for the log-log plot of the results based on CIFAR-10 datasets.

[VggNet16 (Train)]
[VggNet16 (Test)]

[ResNet20 (Train)]
[ResNet20 (Test)]

[ResNet56 (Train)]
[ResNet56 (Test)]

Log-Log Plots: Correlations between the model performance (training and testing accuracy in log range) and the angels (in log range) between model and data $\mathcal{P}$-vectors using CIFAR-10 datasets.

[ResNet20 (Train)]
[ResNet20 (Test)]

[ResNet56 (Train)]
[ResNet56 (Test)]

[ResNet110 (Train)]
[ResNet110 (Test)]

Log-Log Plots: Correlations between the model performance (training and testing accuracy in log range) and the angels (in log range) between model and data $\mathcal{P}$-vectors using CIFAR-100 datasets.

§.§ Explained variance analysis of Top-k singular vector

Explained Variances and the approximation error of the top-$k$ dimensional subspace $E=\|X-U_k \Sigma_k V_k^T\|_F^2$ of top-$k$ singular vectors are shown in Table. <ref>, where $X$ would be the feature matrices, and $U_k$, $\Sigma_k$, $V_k$ are the first  columns in the result of SVD, here we analysis only the changes of top $k=1$ singular vectors through different training checkpoints. The changing procedure for the reconstruction error of approximation is further shown in Figure.<ref>. Also, the variance ratio analysis with various $k$ values is shown in Table.<ref>.

Epoch           $E$      $|X_k|_F^2$     Ratio of $E$

0      246.33         485.18           56.83

20      699.95        1115.90           62.23

40      674.80        1106.20           61.35

60      689.86        1115.35           62.46

80      851.61        1364.77           62.61

100      880.21        1434.66           62.01

120      867.66        1418.27           61.46

140      964.99        1507.03           64.41

160      968.36        1497.79           64.96

180      980.35        1494.50           65.89

200      980.86        1490.68           66.12

Explained variance and reconstruction analysis of top-1 singular vector through the training process.

The reconstruction error of approximation with top-1 singular vector through the training process.

k     Var. ratio(%)     k     Var. ratio(%)

1     56.74      7     88.91

2     64.89      8     92.52

3     70.67      9     95.77

4     75.87      10     98.87

5     80.63      11     98.96

6     84.89      12     99.03

Explained Variances of top-k singular vectors.


# Distributed Ranging SLAM for Multiple Robots with Ultra-WideBand and
Odometry Measurements

Ran Liu, Zhongyuan Deng, Zhiqiang Cao, Muhammad Shalihan, Billy Pik Lik Lau,
Kaixiang Chen,
Kaushik Bhowmik, Chau Yuen, and U-Xuan Tan This work is supported by National
Key R&D Program of China (2019YFB1310805) and National Science Foundation of
China (12175187).R. Liu, M. Shalihan, B. P. L. Lau, K. Bhowmik, C. Yuen, and
U-X. Tan are with the Singapore University of Technology and Design, Singapore
487372. {ran_liu, billy_lau, kaushik_bhowmik, yuenchau,
<EMAIL_ADDRESS>Deng, Z. Cao, and K. Chen are with the Southwest
University of Science and Technology, Mianyang, Sichuan, China 621010.

###### Abstract

To accomplish task efficiently in a multiple robots system, a problem that has
to be addressed is Simultaneous Localization and Mapping (SLAM). LiDAR (Light
Detection and Ranging) has been used for many SLAM solutions due to its superb
accuracy, but its performance degrades in featureless environments, like
tunnels or long corridors. Centralized SLAM solves the problem with a cloud
server, which requires a huge amount of computational resources and lacks
robustness against central node failure. To address these issues, we present a
distributed SLAM solution to estimate the trajectory of a group of robots
using Ultra-WideBand (UWB) ranging and odometry measurements. The proposed
approach distributes the processing among the robot team and significantly
mitigates the computation concern emerged from the centralized SLAM. Our
solution determines the relative pose (also known as loop closure) between two
robots by minimizing the UWB ranging measurements taken at different positions
when the robots are in close proximity. UWB provides a good distance measure
in line-of-sight conditions, but retrieving a precise pose estimation remains
a challenge, due to ranging noise and unpredictable path traveled by the
robot. To deal with the suspicious loop closures, we use Pairwise Consistency
Maximization (PCM) to examine the quality of loop closures and perform outlier
rejections. The filtered loop closures are then fused with odometry in a
distributed pose graph optimization (DPGO) module to recover the full
trajectory of the robot team. Extensive experiments are conducted to validate
the effectiveness of the proposed approach.

##  I Introduction

Multiple robots system collectively accomplishes a task by the collaboration
of a group of robots, which is shown to be efficient and robust when compared
to single robot solution [1]. Localization is crucial for the robot team to
coordinate with each other. The literature shows a number of mature techniques
and implementations for localization in known infrastructures. However, a
prior knowledge of the environment is not always available in some scenarios.
For example, emergency response often requires a group of robots to explore
unknown site within a minimal time.

Therefore, it is crucial for a group of robots to realize Simultaneous
Localization and Mapping (SLAM) in unknown environments. Most multiple robots
SLAM systems are based on centralized or distributed architecture [2] [3].
Centralized solution processes all collected measurements in a server [2]. The
computational cost increases exponentially with the increase of the number of
robots. Distributed system outperforms centralized solution in scalability and
computational efficiency, as the computation is allocated among a group of
robots [3]. We are interested in distributed SLAM solution for multiple robots
during the exploration of unknown environments.

Figure 1: Overview of the proposed distributed ranging SLAM with two robots.
Each robot carries a UWB node for ranging of the nearby robots. We perform
distributed pose estimation with short-term UWB ranging and odometry
measurements (see the two cyan rectangles) when the two robots are in close
proximity. The estimated poses are passed for outlier rejection and then
optimized through a distributed pose graph optimization module.

SLAM has been extensively researched by the robotics community with LiDAR and
visual-based sensors [4], which exhibit good localization accuracy in general,
but it is challenging to deal with featureless scenes (for example tunnels)
and the algorithms often require a lot of computational resources. Due to its
high accuracy and detection without line-of-sight, Ultra-WideBand (UWB) is
widely used for localization [5] [6]. Although a precise distance is provided
by UWB, the lack of bearing information makes it difficult to be applied for
the navigation of autonomous robots. The conventional Angle-of-Arrival
(AoA)-based solution requires a cumbersome setup for phase shift estimation
[7]. We focus on the use of purely UWB ranging and odometry for pose
estimation among a robot team. The proposed solution is supposed to work in
non-line-of-sight conditions, where LiDAR and visual sensors often fail to
perform localization due to the block of the environment.

This paper presents a distributed SLAM for trajectory estimation of a group of
robots based on UWB ranging and odometry measurements. Each robot is equipped
with a UWB node for ranging, as shown in Figure 1. We estimate the pose
(position and orientation) between the robots using short-term UWB ranging and
odometry observations. In particular, we achieve the distributed pose
estimation (i.e., loop closure) through the minimization of ranging
measurements taken at different positions when two robots are close. The UWB
is susceptible to non-line-of-sight (NLOS) propagation [8], which results in
erroneous loop closures, in particular when the path traveled by the robot is
unpredictable [9]. To remove erroneous loop closures, we utilize the
distributed Pairwise Consistency Maximization (PCM) algorithm to verify the
quality of loop closure pairs and perform oulier rejections. The inlier loop
closures are finally incorporated with the odometry through Distributed Pose
Graph Optimization (DPGO) to estimate the trajectory of the robots. We
evaluated the performance of the proposed approach with three robots in an
indoor environment. The contributions of this paper are summarized as follows:

  1. 1.

We propose an approach for distributed pose estimation between two robots
using short-term UWB ranging and odometry measurements. In particular, we
combine coarse search and optimization strategy to avoid the convergence to a
local minimum.

  2. 2.

We utilize distributed outlier rejection to remove the erroneous estimations
and apply distributed pose graph optimization for trajectory estimation. To
the best of our knowledge, this is the first work that uses UWB ranging for
distributed SLAM with multiple robots.

  3. 3.

We perform extensive experiments to evaluate the performance of the proposed
approach. The results show that we achieved a positioning accuracy of 0.45m in
translation and 4.01$\degree$ in rotation by the fusing of UWB and odometry
measurements in an indoor environment with a size of 10m$\times$12m.

We organize the remaining of this paper as follows: Section II introduces the
related work. Section III describes the proposed distributed ranging SLAM.
Section IV shows the experimental setups and the experimental results.
Finally, Section V concludes this paper and discusses the future work.

##  II Related Work

The robotics research community is showing a growing interest in multiple
robots system, due to the efficiency by collaboration of a robot team. The
system is extremely useful when multiple robots with different functionalities
are used for exploration of unknown environments, like underground exploration
in DARPA subterranean challenge111https://www.subtchallenge.com/. SLAM is
considered fundamental to achieve this task due to its capability of
localization and mapping in unknown environments. Tremendous progress has been
made concerning multiple robots SLAM with different sensors and
implementations.

Visual cameras and LiDARs are popular for localization and SLAM in robotics,
due to their high accuracy [4] [10] [11]. Kimera-multi [4] is presented for
robust distributed trajectory estimation and 3D metric-semantic meshing in
multiple robots applications. Ebadi et al. [10] implemented a LiDAR-based SLAM
and tested it for the exploration of underground tunnels with multiple robots.
CamVox [11] is a low-cost and robust SLAM solution that combines Livox LiDAR
and visual features for robust sensor fusion.

However, the performance of LiDAR and visual SLAM decreases in perceptually-
degraded environments, like tunnels and long corridors. Due to detection
without line-of-sight and high accuracy, UWB has been used for localization in
many industrial applications. A sensor fusion scheme is presented in [12] to
combine visual, IMU, and UWB for accurate and low-drift localization. Liu et
al. [13] combined short-range LiDAR and UWB for mapping of mobile robots in
indoor environments. Cao et al. [14] proposed a solution to localize a moving
robot with a single UWB anchor based on the fusion of IMU with an Extended
Kalman Filter. Funabiki et al. [15] presented an approach to combine UWB
beacon and LiDAR-based geometric loop closures for SLAM in large and
perceptually challenging environments.

The implementation of multiple robots SLAM can be briefly classified into
centralized and distributed. The centralized SLAM requires a server for data
collection and sensor fusion. Cao et al. [16] used four UWB nodes with a
square configuration to determine the pose between two robots. But equipping
multiple UWB sensors is costly and is not suitable for small robot
applications. Liu et al. [1] proposed an approach to combine IMU and UWB for
localization of multiple users. Martel [17] et al. presented an approach for
relative pose estimation with six ranging measurements based on linear
algorithm. Nguyen et al. [18] presented an approach to fuse UWB, IMU, and
visual measurements to achieve multi-robot collaborative localization. A
centralized system is not scalable with large robot teams, due to the increase
of computation and communication requirements.

In distributed SLAM, the team works collaboratively to achieve the goal of
SLAM and the computation is distributed to peer robots. The Gauss-Seidel
approach is presented in [19] for distributed trajectory estimation, which is
shown to scale well to large robot teams. Tian et al. [3] presented a
distributed SLAM solution based on sparse semidefinite relaxation, which is
guaranteed to be globally optimal under moderate measurement noise. This
approach is improved by [20] to balance the computational and communication
cost. Lajoie et al. [21] presented DOOR-SLAM, which is a fully distributed
SLAM system based on robust place recognition.

Figure 2: Overview of the modules and their communication requirements for our
distributed ranging SLAM with two robots.

##  III Distributed Ranging SLAM for Multiple Robots

In this section, we formulate the problem of distributed SLAM with UWB ranging
and odometry for a group of robots. Then, we show the details of our proposed
distributed SLAM, which mainly consists of the following three modules (see
Figure 2). Firstly, each robot performs distributed pose estimation to
determine the pose of its peer robots using short-term UWB ranging and
odometry measurements. Secondly, each robot runs a distributed outlier
rejection module to filter out the suspicious estimations. Thirdly, a
distributed SLAM framework is used to fuse the inliers and the odometry to
produce the trajectory in a common reference frame.

The robot team consists of $N$ robots and we denote the pose of robot $\alpha$
($1\leq\alpha\leq N$) at time $t$ as
${{\mathbf{x}}_{\alpha}^{t}}=[\mathbf{R}_{\alpha}^{t},\mathbf{t}_{\alpha}^{t}]\in\text{SE}(2)$,
where $\mathbf{R}_{\alpha}^{t}\in\text{SO}(2)$ represents the rotation (i.e.,
$\theta_{\alpha}^{t}$) and $\mathbf{t}_{\alpha}^{t}\in\mathbb{R}^{2}$
represents the 2D position (i.e., $x_{\alpha}^{t}$ and $y_{\alpha}^{t}$).
${{\mathbf{x}}_{\alpha}^{k,i}}$ denotes the relative pose from time $k$ to
time $i$ measured by odometry of robot $\alpha$. We use
$\mathbf{r}_{\alpha}^{t}=\\{r_{\alpha,\beta}^{t}\\}$ ($\beta\in
N_{\alpha}^{t})$ to denote the UWB ranging measurements received by robot
$\alpha$ at time $t$, where $N_{\alpha}^{t}$ denotes the set of neighboring
robots of robot $\alpha$ and $r_{\alpha,\beta}^{t}$ represents the UWB ranging
between robot $\alpha$ and robot $\beta$ at time $t$. An example of the
problem can be illustrated in Figure 2, where two robots are performing
distributed SLAM in unknown environment. Each robot carries a UWB node for
ranging to its peer robots in close proximity. The robot also offers odometry
to represent its relative movement. The goal is to achieve the trajectory
estimation for multiple mobile robots through the UWB ranging and odometry
measurements without a prior knowledge about the infrastructure. Our approach
works in a distributed fashion, where each robot will perform its pose
estimation by considering the measurements (i.e., odometry and UWB) received
from all neighboring robots.

###  III-A Module 1: Distributed Pose Estimation with Short-term UWB and
Odometry Measurements

The estimation of the relative pose $\overline{\mathbf{x}}_{\alpha,\beta}^{t}$
(i.e., pose of robot $\beta$ with reference to robot $\alpha$ at time $t$) can
be achieved by finding the best pose configuration
$\overline{\mathbf{x}}_{\alpha,\beta}^{t}$ through minimizing the residual
error of UWB ranging taken at different positions (denoted by the odometry)
within a sliding time window $\tau$:

$\small\begin{split}\operatorname*{arg\,min}_{\mathbf{x}_{\alpha,\beta}^{t}}\sum_{i=t-\tau}^{i=t}\underbrace{\mathbf{e}(\mathbf{x}_{\alpha,\beta}^{t},\mathbf{x}_{\alpha}^{t},\mathbf{x}_{\beta}^{t},\mathbf{x}_{\alpha}^{i},\mathbf{x}_{\beta}^{i},r_{\alpha,\beta}^{i})}_{\text{Residual
error of UWB
ranging}}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,&\\\
=\operatorname*{arg\,min}_{\mathbf{x}_{\alpha,\beta}^{t}}\sum_{i=t-\tau}^{i=t}(\underbrace{r_{\alpha,\beta}^{i}}_{\text{Ranging}}{-}||{\underbrace{{(\mathbf{x}_{\alpha}^{t})}^{{\text{-}1}}\cdot\mathbf{x}_{\alpha}^{i}}_{\text{Odom.
of
$\alpha$}},\mathbf{x}_{\alpha,\beta}^{t}\cdot\underbrace{{(\mathbf{x}_{\beta}^{t})}^{{\text{-}1}}\cdot\mathbf{x}_{\beta}^{i}}_{\text{Odom.
of $\beta$}}}||)^{2},\end{split}$ (1)

where $\mathbf{e}(\cdot)$ denotes the residual function that computes the
residual error of a UWB ranging given a pose configuration
$\mathbf{x}_{\alpha,\beta}^{t}$. $||\cdot||$ computes the Euclidean distance
between two positions. $\tau$ denotes the size of the sliding window for pose
estimation. The optimization of Equation 1 can be solved efficiently by state-
of-the-art optimization tools, such as g2o [22], gtsam [23], and Ceres [24],
if a good initial guess of $\mathbf{x}$ is provided. Due to the non-convexity
of Equation 1, the optimization converge to a local minimum without a proper
initial guess and there is no guarantee to find the best solution.

Data: Odometry from robot $\alpha$: $\mathbf{x}_{\alpha}^{i}$, odometry from
robot $\beta$: $\mathbf{x}_{\beta}^{i}$, UWB ranging between $\alpha$ and
$\beta$: $\mathbf{r}_{\alpha,\beta}^{i}$, where $t-\tau\leq i\leq t$

Result: Coarse pose estimation $\widetilde{\mathbf{x}}_{\alpha,\beta}^{t}$

// Record the minimum residual

1 $minimumResidual\leftarrow\infty$

2 for  _$i_{\phi}\leftarrow-w_{\phi}$ to $w_{\phi}$_ do

3 for  _$i_{\theta}\leftarrow-w_{\theta}$ to $w_{\theta}$_ do

4 $\mathbf{x}_{\alpha,\beta}^{\prime}=<r_{\alpha,\beta}^{t}\cos(\delta
i_{\phi}),r_{\alpha,\beta}^{t}\sin(\delta i_{\phi}),\delta i_{\theta}>$

// Initialize the residual to be 0

5 $r\leftarrow 0$

6 for  _$i\leftarrow t-\tau$ to $t$_ do

// Residual with Equation 1

7
$r=r+\mathbf{e}(\mathbf{x}_{\alpha,\beta}^{\prime},\mathbf{x}_{\alpha}^{t},\mathbf{x}_{\beta}^{t},\mathbf{x}_{\alpha}^{i},\mathbf{x}_{\beta}^{i},r_{\alpha,\beta}^{i})$

// Bound residual error for fast computation

8 if  _$r\geq minimumResidual$_ then

9 break

10 end if

11

12 end for

13 if  _$r <minimumResidual$_ then

14 $minimumResidual\leftarrow r$

15
$\widetilde{\mathbf{x}}_{\alpha,\beta}^{t}=<x_{\alpha,\beta}^{\prime},y_{\alpha,\beta}^{\prime},\theta_{\alpha,\beta}^{\prime}>$

16

17 end if

18

19 end for

20

21 end for

return $\widetilde{\mathbf{x}}_{\alpha,\beta}^{t}$

Algorithm 1 Find coarse relative pose
$\widetilde{\mathbf{x}}_{\alpha,\beta}^{t}$ between robot $\alpha$ and robot
$\beta$ with short-term UWB and odom.

(a) Odom. trajectories for Case 1

(b) Odom. trajectories for Case 2

(c) Mean residual for Case 1

(d) Mean residual for Case 2

Figure 3: Trajectories of two test cases with a window size $\tau=50$ and the
mean residual error for different polar angle $\phi$ and relative angle
$\theta$. The global minimal can be found by an arbitrary initial guess for
Case 1, however this is not guaranteed in Case 2.

We propose to use a coarse search strategy to find a decent approximation of
relative pose as initial guess for optimization in Equation 1. Based on coarse
guess, we then perform nonlinear optimization with Equation 1. As UWB ranging
is considered to be accurate, robot $\beta$ will approximately located on a
circle centered at robot $\alpha$ with a radius $r_{\alpha,\beta}^{t}$ (known
by UWB ranging) and an angle $\phi$ (to be estimated) in polar coordinate
system. Therefore, instead of searching the whole 2D space with Equation 1, we
search in a polar space with angle between $-\pi$ and $\pi$ given the known
UWB ranging. The angular searching step size for $\theta$ and $\phi$ is
denoted as $\delta$. The number of searching steps (i.e., $w_{\phi}$ and
$w_{\theta}$) for $\theta$ and $\phi$ are denoted as
$w_{\phi}=w_{\theta}=\lceil\frac{\pi}{\delta}\rceil$. We define the entire
searching space as
$\overline{w}=\\{-w_{\phi},w_{\phi}\\}\times\\{-w_{\theta},w_{\theta}\\}$. The
optimization of Equation 1 is approximated to find the minimum residual given
${\widetilde{\mathbf{x}}_{\alpha,\beta}^{t}}\in\mathbf{w}$:

$\small\begin{split}\mathbf{w}=\\{(\underbrace{r_{\alpha,\beta}^{t}\cos(\delta
i_{\phi})}_{\text{Position $x$}},\underbrace{r_{\alpha,\beta}^{t}\sin(\delta
i_{\phi})}_{\text{Position $y$}},\delta
i_{\theta}):(i_{\phi},i_{\theta})\in\overline{w}\\}\end{split}$ (2)

An overview of the algorithm is shown in Algorithm 1. The residual errors with
respect to different $\phi$ and $\theta$ for two different test cases are
shown in Figure 3. The figure shows the residual follows a multimodal
distribution, and a global minimum might be not found from an arbitrary
initial guess (for example the Case 2 in Figure 3). The coarse solution
${\widetilde{\mathbf{x}}_{\alpha,\beta}^{t}}$ is used as initial guess to
optimize Equation 1 via g2o [22]. The node in the graph is denoted as the pose
to be estimated (i.e., $\overline{\mathbf{x}}_{\alpha,\beta}^{t}$) and the
constraints are represented by the UWB ranging measurements. The algorithm
turns out to find the best configuration of the pose to satisfy the UWB
ranging constraints based on maximum likelihood estimation.

###  III-B Module 2: Distributed Outlier Rejection

Figure 4: Illustration of the distributed PCM consistency check performed by
robot $\alpha$. Robot $\alpha$ estimates the pose of robot $\beta$ at time $k$
and time $i$ using distributed pose estimation described in Section III-A. The
estimated pose and odometry are passed for consistency check.

The pose estimation accuracy in Section III-A highly depends on UWB ranging
accuracy and paths of robots. These factors are not controllable for a
practical applications, especially during autonomous navigation of robot in
unknown environment with non-line-of-sight conditions. Therefore, we use
state-of-the-art outlier rejection algorithm, namely Pairwise Consistency
Maximization (PCM) [25], to deal with suspicious pose estimations (i.e., loop
closures) and tailor it to a distributed and online fashion. Given two loop
closures $\overline{\mathbf{x}}_{\alpha,\beta}^{k}$ and
$\overline{\mathbf{x}}_{\alpha,\beta}^{i}$ between robot $\alpha$ and $\beta$,
the following equation is used to determine if two loop closures are pairwise
consistent:

$\small\begin{split}||\mathbf{x}_{\alpha}^{k,i}\oplus\overline{\mathbf{x}}_{\alpha,\beta}^{i}\oplus\mathbf{x}_{\beta}^{i,k}\ominus\overline{\mathbf{x}}_{\alpha,\beta}^{k}||_{\Sigma}\leq\chi^{2}_{\epsilon},\end{split}$
(3)

where $||\cdot||_{\Sigma}$ represents the Mahalanobis distance. $\oplus$ and
$\ominus$ denote composition and inversion operators of two poses [26]. An
illustration of consistency check is shown in Figure 4. In ideal conditions,
where all pose estimations are without noise, the composition of the poses
along the cycle, which is formed by the two odometries (green lines in Figure
4) and two loop closures (cyan lines in Figure 4) is identify. Equation 3
evaluates if the error accumulated along the cycle is consistent with the
noise covariance $\Sigma$. $\chi^{2}_{\epsilon}$ is the quantile of the Chi-
squared test for a given significance level $\epsilon$ [21].

After consistency check by Equation 3, the PCM builds an adjacency matrix to
record the pairwise consistency of the loop closures. The robot computes the
maximum clique from the matrix, which gives the largest set of pairwise
consistent loop closures (i.e., inliers). To reduce computational cost for
online operations, we use an efficient version of PCM [4], which incrementally
updates the consistent loop closures with a heuristic search strategy to find
the maximum clique. The inliers are passed to the pose graph for optimization,
which will be detailed in Section III-C.

###  III-C Module 3: Distributed Pose Graph Optimization

Early multi-robot SLAM solutions focus on the centralized optimization. The
computation increases with the increase of robots, making this solution
challenging in large environment with a large number of robots. In distributed
SLAM [19], each robot solves its pose graph optimization based on local
information. The robots exchange their estimations to iteratively solve the
global optimization problem. We follow the approach proposed in [3]. Each
robot runs a DPGO module for trajectory optimization based on the odometry and
the inlier loop closures. Particularly, DPGO performs local optimization by
Riemann gradient descent. It solves a rank-restricted relaxation of pose graph
optimization in a distributed manner. During the update, the algorithm only
requires the exchange of the poses connecting the loop closures and therefore
saves a large amount of communication resources by avoiding the exchange of
all pose information.

Communication between robots is crucial in distributed SLAM system. The
exchange of odometry is required in two modules, namely distributed pose
estimation module in Section III-A and distributed outlier rejection in
Section III-B. The former shares short-term odometry measurements when two
robots are in close proximity. The latter uses the relative pose estimation
and its associated odometry for consistency check. Our solution is
particularly suitable for the applications in ad hoc networks where external
communication infrastructure is not deployed. For ad hoc networks, the robots
are moving freely to form dynamic network topology, which leads to a limited
communication range. Our implementation exchanges the updated poses and loop
closures when the robots are in close proximity, while the communication is
not required when the two robots are out of the UWB ranging.

##  IV Experimental Results

###  IV-A Experimental Setups

We present experimental setups to demonstrate the proposed approach with three
robots in an indoor environment with a size of 10m$\times$12m. In our setup,
each robot carried one UWB node (Nooploop
LinkTrack111https://www.nooploop.com/) with a maximum range up to 100 meters.
The sampling rate of UWB is set to 50Hz. The robot outputs the odometry
measurements (differential drive with encoders) with a frequency of 10Hz.
Three robots were manually controlled to move along different paths.

Figure 5: Experiments with three TurtleBots moving in an indoor environment. TABLE I: Relative pose estimation accuracy (translational error in meters and rotational error in degrees) and the time consumption (in millisecond) with various time window $\tau$ (in seconds) and different approaches. $\tau$ |  NLS optimization[27] | Coarse search | Proposed approach (coarse search+opt.)
---|---|---|---
Trans. error | Rot. error | Time | Trans. error | Rot. error | Time | Trans. error | Rot. error | Time
$\tau=$10 |  8.19$\pm$6.85 |  70.71$\pm$53.85 | 2ms |  5.03$\pm$5.19 |  54.56$\pm$50.74 | 4ms |  4.54$\pm$5.24 |  49.34$\pm$50.61 | 6ms
$\tau=$25 |  7.47$\pm$6.89 |  63.48$\pm$56.47 | 4ms |  2.89$\pm$3.38 |  31.82$\pm$41.97 | 5ms |  2.42$\pm$3.37 |  27.08$\pm$41.97 | 9ms
$\tau=$50 |  6.12$\pm$6.43 |  54.11$\pm$56.25 | 6ms |  1.88$\pm$2.21 |  15.58$\pm$23.39 | 7ms |  1.38$\pm$2.01 |  12.66$\pm$24.15 | 13ms
$\tau=$100 |  4.53$\pm$5.74 |  37.16$\pm$48.49 | 10ms |  1.31$\pm$1.45 |  9.23$\pm$11.07 | 9ms |  0.83$\pm$1.17 |  6.74$\pm$11.49 | 19ms
$\tau=$200 |  3.97$\pm$5.49 |  32.56$\pm$46.67 | 15ms |  1.25$\pm$1.49 |  7.90$\pm$9.67 | 14ms |  0.82$\pm$1.24 |  6.29$\pm$11.09 | 29ms

(a) Average translational error under various settings of $\delta$

(b) Average rotational error under various settings of $\delta$

(c) Average time consumption under various settings of $\delta$

Figure 6: Relative pose estimation accuracy and time consumption under
different angular searching step sizes $\delta$.

All modules are running under ROS (Robot Operating System). During the
experiment, the maximum velocity of the robot was 0.2m/s. The total distances
travelled by the three robots were 63.5m, 60.6m, and 57.1m, respectively. To
obtain the ground truth, Hokuyo LiDARs were installed to perform adaptive
Monte Carlo localization222http://wiki.ros.org/amcl given a map created by
GMapping. The number of poses collected by the three robots for evaluation
were 314, 363, and 333, respectively. An overview of the experimental snapshot
is shown in Figure 5. The relative translational and rotational errors are
computed by comparing the estimated poses with the ground truth. The mean
squared error (MSE) in translation and rotation are computed as a metric to
evaluate the positioning accuracy.

(a) Average translational error under various ranging setting and ranging
noise

(b) Average rotational error under various ranging setting and ranging noise

Figure 7: Relative pose estimation accuracy with respect to UWB maximum range
and different UWB ranging noise.

(a) Average translational error under various odometry noise

(b) Average rotational error under various odometry noise

Figure 8: Relative pose estimation accuracy with respect to odometry noise.

(a) Pose estimation of Robot2 and Robot3 with respect to Robot1 using short-
term UWB and odometry measurements

(b) Inliers by performing distributed PCM consistency check

(c) Distributed pose graph optimization after PCM check

Figure 9: Trajectories estimated by different approaches. (a) Pose estimation using short-term UWB and odometry measurements; (b) Filtered poses after performing PCM check; (c) Optimized trajectory using distributed pose graph optimization. TABLE II: Pose estimation accuracy (translational error in meters and rotational error in degrees) with PCM check at different significance levels $\epsilon$ and evaluation of distributed pose graph optimization. $\tau$ | Distributed PCM consistency check | Distributed pose graph optimization
---|---|---
$\epsilon$ in PCM | No. of inliers | Trans. error | Rot. error | Trans. error | Rot. error
$\tau$=10 | Without PCM | 999 |  4.54$\pm$5.24 |  49.34$\pm$50.61 |  2.13$\pm$1.12 |  8.72$\pm$5.27
$\epsilon$=0.01 | 391 |  1.57$\pm$1.31 |  12.15$\pm$9.98 |  0.78$\pm$0.45 |  6.08$\pm$3.60
$\epsilon$=0.05 | 269 |  1.37$\pm$1.07 |  9.59$\pm$8.35 |  0.79$\pm$0.43 |  6.52$\pm$4.07
$\epsilon$=0.1 | 212 |  1.12$\pm$0.82 |  8.09$\pm$6.22 |  0.61$\pm$0.36 |  5.19$\pm$3.18
$\epsilon$=0.2 | 152 |  0.97$\pm$0.82 |  6.47$\pm$4.54 |  0.71$\pm$0.43 |  5.69$\pm$3.56
$\epsilon$=0.5 | 65 |  0.57$\pm$0.42 |  5.68$\pm$3.40 |  0.53$\pm$0.36 |  5.18$\pm$3.38
|  $\epsilon$=0.8 | 25 |  0.73$\pm$0.72 |  7.97$\pm$6.02 |  0.69$\pm$0.38 |  6.60$\pm$3.96
$\tau$=50 | Without PCM | 999 |  1.38$\pm$2.01 |  12.66$\pm$24.15 |  0.74$\pm$0.46 |  5.42$\pm$4.05
$\epsilon$=0.01 | 904 |  0.95$\pm$0.91 |  6.93$\pm$6.75 |  0.63$\pm$0.42 |  4.78$\pm$3.50
$\epsilon$=0.05 | 863 |  0.90$\pm$0.83 |  6.14$\pm$5.16 |  0.60$\pm$0.37 |  4.29$\pm$3.07
$\epsilon$=0.1 | 782 |  0.76$\pm$0.59 |  5.78$\pm$4.30 |  0.52$\pm$0.33 |  4.43$\pm$3.00
$\epsilon$=0.2 | 663 |  0.64$\pm$0.42 |  5.25$\pm$3.71 |  0.49$\pm$0.32 |  4.18$\pm$2.83
$\epsilon$=0.5 | 394 |  0.51$\pm$0.31 |  4.52$\pm$2.99 |  0.45$\pm$0.30 |  4.01$\pm$2.81
|  $\epsilon$=0.8 | 131 |  0.53$\pm$0.36 |  4.35$\pm$3.78 |  0.52$\pm$0.26 |  4.23$\pm$2.75
$\tau$=100 | Without PCM | 999 |  0.83$\pm$1.17 |  6.74$\pm$11.49 |  0.60$\pm$0.41 |  5.05$\pm$3.56
$\epsilon$=0.01 | 987 |  0.75$\pm$0.76 |  5.88$\pm$4.57 |  0.59$\pm$0.40 |  4.96$\pm$3.44
$\epsilon$=0.05 | 974 |  0.74$\pm$0.74 |  5.71$\pm$4.15 |  0.57$\pm$0.37 |  4.88$\pm$3.34
$\epsilon$=0.1 | 929 |  0.69$\pm$0.63 |  5.42$\pm$3.71 |  0.53$\pm$0.34 |  4.71$\pm$3.22
$\epsilon$=0.2 | 736 |  0.64$\pm$0.49 |  5.19$\pm$3.48 |  0.53$\pm$0.34 |  4.69$\pm$3.29
$\epsilon$=0.5 | 702 |  0.51$\pm$0.34 |  5.47$\pm$3.42 |  0.48$\pm$0.32 |  4.92$\pm$3.32
|  $\epsilon$=0.8 | 252 |  0.54$\pm$0.36 |  5.08$\pm$3.10 |  0.49$\pm$0.31 |  4.89$\pm$3.19

###  IV-B Evaluation of Distributed Pose Estimation using Short-term UWB and
Odometry

####  IV-B1 Impact of sliding window size $\tau$

In the first set of experiments, we evaluated the performance of our pose
estimation using UWB ranging and odometry measurements. We fixed the angular
searching step size $\delta$ to 0.1 and examined the positioning accuracy and
average time consumption for one pose estimation under the impact of sliding
window size $\tau$, as shown in Table I. This table also lists a comparison of
the accuracy between nonlinear least square (NLS) optimization [27], coarse
search-based approach, and the combination of coarse search and NLS
optimization. As can be seen from Table I, our approach, which combines the
coarse search and NLS optimization, provides the best accuracy for all
settings of $\tau$. For example, the combined approach gives a localization
accuracy of 0.83m in translation and 6.74$\degree$ in rotation for $\tau$=100,
which lead to improvements of 36.6% and 27.0% when compared with the pure
coarse search-based approach (accuracy of 1.31m in translation and
9.23$\degree$ in rotation). We also observed that the localization accuracy
improves with the increase of the window size $\tau$, but at a cost of more
computational resources.

####  IV-B2 Impact under different angular searching steps $\delta$

We evaluated the performance of relative pose estimation under various angular
searching steps $\delta$. We fixed the sliding window size $\tau$=100. The
positioning accuracy and the average computational time for one pose
estimation are shown in Figure 6. As can be seen from this figure, for all
parameter settings, optimization with a coarse initial guess improves the
localization accuracy when compared to the pure coarse search-based approach.
The localization error is reduced with the decrease of the angular searching
step, while the computational time increases at the same time.

####  IV-B3 Impact of different UWB and odometry settings

In this series of experiments, we set the UWB ranging limit (i.e., 20%, 40%,
60%, and 80% of the maximum UWB range reported during the test) to compare the
accuracy. We also added Gaussian noise to raw UWB ranging to verify the
performance under various noise levels. We fixed the searching step size
$\delta$=0.1 and the sliding window size $\tau$=100. The result is shown in
Figure 7. As can be seen from this figure, the localization accuracy is
reduced with short UWB range. We also observed long UWB range in general helps
to improve the localization accuracy due to the large number of UWB
observations to support pose estimation. However, longer UWB ranging brings
more uncertainty, particularly in non-line-of-sight conditions, which leads to
a decrease of the localization accuracy. For all experimental settings, adding
ranging noise will obviously lead to a decrease of localization accuracy.

We compared the localization accuracy with respect to various odometry errors.
We added Gaussian noises to the position and heading of odometry and show the
results in Figure 8. As can be seen from this figure, the accuracy decreases
with the increase of odometry noise. However, we believe a good odometry can
be provided by the on-board sensors of the robot, in particularly with the
assistance of the state-of-the-art visual odometry or LiDAR odometry
techniques.

TABLE III: Comparison between our approach and the literature (UWB SLAM in [18]). Approach | Pure odom. | Ours |  Literature[18]
---|---|---|---
Tran. error(m) |  0.85$\pm$0.52 |  0.45$\pm$0.30 |  0.71$\pm$0.52
Rot. error(degree) |  8.99$\pm$5.8 |  4.01$\pm$2.81 |  6.76$\pm$4.79

###  IV-C Distributed Pose Graph Optimization with PCM Outlier Rejection

In this series of experiments, we examined the localization accuracy with PCM
check as well as the performance of the distributed pose graph optimization.
We set the angular searching step size $\delta$=0.1. The results are listed in
Table II. As can be seen from this table, PCM filters out the suspicious loop
closures and improves the overall localization accuracy. For a sliding window
size $\tau$=10, we obtain a localization accuracy of 1.37m in translation and
9.59$\degree$ in rotation with significance level $\epsilon=0.05$, which gives
improvements of 69.8% and 80.6% when compared to the case without PCM check
(4.54m in translation and 49.34$\degree$ in rotation). In addition, we observe
that the number of inliers decreases with the increase of the significance
level $\epsilon$ in PCM. A too large $\epsilon$ will harm the online
operations, particularly during the early stage of trajectory estimation:
potential good loop closures are pruned and the pose graph optimization has to
be postponed until a reasonable number of loop closures are reached.
Meanwhile, a small $\epsilon$ is not able to distinguish the erroneous loop
closures and results in an increase of the localization error. The distributed
pose graph optimization further improves the localization accuracy by the
additional incorporation of odometry information. For $\tau$=50 and
$\epsilon$=0.5, we obtained a localization accuracy of 0.45m in translation
and 4.01$\degree$ in rotation (the best accuracy achieved), with improvements
of 11.8% and 11.3% when compared to the results without pose graph
optimization (i.e., 0.51m in translation and 4.52$\degree$ in rotation). The
pose estimated with UWB and odometry, inliers after PCM check, and estimated
trajectory using distributed pose graph optimization are visualized in Figure
9.

Table III compared our approach with the conventional UWB SLAM [18], that uses
UWB ranging for SLAM in a centralized manner. As can be seen from this table,
our approach outperforms the conventional UWB SLAM in localization accuracy.
It is important to point out that the conventional UWB SLAM requires initial
positions of the robots to be known publicly, while our approach discards this
assumption by performing the relative pose estimation using local UWB ranging
and odometry measurements.

###  IV-D Computation and Communication Efficiency

Finally, we evaluated the computational cost and the communication resources
involved in our approach, as listed in Table IV. We choose the window size
$\tau$=50 and the angular searching step size $\delta$=0.1. The significance
level in PCM is set to be $\epsilon$=0.1. Our approach was tested on an AMD
Ryzen 5950X processor with a frequency of 3.4GHz. The computational time was
cumulated over all three robots based on one pose update. The pose is
expressed in double precision (8 Bytes) in our implementation. We show the
total data communicated among all three robots at different updating
frequencies (i.e., 10Hz and 1Hz) of DPGO. As can be seen from this table, the
total computational time required for the three modules is 81ms, which we
believe is fast enough for the online localization and control of mobile
robots. The cost of communication is small, which can be easily deployed in
most robotics applications. As can be also seen from this table, performing
DPGO at a higher frequency costs more communication resources, due to the
frequent exchange of the poses among the robots.

TABLE IV: Computational time (millisecond, ms) and total communication cost (Megabytes, MB) for our approach. Time consumption | Communication
---|---
Pose est. | Dist. PCM | DPGO | Total | 10Hz | 1Hz
39ms | 24ms | 18ms | 81ms | 13.8MB | 4.3MB

##  V Conclusions

We proposed an approach for distributed SLAM using UWB ranging and odometry
measurements for a group of robots in unknown environments. The proposed
method determines the relative pose between robots based on short-term UWB
ranging and odometry measurements. In addition, we presented a distributed
pose graph optimization with pairwise consistency check for trajectory
estimation. The proposed approach was verified with three robots in an indoor
environment with a size of 10m$\times$12m. The results showed that the
proposed distributed ranging SLAM achieved a localization accuracy of 0.45m in
translation and 4.01∘ in rotation. Our approach provides a solution to the
robotics community for the localization of a team of robots without any
knowledge about the infrastructure. In future, we would like to generalize our
approach to 3D, which requires the coarse search of the initial guess in 3D
polar coordinate system. The comparison to the quadratically constrained
quadratic programming-based solution [9] will be also included in our future
work.

## References

  * [1] R. Liu, C. Yuen, T.-N. Do, M. Zhang, Y. L. Guan, and U.-X. Tan, “Cooperative positioning for emergency responders using self imu and peer-to-peer radios measurements,” _Information Fusion_ , vol. 56, pp. 93 – 102, 2020.
  * [2] M. Afrin, J. Jin, A. Rahman, A. Rahman, J. Wan, and E. Hossain, “Resource allocation and service provisioning in multi-agent cloud robotics: A comprehensive survey,” _IEEE Communications Surveys and Tutorials_ , vol. 23, no. 2, pp. 842–870, 2021.
  * [3] Y. Tian, K. Khosoussi, D. M. Rosen, and J. P. How, “Distributed certifiably correct pose-graph optimization,” _IEEE Transactions on Robotics_ , vol. 37, no. 6, pp. 2137–2156, 2021.
  * [4] Y. Tian, Y. Chang, F. H. Arias, C. Nieto-Granda, J. P. How, and L. Carlone, “Kimera-multi: Robust, distributed, dense metric-semantic slam for multi-robot systems,” _IEEE Transactions on Robotics_ , pp. 1–17, 2022.
  * [5] A. Yassin, Y. Nasser, M. Awad, A. Al-Dubai, R. Liu, C. Yuen, R. Raulefs, and E. Aboutanios, “Recent advances in indoor localization: A survey on theoretical approaches and applications,” _IEEE Communications Surveys Tutorials_ , vol. 19, no. 2, pp. 1327–1346, Secondquarter 2017.
  * [6] R. Liu, C. Yuen, T. N. Do, D. Jiao, X. Liu, and U-X. Tan, “Cooperative relative positioning of mobile users by fusing imu inertial and uwb ranging information,” in _the 2017 IEEE International Conference on Robotics and Automation (ICRA 2017)_ , Singapore, May 29-June 3 2017, pp. 5623–5629.
  * [7] N. Rogel, D. Raphaeli, and O. Bialer, “Time of arrival and angle of arrival estimation algorithm in dense multipath,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 5907–5919, 2021.
  * [8] M. Shalihan, R. Liu, and C. Yuen, “NLOS ranging mitigation with neural network model for UWB localization,” _CoRR_ , vol. abs/2206.09607, 2022\. [Online]. Available: https://arxiv.org/abs/2206.09607
  * [9] T. H. Nguyen and L. Xie, “Relative transformation estimation based on fusion of odometry and uwb ranging data,” _CoRR_ , vol. abs/2202.00279, 2022. [Online]. Available: https://arxiv.org/abs/2202.00279
  * [10] K. Ebadi, Y. Chang, M. Palieri, A. Stephens, A. Hatteland, E. Heiden, A. Thakur, N. Funabiki, B. Morrell, S. Wood, L. Carlone, and A.-a. Agha-mohammadi, “Lamp: Large-scale autonomous mapping and positioning for exploration of perceptually-degraded subterranean environments,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , Paris, France, May 31-August 31 2020, pp. 80–86.
  * [11] Y. Zhu, C. Zheng, C. Yuan, X. Huang, and X. Hong, “Camvox: A low-cost and accurate lidar-assisted visual slam system,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , Xi’an, China, May 30-June 5 2021, pp. 5049–5055.
  * [12] T. H. Nguyen, T.-M. Nguyen, and L. Xie, “Range-focused fusion of camera-imu-uwb for accurate and drift-reduced localization,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 1678–1685, 2021.
  * [13] R. Liu, Y. He, C. Yuen, B. P. L. Lau, R. Ali, W. Fu, and Z. Cao, “Cost-effective mapping of mobile robot based on the fusion of uwb and short-range 2-d lidar,” _IEEE/ASME Transactions on Mechatronics_ , pp. 1–1, 2021.
  * [14] Y. Cao, C. Yang, R. Li, A. Knoll, and G. Beltrame, “Accurate position tracking with a single uwb anchor,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , Paris, France, May 31-August 31 2020, pp. 2344–2350.
  * [15] N. Funabiki, B. Morrell, J. Nash, and A.-a. Agha-mohammadi, “Range-aided pose-graph-based slam: Applications of deployable ranging beacons for unknown environment exploration,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 1, pp. 48–55, 2021.
  * [16] Z. Cao, R. Liu, C. Yuen, A. Athukorala, B. K. Kiat Ng, M. Mathanraj, and U.-X. Tan, “Relative localization of mobile robots with multiple ultra-wideband ranging measurements,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , Prague, Czech Republic, September 27-October 1 2021, pp. 5857–5863.
  * [17] F. Molina Martel, J. Sidorenko, C. Bodensteiner, M. Arens, and U. Hugentobler, “Unique 4-dof relative pose estimation with six distances for uwb/v-slam-based devices,” _Sensors_ , vol. 19, no. 20, 2019.
  * [18] T. H. Nguyen, T.-M. Nguyen, and L. Xie, “Flexible and resource-efficient multi-robot collaborative visual-inertial-range localization,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 928–935, 2022.
  * [19] S. Choudhary, L. Carlone, C. Nieto, J. Rogers, H. I. Christensen, and F. Dellaert, “Distributed trajectory estimation with privacy and communication constraints: A two-stage distributed gauss-seidel approach,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_ , Stockholm, Sweden, May 16-May 21 2016, pp. 5261–5268.
  * [20] H. Xu and S. Shen, “BDPGO: balanced distributed pose graph optimization framework for swarm robotics,” _CoRR_ , vol. abs/2109.04502, 2021. [Online]. Available: https://arxiv.org/abs/2109.04502
  * [21] P.-Y. Lajoie, B. Ramtoula, Y. Chang, L. Carlone, and G. Beltrame, “Door-slam: Distributed, online, and outlier resilient slam for robotic teams,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1656–1663, 2020\.
  * [22] R. Kümmerle, G. Grisetti, H. Strasdat, K. Konolige, and W. Burgard, “G2o: A general framework for graph optimization,” in _2011 IEEE International Conference on Robotics and Automation_ , Shanghai, China, May 9-13 2011, pp. 3607–3613.
  * [23] F. Dellaert, “Factor graphs and gtsam: A hands-on introduction,” Technical Report number GT-RIM-CP&R-2012-002, 2012.
  * [24] S. Agarwal _et al._ , “Ceres solver,” http://ceres-solver.org.
  * [25] J. G. Mangelson, D. Dominic, R. M. Eustice, and R. Vasudevan, “Pairwise consistent measurement set maximization for robust multi-robot map merging,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , Brisbane, QLD, Australia, May 21-May 25 2018, pp. 2916–2923.
  * [26] R. Smith, M. Self, and P. Cheeseman, “Estimating uncertain spatial relationships in robotics,” in _1987 IEEE International Conference on Robotics and Automation_ , Raleigh, NC, USA, March 31-April 3 1987, pp. 850–850.
  * [27] T. Ziegler, M. Karrer, P. Schmuck, and M. Chli, “Distributed formation estimation via pairwise distance measurements,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 3017–3024, 2021.