latexYou have requested package capbtabboxtable[][]

11institutetext: University of California Riverside
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# FRAPPE: Fast Rank Approximation with Explainable Features for Tensors

William Shiao  0000-0001-5813-2266 Evangelos E. Papalexakis
0000-0002-3411-8483

###### Abstract

Tensor decompositions have proven to be effective in analyzing the structure
of multidimensional data. However, most of these methods require a key
parameter: the number of desired components. In the case of the
CANDECOMP/PARAFAC decomposition (CPD), the ideal value for the number of
components is known as the canonical rank and greatly affects the quality of
the decomposition results. Existing methods use heuristics or Bayesian methods
to estimate this value by repeatedly calculating the CPD, making them
extremely computationally expensive. In this work, we propose FRAPPE, the
first method to estimate the canonical rank of a tensor without having to
compute the CPD. This method is the result of two key ideas. First, it is much
cheaper to generate synthetic data with known rank compared to computing the
CPD. Second, we can greatly improve the generalization ability and speed of
our model by generating synthetic data that matches a given input tensor in
terms of size and sparsity. We can then train a specialized single-use
regression model on a synthetic set of tensors engineered to match a given
input tensor and use that to estimate the canonical rank of the tensor—all
without computing the expensive CPD. FRAPPE is over 24× faster than the best-
performing baseline, and exhibits a $10\%$ improvement in MAPE on a synthetic
dataset. It also performs as well as or better than the baselines on real-
world datasets.

ank estimation Tensor decomposition Self-supervised CANDECOMP/PARAFAC

###### Keywords:

R

##  1 Introduction

Tensors are widely used in the field of computer science and can be used to
model a variety of datasets. One popular family of methods that can be applied
to tensors is tensor decompositions. The performance of these methods
typically depends on a key parameter—the number of components in the
decomposition, also known as the rank of the decomposition. If our
decomposition rank is too low, we may fail to capture key information, but if
it is too high, we may capture too much noise. In the case of the popular
CANDECOMP/PARAFAC decomposition, (CPD) [9], the ideal value for the rank of
the decomposition is known as the canonical rank. However, finding this is a
challenging task. A similar problem, finding the exact rank of a tensor (the
minimum number of CPD components required to represent the tensor exactly), is
known to be NP-hard over the set of rational numbers [11, 12].

As such, various heuristic-based approaches have been proposed to solve this
task. One such method, proposed by Bro and Kiers [5], uses the core
consistency diagnostic (CORCONDIA) to evaluate the “appropriateness” of a
PARAFAC model. AutoTen [20] uses CORCONDIA to select a suitable number of
components automatically. Normalized SVD (NSVD) [24] also attempts to provide
a heuristic for the rank of a tensor, but it is not fully automatic and
currently requires human supervision. Fernandes et al. [8] propose NORMO,
which computes the CPD of a tensor at various ranks, looks at the mean
correlation between each of its components, and selects the maximum number of
components with no redundancy between the components.

Another set of approaches is the Bayesian approaches to this problem. Zhao et
al. [27] attempt to tackle this problem by using Bayesian inference to
automatically determine the decomposition’s rank. Mørup and Hansen [17] use a
Bayesian framework with Automatic Relevance Determination (ARD) to estimate
the rank of a Tucker decomposition. This can easily be adapted since the CPD
can be formulated as a restricted Tucker decomposition [25]. Zhou et al. [28]
compute the CPD and use a convolutional neural network (CNN) on the factors in
an effort to determine the canonical rank of the tensor. However, they
primarily focus on and evaluate their method within the image domain.

One common theme across all of these methods is that they all require
exhaustively computing the CPD multiple times across all the candidate ranks.
This is computationally expensive and impractical for large tensors, so we
would like to avoid this step. However, doing so poses several challenges: (a)
it is difficult to directly train a supervised machine learning model due to
the scarcity of tensors with known canonical rank (which can typically only be
determined via manual human inspection), and (b) it is difficult to develop a
model capable of handling tensors from various domains and of various sizes.

To solve (a), we propose using synthetic tensors, which we can create with
known canonical rank. To solve (b), we propose using hand-picked size-agnostic
features extracted from each tensor. However, we find that supervised
regression models trained using that data (see Sec. 4.2) fail to generalize
well to real-world tensors. To overcome this limitation, we propose FRAPPE: a
fully self-supervised model that uses synthetic data specifically tailored for
a given input tensor. FRAPPE outperforms existing baselines on synthetic data
and is on par with the best baselines on real-world data. It is also
${\sim}$24× faster than the best-performing baseline (NORMO).

We evaluate the performance of both models on three classes of tensors:
synthetic tensors, real tensors that have been extensively studied in existing
literature to determine their rank, and the weight tensor of a convolutional
neural network (to find the ideal CPD rank for compressing the weights while
maintaining good classification accuracy). Finally, we also use these models
to analyze which features are the best and worst predictors of a tensor’s
rank.

Our contributions include the following:

  * •

Novel method: We propose FRAPPE—the first method to estimate the canonical
rank of generic tensor without having to compute the CPD.

  * •

Speed: We show that FRAPPE offers a ${\sim}$24× improvement in evaluation time
over the best-performing baseline method and is ${\sim}$1.6× faster than the
fastest baseline.

  * •

Extensive Evaluation: We evaluate our methods and baseline methods on
synthetically generated tensors, real tensors with known rank, and the weight
tensor of a convolutional neural network. We show that FRAPPE outperforms the
best baseline by about 10% MAPE on a synthetic dataset and performs as well as
or better than baselines on real-world datasets.

  * •

Explainability: We examine the importance of each of our hand-picked features
to see which features contribute the most to classification results (see Fig.
5).

  * •

Reproducible: We provide the Python code for our method online at
https://anonymous.4open.science/r/frappe-4D01/ (will be posted to publicly
GitHub if accepted). We also provide the code to reproduce our baseline
results.

Table 1: Comparison of different tensor rank estimation methods. All speeds
are relative to the evaluation time of the slowest baseline method. Exact
evaluation times are shown in Fig. 4. †: We did not benchmark against TRN-PD
because it is designed for image data and image denoising, which is not the
focus of this work.

| Proposed | Baseline Methods
---|---|---
|  FRAPPE |  NORMO [8] |  ARD-Tucker [17] |  AutoTen [20] |  TRN-PD [28]
Avoids Computing CPD | ✓ | ✗ | ✗ | ✗ | ✗
Avoids Human Supervision | ✓ | ✓ | ✓ | ✓ | ✓
Avoids Strictly Labelled Dataset | ✓ | ✓ | ✓ | ✓ | ✗
Explainable | ✓ |  $\circ$111NORMO is partially explainable in that higher correlation in the factors indicates a bad choice for a given rank. However, NORMO does not provide any indication as to which attributes in the original tensor contribute to higher/lower predictions. | ✗ | ✗ | ✗
Works for General Tensors | ✓ | ✓ | ✓ | ✓ | ✗222TRN-PD is designed for images and was originally evaluated for CPD-based denoising.
Eval. Speed | $24\times$ | $1\times$ | $12\times$ | $18\times$ | †

##  2 Background

Here we provide an overview of the notation and background for this work. We
denote a set, tensor, matrix, vector, and scalar as $\mathcal{X}$,
$\boldsymbol{\mathscr{X}},\textbf{X}$, x, $x$, respectively. We write the
Frobenius norm of X as $\left\lVert\textbf{X}\right\rVert_{F}$. Let
$\boldsymbol{\mathscr{T}}_{(i,j,k)}$ denote the value in the $i$-th, $j$-th,
and $k$-th entry along the 1st, 2nd, and 3rd modes, respectively, of the 3rd
order tensor $\boldsymbol{\mathscr{T}}$. Let $(:)$ refer to all possible
indexes in the given mode when used in tensor indexing (e.g.,
$\boldsymbol{\mathscr{T}}=\boldsymbol{\mathscr{T}}_{(:,:,:)}$ for a third-
order tensor $\boldsymbol{\mathscr{T}}$).

We also define several functions. Let
$\textsc{Zeros}{}(\boldsymbol{\mathscr{T}})$ return the number of zero values
in a tensor $\boldsymbol{\mathscr{T}}$, and $\textsc{Rand}(i,j)$ return a
random matrix of size $i\times j$. Finally, let
$\boldsymbol{\mathscr{T}}\sim\mathcal{N}_{i\times j\times k}(\mu,\sigma)$
represent sampling a tensor $\boldsymbol{\mathscr{T}}$ of size $i\times
j\times k$ from a normal distribution with a mean of $\mu$ and standard
deviation of $\sigma$. $\circ$ represents the tensor outer product.

###  2.1 CANDECOMP/PARAFAC

One of the most common tensor decompositions is the canonical polyadic or
CANDECOMP/PARAFAC decomposition (CPD) [10, 6]. The CPD decomposes a tensor
$\boldsymbol{\mathscr{T}}$ into the sum of $R$ outer products of vectors.
Formally, for a third-order tensor:

$\boldsymbol{\mathscr{T}}\approx\sum^{R}_{r=1}\mathbf{a}_{r}\circ\mathbf{b}_{r}\circ\mathbf{c}_{r}$
(1)

We can then attempt to iteratively minimize the Frobenius norm of the error:
$\left\lVert\boldsymbol{\mathscr{T}}-\sum^{R}_{r=1}\mathbf{a}_{r}\circ\mathbf{b}_{r}\circ\mathbf{c}_{r}\right\rVert_{F}$.
A common method for this is to use the alternating least squares algorithm. By
convention, we refer to the factors as matrices
$\mathbf{A},\mathbf{B},\mathbf{C}$, which consist of the corresponding vectors
horizontally stacked. For example, the $i$-th column of $\mathbf{A}$ would be
$\mathbf{a}_{i}$. A key component of this method is $R$—the rank of the
decomposition. As described above, choosing a good value for $R$ is extremely
difficult. The ideal value for $R$ is the canonical rank of
$\boldsymbol{\mathscr{T}}$, and estimating that value is the focus of this
work.

##  3 Problem Formulation

We tackle the problem of predicting the canonical rank of a tensor. Formally,
the problem can be stated as follows:

Problem Definition Given a tensor
$\boldsymbol{\mathscr{X}}=\boldsymbol{\mathscr{L}}+\boldsymbol{\mathscr{N}}$,
where $\boldsymbol{\mathscr{L}}$ is a tensor with low rank and
$\boldsymbol{\mathscr{N}}$ is tensor of high-rank noise, find the canonical
rank $R$ of $\boldsymbol{\mathscr{X}}$ such that we minimize
$\left\lVert\boldsymbol{\mathscr{L}}-\sum^{R}_{r=1}\mathbf{a}_{r}\circ\mathbf{b}_{r}\circ\mathbf{c}_{r}\right\rVert_{F}$
where $\textbf{A},\textbf{B},\textbf{C}$ are factors from the CPD of
$\boldsymbol{\mathscr{X}}$.

Note that $R$ is different from the exact rank of the tensor (the minimum
number of components needed to exactly represent $\boldsymbol{\mathscr{X}}$),
which is out of the scope of this work.

Another secondary goal of our method is to identify which features are
important in predicting the rank of a tensor, which could help us
develop/improve rank prediction methods in the future. We define this task as
follows:

Subproblem: Explainability Given a tensor $\boldsymbol{\mathscr{X}}$ and a
predicted rank $\hat{R}$, compute the importance scores t of different
attributes of $\boldsymbol{\mathscr{X}}$ in the regression result.

##  4 Proposed Method

We first describe which hand-picked features we extract from the tensor and
the reasoning behind selecting each one. We then train a regression model on
the dataset and show that it has some key limitations: (a) it performs well on
the synthetic dataset but fails to generalize well to real data, and (b) it
requires a relatively expensive training step. We then propose our method,
FRAPPE, which addresses both of these issues.

###  4.1 Tensor Features

Given a tensor $\boldsymbol{\mathscr{T}}\in\mathbb{R}^{I\times J\times K}$, we
extract 112 features (4D or higher tensors require some adjustment, as
described in Sec. 4.1.3):

  * •

Tensor Dimensions (#1–3): $I,J,K$. We include these features so that the
underlying model can adjust the expectations for other parameters relative to
the size of the tensor. This allows our model to generalize to different input
sizes.

  * •

Non-zeros (#4–13): The total number of non-zero values in the tensor, as well
as the minimum, median, and maximum number of non-zero values over all
possible slices of the tensor. That is, we extract the number of non-zeros
from $\boldsymbol{\mathscr{T}}_{(i,:,:)}\forall i\in[1,I]$,
$\boldsymbol{\mathscr{T}}_{(:,j,:)}\forall j\in[1,J]$, and
$\boldsymbol{\mathscr{T}}_{(:,:,k)}\forall k\in[1,K]$. The intuition behind
this feature is that, generally, the rank of a tensor increases as it becomes
more sparse. Measuring this for all possible slices allows us to essentially
measure the sparsity along each dimension.

  * •

Ranks of Slices (#14-22): The minimum, median, and maximum rank over all
possible slices of the tensor. The intuition behind this feature is that if
slices along any given two dimensions are generally of low rank (implying that
fewer components are needed in the SVD), fewer components may also be required
in the CPD.

  * •

Thresholded Non-zeros (#23–103): The minimum, median, and maximum number of
values in the normalized (we use min-max normalization) tensor less than or
equal to each of the following thresholds: $[0.1,0.2,\ldots,0.9]$ across all
possible slices of the tensor. This is an adaptation of the number of non-
zeros (features #4-13) for dense tensors and allows us to find the number of
low-value entries, which could be a result of noise and be of lesser
importance to the canonical rank of the tensor.

  * •

Correlations (#104–112): The minimum, median, and maximum correlation over all
possible pairs of slices of the tensor. The intuition behind this feature is
that higher correlation between slices generally indicates that fewer
components would be needed along a given axis, resulting in a lower canonical
rank.

Note that the number of features is constant, regardless of the size of the
tensor. This allows it to easily scale to large tensors. We then train a
LightGBM [13] model on these features, with the goal of predicting the rank of
the target tensor.

####  4.1.1 Why Correlation?

It may not be immediately apparent why we choose to use the correlation
between different slices. We choose to do so because the distribution of
pairwise correlation of vectors in a matrix is a good indicator of its rank.
We show this theoretically (under some assumptions) below. We hypothesize that
this remains true in the tensor case; while difficult to formally prove, we
also show that this feature is effective in our model in Fig. 5.

###### Proof

We first define the rank of a matrix ${\mathbf{M}}$ as $k$, where $k$ is the
largest value such that the $k$-th singular value $\sigma_{k}>\epsilon$. Let
us assume that any two vectors $\textbf{x},\textbf{y}$ has an outer product
$\textbf{x}\textbf{y}^{T}$ that can be written as
$\textbf{a}\textbf{b}^{T}+\textbf{N}$, where
$\textbf{N}=\gamma\sigma_{1}\textbf{u}_{1}\textbf{v}_{1}^{T}+\gamma^{2}\sigma_{2}\textbf{u}_{2}\textbf{v}_{2}^{T}+\ldots+\gamma^{k}\sigma_{k}\textbf{u}_{k}\textbf{v}_{k}^{T}$
such that $\sigma_{i}<\epsilon\ \forall i$, where $\epsilon$ is an arbitrarily
small number and $\textsc{Corr}(\textbf{x},\textbf{y})\rightarrow
1\implies\gamma\rightarrow 0$. Let
$c=\mathbb{E}\left[\textsc{Corr}({\mathbf{M}}^{(i)},{\mathbf{M}}^{(j)})\right]\forall({\mathbf{M}}^{(i)},{\mathbf{M}}^{(j)})\in\mathscr{C}$.
Let $\mathbf{M}\in{\mathbb{R}}^{m\times n}$, and ${\mathbf{M}}^{(i)}$ be the
$i$-th column of ${\mathbf{M}}$. For the sake of simplicity, let us also
assume that a set $\mathscr{C}$ contains all $\frac{n}{2}$ non-overlapping
pairs of columns, and all other pairs have a correlation of 0. Then, we can
write ${\mathbf{M}}$ as
$\textbf{a}_{1}\textbf{b}_{1}^{T}+\textbf{a}_{2}\textbf{b}_{2}^{T}+\ldots+\textbf{a}_{n/2}\textbf{b}_{n/2}^{T}+\textbf{N}_{1}+\textbf{N}_{2}+\ldots+\textbf{N}_{n/2}$.
Let
$\textbf{N}_{\Sigma}=\textbf{N}_{1}+\textbf{N}_{2}+\ldots+\textbf{N}_{n/2}$.
Since we have at most $\frac{n}{2}$ duplicate components, the maximum value
for a singular value of N is $\frac{n}{2}\gamma_{c}\epsilon$, which is still
less than $\epsilon$ for $\gamma<\frac{2}{n}$, and subsequently true for
sufficiently large $c$. Therefore, the rank of ${\mathbf{M}}$ decreases as our
mean pairwise correlation $c$ increases.

####  4.1.2 Computational Complexity.

Features #1–13 and #23–103 are extremely quick to calculate, requiring only a
single pass over the tensor. Features #14–22 and #104–112 are significantly
more expensive since they require computing the SVD and the correlation,
respectively, over various slices of the tensor. However, these operations are
still much faster than the baseline methods, which require repeatedly
computing the full CPD over the whole tensor. Fig. 4 compares the speed of our
method against the baselines.

Let $n$ $=\max(I,J,K)$. Features #1-13 take a linear amount of time to
calculate with respect to the number of non-zero values in the tensor, with a
worst-case of ${\mathcal{O}}(n^{3})$ (in the case of a dense tensor). Features
#14-22 take ${\mathcal{O}}((I+J+K)s)={\mathcal{O}}(ns)$ time, where $s$ is the
runtime of sparse SVD on a slice of our tensor. In the case of a fully dense
tensor, this takes ${\mathcal{O}}((I+J+K)n^{3})={\mathcal{O}}(n^{4})$ time to
calculate.

Features #23–103 require the tensor to be min-max normalized first, which can
be done in ${\mathcal{O}}(n^{3})$ time. It also takes ${\mathcal{O}}(n^{3})$
to count the number above/below each threshold. Features #104–112 take
${\mathcal{O}}(n^{2})$ time to calculate over $\frac{n(n-1)}{2}$ pairs,
resulting in a runtime of ${\mathcal{O}}(n^{4})$.

Overall, these features take a maximum of ${\mathcal{O}}(n^{4})$ time to
calculate for dense tensors, where $n$ is the size of the largest dimension in
the tensor. It is worth noting that it would take significantly less time for
sparse tensors (e.g., a multi-view graph, knowledge base, word coccurrence).
This is also a massive improvement compared to the runtime of computing the
CPD once, let alone the repeated runs required by existing methods. We show a
comparison of the runtimes below in Fig. 4.

####  4.1.3 Extension to Higher-Order Tensors

While the features described above work well for 3rd-order tensors, they
require some modification for 4th or higher-order tensors. For example, the
rank features are computed over all slices of the tensor. In the 3rd order
case, these slices are matrices, but in the 4th (or higher) order case, these
slices are 3rd (or higher) order sub-tensors. Finding the rank of a matrix can
be done in polynomial time via the SVD, but finding the rank of a 3D tensor is
the problem we are tackling in this paper. Due to this issue and the
relatively low importance of the rank features in Fig. 5, we omit the rank
features for the 4D model and find that it still performs fairly well. This is
shown empirically in our experiments on a 4D AlexNet weight tensor below in
Sec. 5.4. As with existing work in canonical rank estimation [8, 20], we focus
on primarily on 3rd-order tensors and reserve detailed evaluation for future
work.

###  4.2 A First Attempt

We can now train a regression model using the features defined in Sec. 4.1
above. To do this, we first generate a synthetic dataset with tensors of known
rank (see Sec. 5.2). We then train a LightGBM [13] model on this synthetic
dataset and evaluate the model on each of the evaluation datasets in Sec. 5
below. However, we find that the LightGBM model is unable to generalize well
beyond the synthetic datasets and incorrectly predicts the rank on all of the
real-world datasets (described in Sec. 5.3). This is likely because real-world
tensors come from a variety of different domains.

Another limitation of this approach is that we have to train our model on a
very large dataset to ensure that we account for different ranks, sparsities,
and noise levels. While this makes predicting the rank of a new tensor very
quick, it greatly increases the amount of time required to train a regression
model. To address this, we propose FRAPPE.

Algorithm 1 Pseudocode for FRAPPE: a self-supervised algorithm to estimate the
rank of a given tensor. Note that comments precede the line that they
describe.

1: An input tensor $\boldsymbol{\mathscr{T}}\in\mathbb{R}^{i\times j\times
k}$, $\rho$ (the maximum rank to consider), $\tau$ (the number of times to
sample from each rank), and $\alpha$ (the maximum amount of noise to add in
synthetic samples).

2: $\hat{R}$, an estimate for the rank of the input tensor
$\boldsymbol{\mathscr{T}}$.

3: $\hat{R}$, an estimate for the rank of the input tensor
$\boldsymbol{\mathscr{T}}$. function FRAPPE($\boldsymbol{\mathscr{T}}$,
$\rho,\ \tau,\ \alpha$)

4: $\texttt{data}\leftarrow[\ ]$

5: $\texttt{ranks}\leftarrow[\ ]$

6: $\triangleright$ Calculate the sparsity of our tensor $\triangleleft$

7: $s\leftarrow\frac{\textsc{zeros}(\boldsymbol{\mathscr{T}})}{i\times j\times
k}$

8: $s\leftarrow\frac{\textsc{zeros}(\boldsymbol{\mathscr{T}})}{i\times j\times
k}$ for $x\in[1,\ \tau]$ do

9: for $x\in[1,\ \tau]$ do for $r\in[1,\rho]$ do

10: ${\mathbf{A}}\leftarrow\textsc{Sparsify}(\textsc{Rand}(i,r),s)$

11: ${\mathbf{B}}\leftarrow\textsc{Sparsify}(\textsc{Rand}(j,r),s)$

12: ${\mathbf{C}}\leftarrow\textsc{Sparsify}(\textsc{Rand}(k,r),s)$

13: $\triangleright$ Construct a synthetic tensor $\boldsymbol{\mathscr{S}}$
from ${\mathbf{A}},{\mathbf{B}},{\mathbf{C}}$ $\triangleleft$

14:
$\boldsymbol{\mathscr{S}}\leftarrow\sum^{r}_{y=1}\mathbf{a}_{y}\circ\mathbf{b}_{y}\circ\mathbf{c}_{y}$

15: $\triangleright$ Select a random noise ratio. $\triangleleft$

16: $n\sim[0,\alpha]$

17: $\triangleright$ Sample an $i\times j\times k$ Gaussian noise tensor
$\triangleleft$

18: $\boldsymbol{\mathscr{N}}\sim\mathcal{N}_{i\times j\times k}(0,1)$

19: $\triangleright$ Scale the noise and add it to $\boldsymbol{\mathscr{S}}$
$\triangleleft$

20:
$\boldsymbol{\mathscr{S}}\leftarrow\boldsymbol{\mathscr{S}}+(n\frac{\left\lVert\boldsymbol{\mathscr{S}}\right\rVert_{F}}{\left\lVert\boldsymbol{\mathscr{N}}\right\rVert_{F}})\boldsymbol{\mathscr{N}}$

21: $\triangleright$ Extract features from $\boldsymbol{\mathscr{S}}$, as
described in Sec. 4.1 $\triangleleft$

22: $\texttt{data}.append(\textsc{ExtractFeatures}(\boldsymbol{\mathscr{S}}))$

23: $\texttt{ranks}.append(r)$ $\texttt{ranks}.append(r)$
$\texttt{ranks}.append(r)$ $\texttt{ranks}.append(r)$

24: Train a LightGBM regression model on data and ranks

25: Evaluate the model on $\textsc{ExtractFeatures}(\boldsymbol{\mathscr{T}})$
to get our prediction $\hat{r}$

26: return $\hat{R}$ return $\hat{R}$

1: An input matrix ${\mathbf{M}}\in\mathbb{R}^{a\times b}$ and the desired
sparsity level $s$.

2: The sparsified matrix.

3: The sparsified matrix. function Sparsify(${\mathbf{M}},s$)

4: Select $s\times a\times b$ random indices of ${\mathbf{M}}$

5: Set each of these indices to 0

6: return ${\mathbf{M}}$ return ${\mathbf{M}}$

###  4.3 FRAPPE

The primary issue we need to solve with FRAPPE compared to the approach
described in Sec. 4.2 is the poor generalization of the naive regression
model. This poor generalization is fundamentally caused by a distribution
difference between the training data and the input tensor. We can address this
issue by ensuring that the synthetic training data closely mimics a given
input tensor. We do this by attempting to match the size and sparsity of the
input tensor when generating the synthetic dataset. We can then train a
specialized LightGBM model for each tensor being considered. This specialized
model only performs evaluation on the input tensor once before being
discarded.

Generating training data engineered to mimic an input tensor also solves the
issue of long training times since less data is required to achieve the same
results. This brings down training times significantly, allowing us to train a
disposable regression model for each input tensor while maintaining a speed
advantage over the other baseline models. Note that the “training” happens at
evaluation time, and no external training step is required. FRAPPE is formally
described in detail in Algorithm 1.

Figure 1: Diagram of FRAPPE. We first generate synthetic tensors based on the
size and sparsity of the input tensor, extract features, train the model, and
use that model to predict the rank of the original input tensor.

##  5 Experimental Evaluation

Evaluating the effectiveness of our methods is difficult because there are
very few tensors with known rank so we cannot solely rely on real-world data.
As such, we choose to evaluate our methods on three different types of data:

  1. 1.

Synthetic Tensors: Predicting the rank of tensors that we generate to have
known rank (described in Sec. 5.2).

  2. 2.

Real-world Tensors: Predicting the rank of real-world tensors that have been
previously studied and have known rank (described in Sec. 5.3).

  3. 3.

CNN Weights: Compressing the 4D weights of an AlexNet [16] model’s final layer
(described in Sec. 5.4).

We compare the performance and runtime of our methods to other current state-
of-the-art methods. All of the timing experiments are conducted on a
c3-highcpu-176 Google Cloud Platform instance with 176 CPU cores and 352 GB of
RAM.

###  5.1 Baseline Methods

We compare FRAPPE and FRAPPE against three baselines: NORMO [8], AutoTen [20]
and ARD-Tucker [17]. These methods require a maximum rank to consider. Since
this greatly affects the runtime of the algorithm, we used double the true
rank of the tensor for the two algorithms. This theoretically bounds the
maximum error, but setting it higher would greatly increase both the runtime
and computational cost of running the baselines. We also use this same bound
for FRAPPE.

We had to slightly modify AutoTen to only use the Frobenius-norm-based
evaluation because we have negative values in our dataset, which causes the
KL-divergence-based evaluation to fail. This slightly reduces the accuracy of
AutoTen but greatly improves its runtime. We also had to modify NORMO and
AutoTen to add support for 4D tensors in order to run it on our CNN weight
tensor. We chose to use the suggested value of $\delta=0.7$ for NORMO.

Table 2: Performance of the rank estimation methods on real-world datasets.
The first column is the known rank of the datasets, and correct estimates are
bolded. We can see that FRAPPE and AutoTen generally perform the best overall.

Dataset | Rank | FRAPPE | NORMO | AutoTen | ARD
---|---|---|---|---|---
amino_acid [3] | 3 [3, 14] | 3 | 4 | 3 | 3
dorrit [23] | 4 [23] | 4 | 4 | 4 | 4
sugar [4] | 3 [4] | 3 | 4 | 2 | 3
reality_mining [7] | 6?333The number of communities in this dataset was explored in [19] and estimated to be 6. The number of communities loosely corresponds to the canonical rank of the tensor, but this is not a guarantee like with the fluorescence datasets. [19] | 9 | 2 | 2 | 18
enron [22] | 4/7? [2, 21] | 3 | 1 | 4 | 12

###  5.2 Synthetic Dataset

To generate a synthetic dataset, we first generate the CPD factor matrices
$\mathbf{A},\mathbf{B},\mathbf{C}$ by drawing values from a uniform random
distribution. We ensure that we generate a diverse dataset by generating a set
percentage of tensors from each of the following subsets (each of which
corresponds to a common real-world tensor):

  1. 1.

Sparse tensors with sparse noise: Tensors that have sparsity introduced in the
factors and Gaussian noise added to non-zero elements in the tensor. Example:
a weighted time-evolving graph.

  2. 2.

Sparse tensors with dense noise: Tensors that have sparsity introduced in the
factors and Gaussian noise added to all elements in the tensor.
Example: a time-evolving grid of sensor readings.

  3. 3.

Dense tensors: Tensors that have fully dense tensors and dense noise added.
Example: a tensor of fluorometer readings.

Each of these cases represents a different type of tensor that can be found in
real-world datasets. We add 2%–10% noise in each case and scale the amount of
noise by the norm of the noise tensor compared to the generated tensor. More
formally, where $\alpha$ is the noise ratio, $\boldsymbol{\mathscr{N}}$ is our
noise tensor and $\boldsymbol{\mathscr{G}}$ is our generated tensor:
$\boldsymbol{\mathscr{G}}=\boldsymbol{\mathscr{G}}+\alpha\frac{\left\lVert\boldsymbol{\mathscr{G}}\right\rVert_{F}}{\left\lVert\boldsymbol{\mathscr{N}}\right\rVert_{F}}\boldsymbol{\mathscr{N}}$.
This is similar to the noise addition process in the Tensor Toolbox’s [1]
create_problem function.

$\displaystyle{0}$$\displaystyle{20}$$\displaystyle{40}$$\displaystyle{60}$$\displaystyle{80}$$\displaystyle{100}$Rank$\displaystyle{0.1}$$\displaystyle{0.2}$$\displaystyle{0.3}$$\displaystyle{0.4}$$\displaystyle{0.5}$$\displaystyle{0.6}$$\displaystyle{0.7}$AccuracyFRAPPENORMOAutoTenOriginal
Accuracy

Figure 2: Plot of AlexNet performance with the final CNN layer compressed at
various ranks. Also shown are the rank estimates produced by FRAPPE and the
baselines. NORMO predicts a very low value and AutoTen predicts a very high
value. FRAPPE and FRAPPE pick what appear to be reasonable values. ARD-Tucker
did not terminate in time and is not included.

###  5.3 Real-World Datasets

We evaluate the effectiveness of our methods and baseline methods on datasets
of known rank. These include three fluorescence datasets: amino_acid [3],
dorrit [23], and sugar [4]. These datasets contain naturally tri-linear data
with a known number of chemicals, each of which corresponds to a component in
the CPD. We also evaluate our methods on the reality_mining [7] and enron [22]
datasets. The canonical rank of these tensors is less clear, as we have to
rely on existing literature that has manually analyzed the data and determined
the approximate rank.

From the results displayed in Table 2, we can see FRAPPE is able to perfectly
estimate the rank of the fluorescence datasets. ARD is also able to correctly
estimate the rank of the datasets. NORMO and AutoTen perform fairly well but
are slightly off. Only AutoTen is able to successfully estimate the rank for
the enron dataset, but FRAPPE is also very close. Unfortunately, we are
limited in the amount of data we can use for this evaluation because there are
very few real-world tensors with known ranks. While it is difficult to draw
any concrete conclusions due to the limited number of datasets, we can see
that FRAPPE performs just as well as or slightly better than other existing
rank estimation methods.

###  5.4 CNN-Based Evaluation

We evaluate the performance of our network in selecting the correct number of
components when decomposing a 4D weight tensor for a convolutional neural
network (CNN). We focus specifically on the last layer
(${\mathbb{R}}^{256\times 256\times 3\times 3}$) of an AlexNet [16] model
trained on the CIFAR-10 [15] image dataset.

The goal of this task is to select the correct number of components such that
we compress the weights as much as possible and maintain as high accuracy as
possible. This rank would be the “knee” of the rank-accuracy curve. A plot of
this curve with the values selected by the various rank estimation methods is
shown in Fig. 2.

NORMO predicts a low estimate of 3 for the rank of the weight tensor, and
AutoTen predicts the maximum estimate of 100. Both of these appear to be poor
choices for this task since the accuracy is very low at rank 3 and the gains
are heavily diminished by rank 100. FRAPPE estimates the rank to be 40, which
is not perfect, but appears to be fairly reasonable estimates near the “knee”
of the plot. ARD-Tucker did not terminate within 24 hours and is therefore
excluded from the plot.

Figure 3: Runtime/performance comparison between the different rank estimation
methods on a synthetic dataset. The closer to the lower left, the better.

Model | MAPE | MAE | MSE | Runtime
---|---|---|---|---
AutoTen444Note that we only used the Frobenius-norm-based AutoTen for these synthetic examples due to the presence of negative values in the tensors (causing KL-divergence-based AutoTen to fail). This could result in lower accuracy but also decrease the evaluation times. [20] | $36.99\%$ | 11.55 | 314.80 |  12.77s
ARD [17] | $43.06\%$ | 13.52 | 414.70 |  19.72s
NORMO [8] | $34.36\%$ | 10.44 | 282.38 | 229.34s
FRAPPE | 23.99% |  6.43 |  81.40 |  9.34s

Figure 4: Mean Absolute Percentage Error (MAPE), Mean Absolute Error (MAE),
and Mean Square Error (MSE), and mean runtime of different models on a
synthetic dataset.

###  5.5 Synthetic Dataset Results

Fig. 4 shows that FRAPPE outperforms the baseline methods, with a ${\sim}$10%
improvement in Mean Absolute Percentage Error (MAPE) over NORMO (the best-
performing baseline), ${\sim}$4 rank improvement in Mean Absolute Error (MAE),
and an improvement of 200 in Mean Squared Error (MSE). It is also
${\sim}$24$\times$ faster than NORMO, and ${\sim}$1.6× faster than AutoTen
(the fastest baseline).

###  5.6 Feature Importance Results

In Fig. 5, we plot the importance of the different feature groups (described
in Sec. 4.1 above). We compute the importance of a feature based on the number
of times it is used in the split of a LightGBM tree.

We can see that the most important features are the correlation over different
slices. This makes intuitive sense because higher correlation between slices
along a given mode means that components are more related along that mode,
leading to fewer components being required. This also happens to be the
intuition behind NORMO, one of the baseline methods, which looks instead at
the correlation between different factors of the CPD.

Another interesting set of features are the thresholded non-zero counts. In
Fig. 5, we can see that the importance increases from a threshold of 0.1 to
0.2, but steadily decreases after that until 0.9, which is almost irrelevant.
This is likely because lower relative values (since we normalize the values)
are more likely to be noise artifacts or simply unimportant values in the
decomposition. Lower values are also more likely to be excluded in an $R$-rank
decomposition, simply because the incremental error from excluding that value
is smaller than it would be for a larger-valued element.

$\displaystyle{0}$$\displaystyle{1000}$Mean ImportanceTensor DimensionsNon-
ZerosRanks of SlicesNon-Zeros w/ Threshold of 0.1Non-Zeros w/ Threshold of
0.2Non-Zeros w/ Threshold of 0.3Non-Zeros w/ Threshold of 0.4Non-Zeros w/
Threshold of 0.5Non-Zeros w/ Threshold of 0.6Non-Zeros w/ Threshold of 0.7Non-
Zeros w/ Threshold of 0.8Non-Zeros w/ Threshold of 0.9Correlation over
SlicesFeature Group

Figure 5: Chart of feature group importance (higher is more important), as
output by a LightGBM model trained on a synthetic dataset. The most important
feature group is the correlation over slices of the tensor, followed by the
number of non-zeros thresholded at 0.2. We calculate feature group importance
by how often a feature is used in splits. The feature groups are described in
more detail in Sec. 4.1.

##  6 Additional Related Work

On top of the CPD rank estimation methods described above, there have also
been similar methods developed for different tensor factorizations. Vesselinov
et al. [26] propose NTF$k$, which is an unsupervised method to perform Non-
negative Tensor Factorization (NTF). NTF$k$ uses a custom clustering algorithm
(similar to $k$-means) to cluster the columns of a factor matrix and evaluates
their quality with silhouette score. It then uses a combination of those
scores and the reconstruction error to choose a good target for the number of
components. Nebgen et al. [18] attempt to solve a similar problem for non-
negative matrix factorization (NMF), which decomposes a matrix $\mathbf{V}$
into $\mathbf{V}=\mathbf{W}\mathbf{H}$.

They first repeatedly calculate the NMF on data perturbed by noise from a
random uniform distribution. Then, they use $k$-means to cluster the columns
of $\mathbf{W}$. They then use the Akaike information criterion (AIC) and
silhouette score to calculate the quality of the clusters. Finally, they use a
neural network over a sliding window of the AIC and cluster silhouette values
to determine the ideal number of components.

##  7 Conclusion

In this work, we propose FRAPPE, a novel method for estimating the canonical
rank of a tensor in a self-supervised manner without having to compute the
expensive CPD. The key insight of this approach is that it is much cheaper to
generate synthetic data of known rank than to compute the actual CPD.
Following this idea, we also propose a set of hand-crafted features that can
effectively determine a given tensor’s rank.

We find that naively training a supervised regression model on this data
results in poor generalization to real-world tensors and requires a long
training step. To alleviate this issue, we propose FRAPPE, which first
generates synthetic tensors (with known rank) from an input tensor and trains
a new model on these synthetic tensors, improving the generalization ability
of the model, reducing the amount of training data required, and speeding up
the model training.

We extensively evaluate FRAPPE on three different types of data: synthetic
tensors, real-world tensors, and the weight tensor of a CNN. We find that our
model generally outperforms existing baselines across all three datasets, is
over 24× faster than the best-performing baseline (NORMO), and ${\sim}$1.6×
faster than the fastest baseline (AutoTen).

Finally, we examine the features we used to make these predictions and why we
believe these features are important in estimating a tensor’s rank. We provide
both theoretical reasoning and empirical evidence for why the pair-wise
correlation between slices is a good indicator of the tensor’s rank. We
believe this opens up an interesting line of research about estimating the
number of components of a tensor decomposition without ever having to compute
a single decomposition.

##  8 Ethical Statement

One important ethical consideration regarding this work can arise if it is
used on sensitive information like financial or healthcare data. As with all
machine learning methods, our method is not infallible and may produce
incorrect results. In certain mission-critical systems, that can potentially
lead to severe consequences if care is not taken to validate model results.

Another ethical concern is the possibility of adversarial attacks on the
model. In our work, we do not consider the case of adversarially crafted input
tensors designed to trick the model into performing poorly or predicting an
extremely high rank, which could potentially result in maliciously altered
results.

expansion=sloppy

## References

  * Bader and Kolda [2007] Bader, B., Kolda, T.: Matlab tensor toolbox version 2.2. Albuquerque, NM, USA: Sandia National Laboratories (2007)
  * Bader et al. [2006] Bader, B.W., Harshman, R.A., Kolda, T.G.: Pattern analysis of directed graphs using dedicom: an application to enron email. Tech. rep., Sandia National Laboratories (SNL), Albuquerque, NM, and Livermore, CA … (2006)
  * Bro [1997] Bro, R.: Parafac. tutorial and applications. Chemometrics and Intelligent Laboratory Systems 38(2), 149–171 (1997), ISSN 0169-7439, doi:https://doi.org/10.1016/S0169-7439(97)00032-4
  * Bro [1999] Bro, R.: Exploratory study of sugar production using fluorescence spectroscopy and multi-way analysis. Chemometrics and Intelligent Laboratory Systems 46, 133–147 (1999)
  * Bro and Kiers [2003] Bro, R., Kiers, H.A.: A new efficient method for determining the number of components in parafac models. Journal of chemometrics 17(5), 274–286 (2003)
  * Carroll and Chang [1970] Carroll, J.D., Chang, J.J.: Analysis of individual differences in multidimensional scaling via an n-way generalization of “eckart-young”decomposition. Psychometrika 35(3), 283–319 (1970), doi:10.1007/BF02310791
  * Eagle and [Sandy] Eagle, N., (Sandy) Pentland, A.: Reality mining: Sensing complex social systems. Personal Ubiquitous Comput. 10(4), 255–268 (mar 2006), ISSN 1617-4909, doi:10.1007/s00779-005-0046-3
  * Fernandes et al. [2020] Fernandes, S., Fanaee-T, H., Gama, J.: Normo: A new method for estimating the number of components in cp tensor decomposition. Engineering Applications of Artificial Intelligence 96, 103926 (2020), ISSN 0952-1976, doi:https://doi.org/10.1016/j.engappai.2020.103926
  * Harshman [1970] Harshman, R.: Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis (1970)
  * Harshman et al. [1970] Harshman, R.A., et al.: Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis. UCLA Working Papers in Phonetics (1970)
  * Håstad [1990] Håstad, J.: Tensor rank is np-complete. Journal of Algorithms 11(4), 644–654 (1990)
  * Hillar and Lim [2013] Hillar, C.J., Lim, L.H.: Most tensor problems are np-hard. Journal of the ACM (JACM) 60(6), 45 (2013)
  * Ke et al. [2017] Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., Liu, T.Y.: Lightgbm: A highly efficient gradient boosting decision tree. NeurIPS 30 (2017)
  * Kiers [1998] Kiers, H.: A three-step algorithm for candecomp/parafac analysis of large data sets with multicollinearity. Journal of Chemometrics - J CHEMOMETR 12, 155–171 (05 1998), doi:10.1002/(SICI)1099-128X(199805/06)12:33.3.CO;2-X
  * Krizhevsky et al. [2009] Krizhevsky, A., Hinton, G., et al.: Learning multiple layers of features from tiny images. Citeseer (2009)
  * Krizhevsky et al. [2012] Krizhevsky, A., Sutskever, I., Hinton, G.E.: Imagenet classification with deep convolutional neural networks. In: Pereira, F., Burges, C.J.C., Bottou, L., Weinberger, K.Q. (eds.) Advances in Neural Information Processing Systems, vol. 25, Curran Associates, Inc. (2012)
  * Mørup and Hansen [2009] Mørup, M., Hansen, L.K.: Automatic relevance determination for multi-way models. Journal of Chemometrics 23(7-8), 352–363 (2009)
  * Nebgen et al. [2020] Nebgen, B.T., Vangara, R., Hombrados-Herrera, M.A., Kuksova, S., Alexandrov, B.S.: A neural network for determination of latent dimensionality in non-negative matrix factorization. Machine Learning: Science and Technology 2, 25012 (2020), doi:10.1088/2632-2153/aba372, URL https://doi.org/10.1088/2632-2153/aba372
  * Papalexakis et al. [2013] Papalexakis, E.E., Akoglu, L., Ience, D.: Do more views of a graph help? community detection and clustering in multi-graphs. In: FUSION ’16, pp. 899–905 (2013)
  * Papalexakis, Evangelos E [2016] Papalexakis, Evangelos E: Automatic unsupervised tensor mining with quality assessment. In: SIAM SDM (2016)
  * Pasricha et al. [2018] Pasricha, R., Gujral, E., Papalexakis, E.E.: Identifying and alleviating concept drift in streaming tensor decomposition. CoRR abs/1804.09619 (2018), URL http://arxiv.org/abs/1804.09619
  * Priebe et al. [2005] Priebe, C.E., Conroy, J.M., Marchette, D.J., Park, Y.: Scan statistics on enron graphs. Computational & Mathematical Organization Theory 11(3), 229–247 (Oct 2005), ISSN 1572-9346, doi:10.1007/s10588-005-5378-z
  * Riu and Bro [2003] Riu, J., Bro, R.: Jack-knife technique for outlier detection and estimation of standard errors in parafac models. Chemometrics and Intelligent Laboratory Systems 65(1), 35–49 (2003), ISSN 0169-7439, doi:https://doi.org/10.1016/S0169-7439(02)00090-4
  * Tsitsikas and Papalexakis [2020] Tsitsikas, Y., Papalexakis, E.E.: NSVD : N ormalized S ingular V alue D eviation Reveals Number of Latent Factors in Tensor Decomposition . In: Proceedings of the 2020 SIAM International Conference on Data Mining, pp. 667–675, Society for Industrial and Applied Mathematics (1 2020), doi:10.1137/1.9781611976236.75
  * Tucker [1966] Tucker, L.: Some mathematical notes on three-mode factor analysis. Psychometrika 31(3), 279–311 (1966)
  * Vesselinov et al. [2018] Vesselinov, V.V., Mudunuru, M.K., Karra, S., O’Malley, D., Alexandrov, B.S.: Unsupervised machine learning based on non-negative tensor factorization for analyzing reactive-mixing. CoRR abs/1805.06454 (2018), URL http://arxiv.org/abs/1805.06454
  * Zhao et al. [2014] Zhao, Q., Zhang, L., Cichocki, A.: Bayesian CP factorization of incomplete tensors with automatic rank determination. CoRR abs/1401.6497 (2014), URL http://arxiv.org/abs/1401.6497
  * Zhou et al. [2019] Zhou, M., Liu, Y., Long, Z., Chen, L., Zhu, C.: Tensor rank learning in cp decomposition via convolutional neural network. Signal Processing: Image Communication 73, 12–21 (2019), ISSN 0923-5965, doi:10.1016/j.image.2018.03.017, tensor Image Processing


2021

[1]Thai Binh Nguyen

1]Institute of Innovation, Science and Sustainability, Federation University
Australia, Northways road, Churchill, 3842, Victoria, Australia 2]School of
Information Technology, Faculty of Science, Engineering and Built Environment,
Deakin University, 221 Burwood Highway, Melbourne, 3125, Victoria, Australia

# Depth-based Sampling and Steering Constraints for Memoryless Local Planners

<EMAIL_ADDRESS>Linh Nguyen<EMAIL_ADDRESS>Tanveer Choudhury<EMAIL_ADDRESS>Kathleen Keogh
<EMAIL_ADDRESS>Manzur Murshed<EMAIL_ADDRESS>[  [

###### Abstract

By utilizing only depth information, the paper introduces a novel but
efficient local planning approach that enhances not only computational
efficiency but also planning performances for memoryless local planners. The
sampling is first proposed to be based on the depth data which can identify
and eliminate a specific type of in-collision trajectories in the sampled
motion primitive library. More specifically, all the obscured primitives’
endpoints are found through querying the depth values and excluded from the
sampled set, which can significantly reduce the computational workload
required in collision checking. On the other hand, we furthermore propose a
steering mechanism also based on the depth information to effectively prevent
an autonomous vehicle from getting stuck when facing a large convex obstacle,
providing a higher level of autonomy for a planning system. Our steering
technique is theoretically proved to be complete in scenarios of convex
obstacles. To evaluate effectiveness of the proposed DEpth based both Sampling
and Steering (DESS) methods, we implemented them in the synthetic environments
where a quadrotor was simulated flying through a cluttered region with
multiple size-different obstacles. The obtained results demonstrate that the
proposed approach can considerably decrease computing time in local planners,
where more trajectories can be evaluated while the best path with much lower
cost can be found. More importantly, the success rates calculated by the fact
that the robot successfully navigated to the destinations in different testing
scenarios are always higher than 99.6% on average.

###### keywords:

quadrotor, depth image, steering, sampling, local planning

##  1 Introduction

In recent years, depth cameras have become available and more affordable,
increasing their usage for robotic applications in the industrial and the
research community Popovic2021 . One interesting use case is on quadrotors,
where onboard depth sensors generate data for obstacle avoidance algorithms
(path planners) to keep the drone safe. Many research groups have employed
depth data in their studies on path planning for quadrotors. However, those
depth data remain underutilised in existing works, giving us an opportunity to
further exploit them to increase sampling efficiency and path planning
performance.

The sampling-based technique is ubiquitously used as a trajectory planning
approach for quadrotors in unstructured environments Quan2020 ; however, the
number of sampled collision-free trajectories is limited. A sampled set with a
large percentage of in-collision candidates will require a significant amount
of computational power on collision checking. For instance, Florence $et\ al$.
Florence2020  exploit their piecewise triple-double integrator in time for
state modelling to produce a closed-form future maneuver library. Although
their probabilistic maneuver library can be quickly propagated by uniformly
sampling in the state space, it can also be rich in in-collision candidates.
Since they utilise a k-d tree representation of local structures for the
collision checking process, it requires up-front computation cost to construct
that 3D structure. Also using k-d trees for collision checking, Lopez $et\
al$. Lopez2017  generate a set of motion primitives specified by desired
speeds and heading angles. They uniformly sample possible headings from within
the field of view (FoV) while keeping the forward speed constant. Those motion
primitives can easily satisfy state constraints, such as remaining within the
FoV of the sensor. This approach may work effectively in environments with
sparse obstacles. However, in an obstacle-dense area that would increase the
binary search of a single k-d tree query, an unrefined motion primitive set
can cause an exploded computation load for collision checking. Additionally,
in Ryll2019 , Ryll $et\ al$. uniformly sampled over yaw-heading and final
Euclidean distance from the vehicle to generate a spatially and dynamically
diverse set of motion primitives, which will also suffer from the same
potential issue of heavy workload. While the works mentioned above fuse recent
depth images into a local 3D structure, RAPPIDS - Rectangular Pyramid
Partitioning using Integrated Depth Sensors Bucki2020  is more lightweight
since it exploits only the latest depth image for directly partitioning the
free space using a number of rectangular pyramids. However, RAPPIDS explores
the space by randomly sampling trajectories’ endpoints over a given range in
the sensor’s FoV (sampling space), which is not carefully constrained before
sampling. This approach results in a low-quality sampled set that includes
many in-collision trajectories when the quadrotor flies in an obstacle-filled
area, burdening the collision checking.

Besides computational efficiency, planning performance can also be improved by
directly utilising depth data. All those previously mentioned planners are
classified as mapless and memoryless algorithms. Mapless planners Florence2020
; Lopez2017 ; Ryll2019 ; Florence2018  only employ a local map that is just
enough for local planning, and a memoryless one Bucki2020  plans directly on
sensor data. Distinguishing from those, map-based planning systems Bircher2016
; Oleynikova2016  generally integrate global maps Hornung2013 ; Voxbloxs  $a\
priori$ and a global planning algorithm such as RRT* Karaman2011  to guarantee
planning completeness in tasks such as exploration, Cieslewski2017 ;
Oleynikova2018 ; Lu2022 ; Wagner2022 ; Faria2018  or navigating toward a goal
Cai2021 ; Grando2022 ; Jose2019 . All mapless and memoryless planners cannot
access any map but use sensor data. They are likely to fail to reach the goal
in cluttered environments. For example, Lee2021  introduces an autonomous
system employing RAPPIDS to planning trajectories toward a goal in a cluttered
environment. Their algorithm would stop planning when facing a large obstacle
such as a wall. Therefore, there is still room for enhancement of the
memoryless planner’s performance and ability to autonomously plan and find a
solution.

In those mapless and memoryless algorithms, motion primitives are sampled
without considering available depth information. We can utilise available
depth data to refine the motion primitive set before it comes to collision
checking. This is where our depth-based sampling method contributes.
Specifically, we find and reject obscured primitives’ endpoints or constrain
the sampling space only by querying depth values, automatically excluding many
meaningless (in-collision) candidates in the sampled set. Since this procedure
can be done in constant time, our approach will reduce the computational
workload required in collision checking. Additionally, we will promote the
planning autonomy level by proposing a steering mechanism, increasing the
system’s success rate. Our steering mechanism uses a heuristic algorithm
coupling with a direction utility function to offer existing memoryless
planners a higher autonomy in scenarios of convex obstacles.

To evaluate our proposed approaches, we build a complete planning system that
is capable of autonomously navigating a quadrotor through a cluttered
environment to a given destination. The planning system consists of only a
memoryless local planner, and no explicit global planners are required.
Results show that our depth-based steering mechanism can navigate through
large-obstacle scenarios in which existing depth-based local planners will get
stuck. That is, our approach enables successful autonomous flights toward a
goal without the integration of high-level global planners. Another simulation
benchmark also demonstrates that our depth-based sampling method will reach a
considerably higher computational efficiency than the original uniform
sampling technique used in previous works. To the best of our knowledge, this
is the first work that integrates a depth-based sampling technique and an
autonomous steering mechanism within a memoryless planner framework for
navigating toward a goal in unknown cluttered environments.

Our contributions can be summarized as follows:

1) We propose a depth-based sampling technique for memoryless planners to
improve their computational efficiency using depth sensors.

2) We propose a steering mechanism for an autonomous planning system. The
planning system is capable of autonomously navigating aerial robots through
cluttered environments to a goal, without employing any external global
planners. Our system outperforms the previous memoryless planner Lee2021  in
terms of escaping large obstacles (larger than the FoV).

##  2 Memoryless Local Planning Problem

###  2.1 Memoryless Local Planning Pipeline

Mapless and memoryless planners only use a very limited dynamic memory for
collision checking. Mapless planners Florence2020 ; Lopez2017 ; Ryll2019 ;
Florence2018  fuse the most recent historical sensor measurements into a local
3D data structure (e.g., k-d tree) but not an entire global map for collision
checking. A memoryless planner Bucki2020  exploits only the latest sensor
data. These types of planners lend themselves to local planning purposes.

Figure 1: Depth-based sampling for local path planners.

Fig. 1 describes the workflow of a memoryless planner using depth sensors,
where the positions of collision checking and cost evaluating steps are
interchangeable. Each step would use different computational resources,
depending on the types of planners. For example, memoryless or mapless
planners spend most of their resources on collision checking, while map-based
planners may cost substantial computational power in both evaluating costs
(e.g., volumetric exploration gain Dharmadhikari2020 ) and collision checking
(updating and querying a map). In both cases, motion primitive sampling is
always the first step in the flow. Thus, its efficiency affects the
performance of the entire workflow.

###  2.2 Current Memoryless Local Planning Problem

The uniform sampling technique Florence2020 ; Lopez2017 ; Ryll2019 ;
Florence2018 ; Bucki2020  only offers an unrefined sample set, which might
contain a number of meaningless (in-collision) candidates. Such a low-quality
sampled set (low in collision-free samples) will burden the collision-checking
stage, thus it can be a bottleneck to the system’s performance. Checking for
collision is the most computationally expensive task in the pipeline as we
have to manipulate memory-consuming data structures such as depth images, k-d
trees or Euclidean signed distance field (ESDF) Voxbloxs . If the sampling
space contains too many in-collision candidates, planners would cost much more
computational power and time for the same given output in a planning sequence.
As a result, many planning sequences would have been done for nothing, e.g.,
wasting power and time in evaluating many trajectories that would result in a
collision. Therefore, refining the sampling space can be a promising
technique.

To lower the requirement for computing resources, we need to eliminate a
certain number of in-collision candidates in the sample set by considering
depth data. We can determine a particular type of in-collision trajectory by
comparing the trajectory endpoint’s position with its projection in the depth
matrix. If the endpoint is farther from the robot than its projection’s depth,
it is considered inside obstacles, and its trajectory is considered in-
collision, and vice versa. By analysing the relative position between
trajectories and the occupied space, we categorized them into four types:
(type (0)) collision-free trajectories have their entire path outside the
occupied space, (type (1)) trajectories have their endpoints too close to the
occupied space (the robot will collide when getting to the endpoint), (type
(2)) trajectories have their middle sections inside the occupied space and
their endpoints outside (the robot will collide before getting to the
endpoint) and (type (3)) trajectories have their endpoints inside the occupied
space as in Fig. 2. Although the collision-checking process can reject all
types of in-collision trajectories, the fewer in-collision candidates that
need to be checked, the fewer computing resources will be spent. Therefore, in
this work, we will systematically exclude the type (3), obscured-endpoint
candidates, by constraining the sampling space using raw depth data. Type (1)
and type (2) will be rejected by the collision-checking step afterwards.
Reducing the number of in-collision candidates that require checking will
reduce needed computing resources.

Figure 2: Three types of in-collision polynomial trajectories.

##  3 Depth-based Sampling Constraint for Local Planning

###  3.1 Constrained Sampling Space

For local planning, motion primitives can generally be sampled in input space
or state space. When sampling in input space, a dynamic model will be utilised
to propagate a trajectory, e.g., constant-acceleration point-mass as in
Florence2020 . These propagated trajectories are easily dynamically feasible
as dynamic constraints are already applied to the input space. However,
uniform sampling in input space is not necessarily uniform in configuration
space, where collision checking will be done. On the other hand, the sampling
in state space can guarantee probabilistic optimality, but checking for
dynamic feasibility has to be done afterwards Bucki2020 . Both two mentioned
methods do not consider depth perception before collision checking. The denser
the environment is, the more in-collision trajectories are in the sampled set,
wasting much computational resource in checking many in-collision candidates.
Hence, we propose to utilise available depth data to preprocess sampling space
or sampled set for both mentioned methods in different pipelines shown in Fig.
3a, b.

For sampling in input space as outlined in Fig. 3a, propagated trajectories’
endpoints will then be projected into the latest image plane. We consider a
coordinate system with the axis Ox and Oz parallel to the ground and Oy
pointing downwards such that the quadrotor is at [0,0,0] and the Oxy plane is
parallel to the image plane as sketched in Fig. 2. The axis Oz is
perpendicular to the Oxy and the image plane; thus, it points the depth
direction. Note that we always have available depth matrix outputs directly
from the depth sensor $\mathcal{D}_{w\times h}(\Re)$ where $w$, $h$ are the
width and the height (in pixels) of the input depth image. We also have
$(x,y)\ \forall\ x\in[0,w],y\in[0,h]$ are pixel coordinates of the endpoint’s
projection in the image plane and the associated depth value $d_{x,y}$ by
querying the depth matrix, $d_{o}$ and $d_{p}$ is the distance from the
sampled endpoint to the robot in the Oz axis of the local frame before and
after constraining, respectively. We then compare $d_{x,y}$ with $d_{o}$: if
$d_{o}>d_{x,y}$, the current endpoint is obscured, and the corresponding
trajectory is discarded; otherwise, we accept the sampled endpoint and the
corresponding trajectory.

Figure 3: Sampling space preprocess pipelines for sampling in input space and
state space.

For sampling in state space as represented in Fig. 3b, we propose a technique
for excluding these type (3) trajectories as in Fig. 2. Firstly, we decouple
sampling the endpoints’ position into depth direction and Oxy plane (image
plane). We then adapt a constraint process to the depth sampling space (sample
region). Specifically, we shrink the depth sample region (for each sampled
pixel) to fit into an unoccupied region, pulling all the sampled endpoints
into the free space. We directly utilise available raw data from a depth
sensor to determine the unoccupied region for each sampled pixel $(x,y)$.

Figure 4: Depth-based sampling. Black crosses are initial sampled points; they
can fall into free space or occupied space (gradient dark blue region) to
become obscured sampled points (dash black crosses). After the constraining
process, these obscured ones become depth-based constrained sampled points
(solid blue crosses), and the depth sample region will be shrunk into a
constrained depth sample region (solid blue segment).

In sensor-based collision checking, we only sample in a reliable sensor range.
Therefore, a predefined upper bound $u$ and lower bound $l$ will limit the
sample region to the sensor’s reliable range $[l,u]$. Within this work, we
choose $u=3m$ and $l=1m$, suitable for most popular depth sensors on the
market, such as Intel® RealSense™ D435. If the associated depth value of the
sampled pixel $d_{x,y}$ falls inside the reliable range $[l,u]$
($d_{x,y}\in[l,u]$), we scale down the initial depth sample region $[l,u]$ to
the constrained depth sample region $[l,d_{x,y}]$ for each sampled pixel
$(x,y)$. This manipulation will convert obscured sampled points into depth-
based constrained sampled points, lying in the free space, as illustrated in
Fig. 4.

If the associated depth value exceeds the upper bound ($d_{x,y}>u$), in-range
sampled points will always lie in the free space. If the associated depth is
smaller than the lower bound ($d_{x,y}<l$), the generated trajectory will
collide with obstacles. However, in this case, we keep the in-range sampled
point and let the collision-checking process rejects this in-collision
candidate. We do not need to constrain the depth sample region in these two
cases: $d_{p}=d_{o}$. That means we are not directly rejecting any candidate
in sampling but constraining the sampling space, maintaining its continuity
and probabilistic optimality. The sampling space for trajectories’ endpoints
before and after constraining is shown:

$\begin{split}\Omega_{o}=&\\{(x,y,d_{o})|x\in[0,w],y\in[0,h],d_{o}\in[l,u]\\}\mathchar
24635\relax\;\\\ \Omega_{p}=&\\{(x,y,d_{p})|x\in[0,w],y\in[0,h],\\\
&d_{p}\in[l,u]\setminus[d_{x,y},u]\\}\mathchar 24635\relax\;\\\ \end{split}$
(1)

where $\Omega_{o}$, $\Omega_{p}$, $(x,y,d_{o})$, $(x,y,d_{p})$ are original
space, constrained space, and their corresponding sample points, respectively.

$\begin{split}d_{p}=\begin{cases}(d_{o}-l)\frac{d_{x,y}-l}{u-l}+l,&\text{if
}d_{x,y}\in[l,u]\\\ d_{o},&\text{otherwise}\end{cases}\end{split}$ (2)

Details of the depth-based sampling method are described in Algorithm 1.

Algorithm 1 Depth-based sampling

Input: Latest depth image $\mathcal{D}$, user-defined range: lower bound $l$,
upper bound $u$

Output: Sampled point $P$ as the endpoint’s position of the next trajectory
for evaluating

1:function depth_based_sample()

2: Randomly sample a pixel $(x,y)$ from image frame;

3: Randomly sample $d_{o}$ $\in$ [$l$,$u$];

4: Query $d_{x,y}$ value in $\mathcal{D}$;

5: if $d_{x,y}\in[l,u]$ then

6: $d_{p}=(d_{o}-l)\frac{d_{x,y}-l}{u-l}+l$;

7: else

8: $d_{p}=d_{o}$;

9: end if

10: Deproject current candidate ($x$,$y$,$d_{p}$) to a real-world endpoint
$P$;

11:end function

###  3.2 Planning Algorithm

We propose the path planning algorithm described in Fig. 5 in this work. The
local planner aims to find the lowest-cost collision-free trajectory over a
given user-defined computation time. The planner will run many sequences until
the time limit is reached. For each sequence, we generate a minimum-average-
jerk polynomial trajectory Mueller2015  using our depth-based sampling method
described in 3.1. The collision checking function follows the algorithm
RAPPIDS in Bucki2020 . This is the most computationally demanding step and is
only called when the sampled trajectory’s cost is lower than the current best
cost. This way, the best trajectory’s cost converges to the lowest over time.

Figure 5: Planning pipeline.

We propose a direction function $\mathcal{J}(p)$ for evaluating cost
trajectories as follows:

$\begin{split}O&=p(0)\mathchar 24635\relax\;P=p(T)\\\
\overrightarrow{d}&=\frac{\overrightarrow{OG}}{\|\overrightarrow{OG}\|}\\\
\mathcal{J}(p)&=\frac{\overrightarrow{d}.\overrightarrow{PO}}{\|\overrightarrow{PO}\|}\end{split}$
(3)

where $p$ is the polynomial trajectory to be evaluated; $O$, $G$ and $P$ are
the world-frame positions of the vehicle, the goal and the trajectory’s
endpoint, respectively; $T$ is the execution time of the trajectory.
$\overrightarrow{d}$ is the vector of exploration direction. This vector-based
function corrects the velocity vector to the desired direction. The works
Lee2021  and Bucki2020  employ a function prioritising the average velocity to
the goal, allowing the lateral motion to be free. Thus it causes a more
diverged and fluctuated path, even in completely free space. In our proposed
system, $\overrightarrow{d}$ will be specified as the unit vector from the
vehicle’s current position to the goal. The cost is simply the dot product of
the exploration vector $\overrightarrow{d}$ and the negative unit vector of
the motion primitive. That is, the direction of a better trajectory will align
closer to the goal direction. This feature also paves a platform for the
autonomous steering algorithm, which we will discuss later in section 4.

###  3.3 Algorithm Evaluation

We compared the planner’s performance in terms of the cost of the best
candidate found. We applied multiple Monte Carlo simulations111An
implementation of the algorithm and the simulation source code will be
available online at https://github.com/thethaibinh/dess as described in
Bucki2020  on the two approaches: uniform sampling and depth-based sampling.
The simulations were based on 1000 uniformly synthetic obstacle scenarios and
robot states. We ran the two mentioned planners 1000 times for each of the 15
allocated computation time values, ranging from 0 to 20 milliseconds. We
repeated the same simulation setup on two hardware platforms: an i7-1185G7
laptop and a Jetson Nano embedded computer. We divided computation time to
have higher resolution at shorter computation time regions. The results
depicted in Fig. 6 show that the depth-based sampling planner finds a
considerably lower cost trajectory than the uniform sampling one, especially
in short computation time regions or a less powerful platform. Since we
constrained depth sample regions, it is a higher chance for us to sample a
collision-free candidate. Therefore, there is a higher ratio of collision-free
candidates in the sampled set, the computation power spent on collision
checking decreases, resulting in more sampled trajectories being evaluated as
in Fig. 7. Therefore, the planner then again has spare time to sample more
meaningful trajectories in the space, which include more potential collision-
free trajectories as in Fig. 8.

Figure 6: Best cost values for proposed depth-based technique and uniform
sampling used in RAPPIDS Bucki2020 .

Figure 7: Comparing number of evaluated trajectories for proposed depth-based
technique and uniform sampling used in RAPPIDS Bucki2020 .

Figure 8: Comparing number of generated collision-free trajectories for
proposed depth-based technique and uniform sampling used in RAPPIDS Bucki2020
.

###  3.4 Performance Analysis

We evaluate the sampling methods’ efficiency by doing a time complexity
analysis. Given a depth image containing $n$ pixels and a computation time,
accessing a value from a depth image can be done in constant time. The
percentage of obscured-endpoint candidates in the uniformly sampled set and
the depth-based sampled set is $RD_{aic}$ and $DB_{aic}$ respectively, and the
average number of elements $N_{avg}$ in sampled sets in both cases are equal.
Checking an obscured endpoint trajectory requires the planner to inflate a new
pyramid, which also cost an $O(n)$ operation. The depth-based technique
excluded all obscured-endpoint candidates, thus $DB_{aic}$ is zero. Since
random generators are uniform distributions, $RD_{aic}$ roughly equals the
percentage of occupied volume out of the total sample range volume (in the
FoV) as illustrated in Fig. 4. In other words, we can call it the level of
density $L_{dens}$, thus $RD_{aic}$ = $L_{dens}$. For both uniform and depth-
based pipelines, we assume that the ratio $K_{c}$ of the number of candidates
that pass the better-cost checking is equal, e.g., 100 trajectories pass the
cost checking to come for collision checking out of 1000 sampled trajectories,
following the planning pipeline sketched in Fig 5. Therefore, the depth-based
pipeline will save us an operation of $O(L_{dens}N_{avg}K_{c}n)$ for each
sequence. In free space, uniform and depth-based sampling would perform
equally because $L_{dens}=0$. But the depth-based method would spare us a
substantial amount of computational power in densely occupied structures where
$L_{dens}$ is considerable.

##  4 Depth-based Steering for Continuous Planning

###  4.1 Trajectory Generation for Stuck Motions

We propose a steering command for the local planner to autonomously generate
feasible trajectories toward a goal. The procedure of generating a normal
trajectory follows the pipeline described in subsection 3.2. We request
steering commands when the local planner cannot find any feasible conventional
trajectory for one second (vehicle to rest at the endpoint of the last
executed trajectory). The steering command is a motion primitive with the
position as the latest executed trajectory’s endpoint and an additional yaw
orientation. New depth data will come due to changes in the FoV as the
quadrotor executes steering commands. The planner will consistently send
steering commands until it finds at least a feasible trajectory from this
latest new depth frame. This trajectory generation scheme is summarised in
Algorithm 2.

Algorithm 2 Steering command generation pipeline

1:while not at goal do

2: generate normal trajectory;

3: if cannot generate any feasible normal trajectory for more than one second
then

4: generate and send steering commands to the quadrotor;

5: else

6: send the best-direction trajectory found to the quadrotor;

7: end if

8:end while

We propose a depth-based mechanism for generating these steering commands.
When the vehicle is stuck, we scan the latest depth image to locate the
nearest point in the space, which means the closest obstacle. We divide space
in half by a vertical plane that passes through the vertical centreline of the
image. On which half-space the nearest point lies, the quadrotor steers toward
the other half as visualised in Fig. 9. The desired yaw angle will be
accumulated by steering value over time on the controller side and calculated
as follows,

$\gamma=\begin{cases}1,&\text{if }x_{c}<(w/2)\\\
-1,&\text{otherwise}\end{cases}$ (4) $\psi_{r}=\psi_{f}+\gamma*K_{p}*\Delta t$
(5)

where:

  * •

$\gamma$ is the sign value representing steering direction, $\gamma=1$ means
steering right, $\gamma=-1$ means steering left.

  * •

$x_{c}$ is the pixel location of the closest point along the x-axis of the
image frame.

  * •

$w$ is the width of the image frame in pixel unit.

  * •

$\psi_{r}$ is the desired reference yaw angle for yaw control when executing a
steering command.

  * •

$\psi_{f}$ is the feedback of the current yaw angle.

  * •

$K_{p}$ is the accumulation gain for the yaw angle.

  * •

$\Delta t$ is the cycle time of a control loop.

Figure 9: Steering value calculation.

The mechanism is simple yet effective with the type of geometry-based planners
since it is hungry for new input. When a sufficiently clear space view comes,
the planner will generate feasible trajectories again, pulling the quadrotor
out of infeasible configurations, especially when facing a massive obstacle
that is even bigger than the FoV. We can prove by geometry that it is feasible
to navigate the robot over large convex obstacles.

###  4.2 Yaw Control in Obstacle-dense Areas

When it comes to travelling long distances over many obstacles, we need
another mechanism to integrate the best direction trajectory and depth-based
steering to control. Instead of applying a fixed yaw control that faces the
camera toward the goal Lee2021 , we adopt a yaw control scheme that faces the
camera to the local goal (the latest feasible trajectory’s endpoint) by
setting the desired control yaw angle $\psi_{d}$ equal to the bearing of the
vehicle-goal vector $\psi_{b}$. We only do that when the position error
vector’s length $\|\overrightarrow{r_{err}}\|$ is greater than one meter.
Otherwise, we set $\psi_{d}$ equal to either $\psi_{r}$ when the robot needs
to do a steering command or the latest executed control yaw value
$\psi_{d}^{-}$ when it does not (Algorithm 3). In this way, the system
automatically prioritizes the exploration direction that navigates toward the
free space, steering away from an obstacle before bumping into it. As a
result, the quadrotor would autonomously find a clearer path, escaping
cluttered areas faster and safer, avoiding more potential ”stuck” situations
(need to steer).

Algorithm 3 Yaw control

1:if $\|\overrightarrow{r_{err}}\|>1m$ then

2: $\psi_{d}=\psi_{b}$

3:else

4: if is steering then

5: $\psi_{d}=\psi_{r}$;

6: else

7: $\psi_{d}=\psi_{d}^{-}$;

8: end if

9:end if

10:$\psi_{d}^{-}=\psi_{d}$;

###  4.3 Completeness

Figure 10: Steering through obstacles scenarios, including the wall and the
spheroid. When the vehicle is far away from any obstacle, it will perfectly
track the vehicle-goal line to hit the steering point $p_{steer}$ on the
boundary of obstacles where we need to steer for a more open field of view.
This phase is illustrated by a solid green arrow connecting the robot’s
current position to the steering point $p_{steer}$. After that steering, the
next consecutive best-cost trajectories (solid green arrows) would form a line
parallel to the wall’s peripheral until they end up at a visible point
$p_{visible}$ where the goal is visible to the vehicle. As formulated in 3,
best-cost candidates are collision-free and form the smallest angle with
vehicle-goal vector (dash green angle). The entire maneuver is visualized in
sub-figure a). The drone needs time to change its heading to the desired
angle. Thus, the practical trajectory will be like the one connected by best-
direction trajectories, which is illustrated by curved solid arrows in sub-
figure b) and c). The remaining black dashed line arrows are rejected
candidates in the sampled set. The solid red arrows denote the steering
direction, following the steering mechanism proposed in 4.1. From the analysis
of sub-figure a), we can also infer that when the vehicle is close to the
obstacle, the best path will align and snap to the obstacle’s peripheral as
shown in sub-figure b) and c)).

Here we follow up on our previous discussions and formally prove the
completeness of our proposed algorithm. In this work, we only consider convex
obstacles.

We assume there is no uncertainty in the control process, and the planning is
locally optimal. That is, the robot can flawlessly track the given path, and
the local planner always finds the best-cost trajectory, which is closest to
the goal direction and collision-free.

On the other hand, we define the vehicle-goal line as the unique line that
passes through the robot’s current position and goal point. The robot may need
to execute a steering command (we now call the act of executing a steering
command to steer for brevity) only when obstacles lie in the vehicle-goal
line.

######  Lemma 1.

When the vehicle faces a convex obstacle, it will reach a visible point
$p_{visible}$, where the goal is visible in the FoV, after only the first
steer.

Proof: All convex obstacles can be represented by the wall and spheroid shapes
as illustrated in Fig. 10. We consider scenarios when a given goal lying on
the other side of the obstacle with the robot. When the vehicle is away from a
wall, it will perfectly track the vehicle-goal line to hit the steering point
$p_{steer}$ on the wall’s boundary where we need to steer for a more open
field of view. After that steering, the first feasible trajectory’s direction
will be parallel to the current side of the wall. And the subsequent
consecutive best-cost trajectories would form a line parallel to the wall’s
peripheral until they end up at a visible point $p_{visible}$ where the goal
is visible to the vehicle (as shown in Fig. 10a)). Indeed, these trajectories
have the best direction costs as they are the minus dot product of two unit
vectors of the vehicle-goal line and trajectory’s direction as formulated in
3. Spheroids and other convex obstacle scenarios can be proved by following
the same geometrical approach as visualised in Fig. 10b,c. Thus, when the
vehicle faces a convex obstacle, it will reach a visible point $p_{visible}$
after only the first steer.

######  Theorem 1.

The steering algorithm is complete in all scenarios of convex obstacles.

Proof: Let us prove the theorem by contradiction. That is, we assume that the
steering algorithm is incomplete. Consider situations when the robot faces a
wall and a large spheroid on the way to the goal. If the steering algorithm
does not find the path, then necessarily it will follow an infinite path or
have infinite times of steer to get to the goal. Suppose the vehicle has to
travel an infinite path to get to the goal. The path will either diverge from
the goal or be a local cyclic course in free space. This conflicts with our
previous assumptions. On the other hand, if the vehicle has to steer forever,
the goal would never be visible to the robot. However, Lemma 1 shows that when
facing a large wall or spheroid, the robot will take a maximum of one steer to
have the goal visible in the FoV, thus the goal cannot be obscured forever.
Therefore, the proposed algorithm always be able to find a feasible path to
the goal.

##  5 Implementations

###  5.1 Setup

Figure 11: Planning system overview based on quadrotor simulator from
DodgeDrone Challenge 2022. The mission is to navigate a simulated quadrotor to
fly through obstacle-dense environments. The obstacles are different in size
and randomly located in the space.

We compare planning performance between two policies: the fixed yawing
obstacle avoidance method in Lee2021  and our proposed planning system, which
we now call DESS \- DEpth-based Sampling and Steering for brevity. The fixed
yawing policy utilises a utility function of average velocity to the goal, and
a collision-checking algorithm with uniform sampling Bucki2020 , while DESS
manipulates a direction function, collision checking with depth-based sampling
and an integrated steering mechanism for stuck motions. We evaluate their
success rate in navigating the quadrotor to reach a goal.

We conduct experiments for that flight task on a simulation system222Our
simulation system source code is available online at
https://github.com/thethaibinh/agile_flight of a quadrotor flying through an
obstacle-dense area from DodgeDrone Challenge 2022 DodgeDrone . The challenge
is to navigate a quadrotor through the obstacle-dense area to a 65m far ahead
goal. For the convenience of bulk autotests, we decrease the goal distance to
15m. The obstacle-dense area is a 3D bounding box with dimensions of [0, 15,
-5, 5, 0, 10] (meters) in a coordinate system with the axis Ox and Oy parallel
to the ground and Oz pointing upwards, in which the quadrotor is at [0, 0, 0]
and the goal is at [17, 0, 5]. Obstacles are solid spheres of different
diameters ranging from 0.1 to 4.0 meters, stationary, and randomly generated
inside the bounding box. We generate three obstacle-dense scenarios based on
predefined density levels: easy, medium and hard so that there are 29, 51 and
67 obstacles in the bounding box for each scenario, respectively. Fig. 11b)
and c) capture a typical scenario with additional visualisation of point cloud
ground truth. All obstacles of the easy scenario are already included in the
medium scenario, and the hard scenario has all of the medium obstacles in its
configuration. A trial will be labelled a success when the robot reaches the
goal’s position before the timeout and does not collide with any obstacle. We
integrated the planner with classical nested feedback loops Michael2010  as a
low-level controller who received position guidance commands in a receding
horizon manner. We limit the maximum desired velocity in the position feedback
controller to 1m/s, mitigating control uncertainty because we only focus on
planning evaluation. Whenever a new depth image arrives, if the planner finds
a collision-free trajectory, the position controller drops the previous
trajectory and starts tracking the newfound trajectory. Fig. 11a) describes an
overview of the system. We will run the entire system 1000 times for each
policy (fixed yawing, DESS) in each scenario (easy, medium, hard) on an
i7-1185G7 laptop; thus, we will have 6000 flight simulations.

Table 1: Success rate of finishing trials by policies | Fixed yawing Lee2021  | DESS (proposed)
---|---|---
Easy | 97.6% | 99.7%
Medium | 69.3% | 99.9%
Hard | 55.19% | 99.6%

Figure 12: Traveling distances to the goal, summarizing only successful
trials. For the fixed yawing method: 976 successes in the easy, 693 successes
in the medium and 552 successes in the hard out of 1000 runs on each. For our
DESS method: 997 successes in the easy, 999 successes in the medium and 996
successes in the hard out of 1000 runs on each.

Figure 13: Traveling time to the goal, summarizing only successful trials. For
the fixed yawing method: 976 successes in the easy, 693 successes in the
medium and 552 successes in the hard out of 1000 runs on each. For our DESS
method: 997 successes in the easy, 999 successes in the medium and 996
successes in the hard out of 1000 runs on each.

###  5.2 Results

Whenever the vehicle meets the situation of not finding any feasible
trajectory (e.g., stuck in front of a huge obstacle or a wall), the DESS
planner generates a steering command. Following a steering command, the robot
will steer consistently until its FoV is clear enough for a normal collision-
free trajectory to be generated. Thus, it continues on its mission to the
goal. Being equipped with our DESS method, the drone always finds a feasible
path with success rates for three scenarios are 99.7%, 99.9%, and 99.6%,
respectively, while simulations using the fixed yawing method will fail more
often when scenarios get denser. A few failed corner cases of DESS appear only
due to uncertainties of the position tracking controller. The results are
summarized in Table 1, reflecting theoretical analysis of DESS’s completeness
on convex obstacles.

On the other hand, we will evaluate the finishing time and travel distance for
successful trials of the two methods. Fig. 12 shows that the average travelled
distances and also their standard deviations of our DESS are approximately
equal to those of fixed yawing’s successful trials. This result demonstrates
that DESS guarantees robustness while maintaining paths converged. However, it
is noted that our method produced much more successful trials than the fixed
yawing technique. Fig. 13 illustrates that the mean values of DESS’s finishing
times are slightly lower than those of the fixed yawing’s successful cases. It
also can be noted that standard deviations of DESS are 1.7, 2.0 and 2.7
seconds for the easy, medium and hard, respectively, significantly smaller
than those of the fixed yawing, which are 3.1, 4.9 and 5.7 seconds,
correspondingly. This comparison shows that there are many trials where the
fixed yawing takes a long time to complete, while DESS consistently reaches
the goal in more steady amounts of time. That is, DESS is more robust and
performs stably in environments of different dense levels. The depth-based
sampling gives a head-start in terms of efficiency, allowing DESS to find the
best-direction trajectory by evaluating the direction function. Therefore,
DESS can find a better direction trajectory, avoiding fluctuated paths caused
by lateral motions from the average velocity function. This direction function
also paves a concrete platform for the depth-based steering algorithm to
manipulate, enabling a higher level of autonomy and a better planning
performance for our proposed planner compared to the fixed yawing approach.

##  6 Conclusion

In this paper, we have presented the depth-based sampling and steering
technique that improves memoryless planners’ efficiency. Our proposed sampling
method exploited available depth data to constrain sampling space, excluding
obscured endpoint samples, which are candidates that are not meaningful (in-
collision) for planning. Thus, the enhanced planner generates a better-cost
trajectory over the same amount of allocated computation time compared with
the same planner that uses the uniform sampling technique. Additionally, we
built a complete planning system that adopted the proposed depth-based
steering algorithm based on that advantage. The steering technique equips the
local planner with the capability of autonomously escaping potential stuck
motions such as large convex obstacles, enabling a higher success rate for the
planning system in obstacle-dense scenarios. The simulation results have
demonstrated that the proposed approach outperforms the state-of-the-art
technique. Since the adopted planner only guarantees instantaneous best
direction at every planning sequence, extensions to this work may include a
mechanism to translate it into a local optimal such as the shortest distance
travelled or fastest time to escape one obstacle when integrating with other
techniques. Acknowledgments This research was supported by the Henry Sutton
scholarship from Federation University Australia.

## Declarations

  * •

Funding: Federation University Australia funded the research via the Henry
Sutton scholarship - Application ID: 3056759.

  * •

Conflict of interest: The authors have no relevant financial or non-financial
interests to disclose

  * •

Ethics approval: No ethical approval is required by this research.

  * •

Consent to participate: Not applicable

  * •

Consent for publication: This paper does not require any consent for
publication.

  * •

Data and/or Code availability: All generated data and implementation codes for
simulations will be available and maintained in the author’s online repository
https://github.com/thethaibinh. If the reader has further needs and questions,
do not hesitate to contact the corresponding author.

## References

  *   * (1) Popovic, M., Thomas, F., Papatheodorou, S., Funk, N., Vidal-Calleja, T., Leutenegger, S.: Volumetric occupancy mapping with probabilistic depth completion for robotic navigation. IEEE Robotics and Automation Letters 6, 5072–5079 (2021). https://doi.org/10.1109/LRA.2021.3070308
  * (2) Quan, L., Han, L., Zhou, B., Shen, S., Gao, F.: Survey of uav motion planning. IET Cyber-Systems and Robotics 2, 14–21 (2020). https://doi.org/10.1049/IET-CSR.2020.0004
  * (3) Florence, P., Carter, J., Tedrake, R.: Integrated perception and control at high speed: Evaluating collision avoidance maneuvers without maps. Springer Proceedings in Advanced Robotics 13, 304–319 (2020). https://doi.org/10.1007/978-3-030-43089-4_20
  * (4) Lopez, B.T., How, J.P.: Aggressive 3-d collision avoidance for high-speed navigation. Proceedings - IEEE International Conference on Robotics and Automation, 5759–5765 (2017). https://doi.org/10.1109/ICRA.2017.7989677
  * (5) Ryll, M., Ware, J., Carter, J., Roy, N.: Efficient trajectory planning for high speed flight in unknown environments. Proceedings - IEEE International Conference on Robotics and Automation 2019-May, 732–738 (2019). https://doi.org/%****␣sn-article.bbl␣Line␣125␣****10.1109/ICRA.2019.8793930
  * (6) Bucki, N., Lee, J., Mueller, M.W.: Rectangular pyramid partitioning using integrated depth sensors (rappids): A fast planner for multicopter navigation. IEEE Robotics and Automation Letters 5, 4626–4633 (2020). https://doi.org/10.1109/LRA.2020.3003277
  * (7) Florence, P.R., Carter, J., Ware, J., Tedrake, R.: Nanomap: Fast, uncertainty-aware proximity queries with lazy search over local 3d data. Proceedings - IEEE International Conference on Robotics and Automation, 7631–7638 (2018). https://doi.org/10.48550/arxiv.1802.09076
  * (8) Bircher, A., Kamel, M., Alexis, K., Oleynikova, H., Siegwart, R.: Receding horizon next-best-view planner for 3d exploration. Proceedings - IEEE International Conference on Robotics and Automation 2016-June, 1462–1468 (2016). https://doi.org/%****␣sn-article.bbl␣Line␣175␣****10.1109/ICRA.2016.7487281
  * (9) Oleynikova, H., Burri, M., Taylor, Z., Nieto, J., Siegwart, R., Galceran, E.: Continuous-time trajectory optimization for online uav replanning. IEEE International Conference on Intelligent Robots and Systems 2016-November, 5332–5339 (2016). https://doi.org/10.1109/IROS.2016.7759784
  * (10) Hornung, A., Wurm, K.M., Bennewitz, M., Stachniss, C., Burgard, W.: Octomap: An efficient probabilistic 3d mapping framework based on octrees. Autonomous Robots 34, 189–206 (2013). https://doi.org/10.1007/S10514-012-9321-0/FIGURES/18
  * (11) Oleynikova, H., Taylor, Z., Fehr, M., Nieto, J., Siegwart, R.: Voxblox: Incremental 3d euclidean signed distance fields for on-board mav planning. IEEE International Conference on Intelligent Robots and Systems 2017-September, 1366–1373 (2016). https://doi.org/10.1109/IROS.2017.8202315
  * (12) Karaman, S., Frazzoli, E.: Sampling-based algorithms for optimal motion planning (2011). https://doi.org/10.48550/arxiv.1105.1186
  * (13) Cieslewski, T., Kaufmann, E., Scaramuzza, D.: Rapid exploration with multi-rotors: A frontier selection method for high speed flight. IEEE International Conference on Intelligent Robots and Systems 2017-September, 2135–2142 (2017). https://doi.org/10.1109/IROS.2017.8206030
  * (14) Oleynikova, H., Lanegger, C., Taylor, Z., Pantic, M., Millane, A., Siegwart, R., Nieto, J.: An open-source system for vision-based micro-aerial vehicle mapping, planning, and flight in cluttered environments. Journal of Field Robotics 37, 642–666 (2018). https://doi.org/10.1002/rob.21950
  * (15) Lu, B.X., Tseng, K.S.: 3d map exploration using topological fourier sparse set. Journal of Intelligent & Robotic Systems 104, 1–22 (2022). https://doi.org/10.1007/S10846-021-01565-1
  * (16) Wagner, A., Peterson, J., Donnelly, J., Chourey, S., Kochersberger, K.: Online aerial 2.5d terrain mapping and active aerial vehicle exploration for ground robot navigation. Journal of Intelligent & Robotic Systems 106, 1–18 (2022). https://doi.org/10.1007/S10846-022-01751-9
  * (17) Faria, M., Maza, I., Viguria, A.: Applying frontier cells based exploration and lazy theta* path planning over single grid-based world representation for autonomous inspection of large 3d structures with an uas. Journal of Intelligent & Robotic Systems 93, 113–133 (2018). https://doi.org/%****␣sn-article.bbl␣Line␣325␣****10.1007/S10846-018-0798-4
  * (18) Cai, K., Wang, C., Song, S., Chen, H., Meng, M.Q.H.: Risk-aware path planning under uncertainty in dynamic environments. Journal of Intelligent & Robotic Systems 101, 1–15 (2021). https://doi.org/10.1007/S10846-021-01323-3
  * (19) Grando, R.B., de Jesus, J.C., Kich, V.A., Kolling, A.H., Drews-Jr, P.L.J.: Double critic deep reinforcement learning for mapless 3d navigation of unmanned aerial vehicles. Journal of Intelligent & Robotic Systems 104, 1–14 (2022). https://doi.org/10.1007/S10846-021-01568-Y
  * (20) Sanchez-Lopez, J.L., Wang, M., Olivares-Mendez, M.A., Molina, M., Voos, H.: A real-time 3d path planning solution for collision-free navigation of multirotor aerial robots in dynamic environments. Journal of Intelligent & Robotic Systems 93, 33–53 (2019). https://doi.org/10.1007/s10846-018-0809-5
  * (21) Lee, J., Wu, X., Lee, S.J., Mueller, M.W.: Autonomous flight through cluttered outdoor environments using a memoryless planner. 2021 International Conference on Unmanned Aircraft Systems, ICUAS 2021, 1131–1138 (2021). https://doi.org/10.1109/ICUAS51884.2021.9476874
  * (22) Dharmadhikari, M., Dang, T., Solanka, L., Loje, J., Nguyen, H., Khedekar, N., Alexis, K.: Motion primitives-based path planning for fast and agile exploration using aerial robots. 2020 IEEE International Conference on Robotics and Automation (ICRA), 179–185 (2020). https://doi.org/10.1109/ICRA40945.2020.9196964
  * (23) Mueller, M.W., Hehn, M., D’Andrea, R.: A computationally efficient motion primitive for quadrocopter trajectory generation. IEEE Transactions on Robotics 31, 1294–1310 (2015). https://doi.org/10.1109/TRO.2015.2479878
  * (24) Song, Y., Kaufmann, E., Bauersfeld, L., Loquercio, A., Scaramuzza, D.: ICRA 2022 DodgeDrone Challenge: Vision-based Agile Drone Flight. https://uzh-rpg.github.io/icra2022-dodgedrone/
  * (25) Michael, N., Mellinger, D., Lindsey, Q., Kumar, V.: The grasp multiple micro-uav testbed. IEEE Robotics and Automation Magazine 17, 56–65 (2010). https://doi.org/10.1109/MRA.2010.937855