# A new form of the Hahn - Banach Theorem

Sokol Bush Kaliaj  Mathematics Department, Science Natural Faculty, University
of Elbasan, Elbasan, Albania<EMAIL_ADDRESS>
###### Abstract.

In this paper, we present a new form of the Hahn - Banach Theorem in terms of
the sub-additive convex functions.

###### Key words and phrases:

Hahn-Banach theorem, sub-additive convex functionals.

###### 2010 Mathematics Subject Classification:

Primary 46A22; Secondary 46A32

##  1\. Introduction and Preliminaries

The Hahn-Banach theorem is one of the central results in functional analysis.
This theorem originated from Hahn [5] and Banach [1] is of basic importance in
the analysis of problems concerning the existence of continuous linear
functionals. Its principal formulations are as a dominated extension theorem
(analytic form) and as a separation theorem (geometric form), c.f. [12,
Theorem II, 3.2] and [12, Theorem II, 3.1]. It is known that [12, Theorem II,
3.2] and [12, Theorem II, 3.1] imply each other.

Throughout this paper, $L$ denotes a vector space over $\Phi$. The scalar
field $\Phi$ is the real field $\mathbb{R}$ or the complex field $\mathbb{C}$.
We use the terminologies in the book [12]. A function $p:L\to\mathbb{R}$ is
said to be a sub-linear function if for every $x,y\in L$ and for every
$\lambda\geq 0$, we have

  * •

$p(x+y)\leq p(x)+p(y)$,

  * •

$p(\lambda x)=\lambda p(x)$.

The function $p$ is said to be a semi-norm if $p$ is a sub-linear function
such that for every $x\in L$ and for every $\lambda\in\Phi$, we have

$p(\lambda x)=|\lambda|p(x).$

A function $\varphi:L\to\mathbb{R}$ is said to be a convex function if for
every $x,y\in L$ and for every $\lambda\in\mathbb{R}$, we have

$0<\lambda<1\Rightarrow\varphi(\lambda
x+(1-\lambda)y)\leq\lambda\varphi(x)+(1-\lambda)\varphi(y).$

We say that $\varphi$ is a sub-additive function if

$\varphi(x+y)\leq\varphi(x)+\varphi(y),\text{ for all }x,y\in L,$

and $\varphi(\theta)=0$, where $\theta$ is the zero vector. The function
$\varphi$ is called a sub-additive convex function if $\varphi$ is a sub-
additive and convex function.

Let $L$ be a topological vector space with $0$-neighborhood base
$\mathfrak{B}$ and let $\varphi:L\to\mathbb{R}$ be a convex function. The
function $\varphi$ is said to be upper semi-continuous at $x_{0}\in L$ if
given $\varepsilon>0$ there exists a neighborhood
$U_{\varepsilon}\in\mathfrak{B}$ such that

$x\in x_{0}+U_{\varepsilon}\Rightarrow\varphi(x)<\varphi(x_{0})+\varepsilon;$

$\varphi$ is said to be upper semi-continuous on a subset $G\subset L$ if for
each $x\in G$, $\varphi$ is upper semi-continuous at $x$.

The following version of the analytic Hahn-Banach theorem first appears in
Banach [1] (see also [3, Theorem II.3.10]).

######  Theorem 1.1.

Let $\mathfrak{Y}$ be a vector subspace of the real vector space
$\mathfrak{X}$ and let $p:\mathfrak{X}\to\mathbb{R}$ be a sub-linear function.
Let $f:\mathfrak{Y}\to\mathbb{R}$ be a linear function with

$f(x)\leq p(x),\quad x\in\mathfrak{Y}.$

Then there exists a linear function $F:\mathfrak{X}\to\mathbb{R}$ for which

$F(x)=f(x),\quad x\in\mathfrak{Y};\qquad F(x)\leq p(x),\quad
x\in\mathfrak{X}.$

In this paper, we present a new form of the analytic Hahn-Banach theorem in
terms of the sub-additive convex functions, see Theorem 2.3. Since each sub-
linear function is sub-additive convex, Theorem 2.3 implies Theorem 1.1. It is
also shown that Theorem 2.3 implies [12, Theorem II, 3.2], see Corollary 2.4.

In the papers [8], [9] and [4] are presented generalizations of the Hahn -
Banach theorem in terms of the convex functionals. There are also
generalizations of the Hahn - Banach theorem in the papers [6], [2], [10],
[11] and [13]. Kakutani [7] gave a proof of the analytic Hahn-Banach theorem
by using the Markov-Kakutani fixed-point theorem.

##  2\. The Main Result

The main result is Theorem 2.3. Let us start with two auxiliary lemmas.

######  Lemma 2.1.

Let $L$ be a topological vector space over $\mathbb{R}$, let $\mathfrak{B}$ be
the family of all circled $0$-neighborhood and let $\varphi:L\to\mathbb{R}$ be
a convex function. If $G$ is a non-empty, convex, open subset of $L$ on which
$\varphi$ is bounded above, then $\varphi$ is upper semi-continuous on $G$.

###### Proof.

Let $x_{0}\in G$ be and arbitrary vector and let $0<\varepsilon<1$. Then by
hypothesis there exists $M>0$ such that

$\varphi(x)\leq M,\text{ for all }x\in G$

and there exists $B\in\mathfrak{B}$ such that $x_{0}+B\subset G$. Thus,
$\varphi(x)\leq M$ for all $x\in x_{0}+B$. If we choose

$\eta=\frac{\varepsilon}{3(1+M)(1+|\varphi(x_{0})|)}$

then $\eta B\subset B$, since $B$ is circled. For every $u\in B$ the equality

$x_{0}+\eta u=\eta(x_{0}+u)+(1-\eta)x_{0}$

yields

$\begin{split}\varphi(x_{0}+\eta
u)&\leq\eta\varphi(x_{0}+u)+(1-\eta)\varphi(x_{0})\\\
&=\varphi(x_{0})+\eta\varphi(x_{0}+u)-\eta\varphi(x_{0})\leq\varphi(x_{0})+\eta
M-\eta\varphi(x_{0})\\\ &\leq\varphi(x_{0})+\frac{\varepsilon
M}{3(1+M)(1+|\varphi(x_{0})|)}+\frac{\varepsilon(-\varphi(x_{0}))}{3(1+M)(1+|\varphi(x_{0})|)}\\\
&\leq\varphi(x_{0})+\frac{\varepsilon}{3}+\frac{\varepsilon}{3}<\varphi(x_{0})+\varepsilon.\end{split}$

Thus, for every $x\in x_{0}+\eta B$, we have

$\varphi(x)<\varphi(x_{0})+\varepsilon.$

This means that $\varphi$ is upper semi-continuous at $x_{0}\in G$, and since
$x_{0}$ was arbitrary it follows that $\varphi$ is upper semi-continuous on
$G$ and the proof is finished. ∎

######  Lemma 2.2.

Let $L$ be a topological vector space over $\mathbb{R}$, let $\mathfrak{B}$ be
the family of all circled $0$-neighborhood and let $\varphi:L\to\mathbb{R}$ be
a convex function. Then the following statements are equivalent:

  * (i)

$\varphi$ is continuous on $L$.

  * (ii)

$\varphi$ is upper semi-continuous on $L$.

  * (iii)

There exists a non-empty, convex, open subset of $L$ on which $\varphi$ is
bounded above.

###### Proof.

$(i)\Leftrightarrow(ii)$ Clearly, if $\varphi$ is continuous on $L$, then
$\varphi$ is upper semi-continuous on $L$. Conversely, assume that $\varphi$
is upper semi-continuous at $x_{0}\in L$. Then, given $\varepsilon>0$ there
exists $U_{\varepsilon}\in\mathfrak{B}$ such that

(2.1) $u\in
U_{\varepsilon}=-U_{\varepsilon}\Rightarrow\max\\{\varphi(x_{0}+u),\varphi(x_{0}-u)\\}<\varphi(x_{0})+\varepsilon$

and since

$2\varphi(x_{0})\leq\varphi(x_{0}+u)+\varphi(x_{0}-u)$

it follows that

(2.2) $u\in U_{\varepsilon}\Rightarrow
2\varphi(x_{0})\leq\varphi(x_{0}+u)+\varphi(x_{0}-u)<2\varphi(x_{0})+2\varepsilon.$

The last result together with (2.1) yields

$u\in
U_{\varepsilon}\Rightarrow\varphi(x_{0})-\varepsilon<\varphi(x_{0}+u)<\varphi(x_{0})+\varepsilon,$

since $\varphi(x_{0}+u)\leq\varphi(x_{0})-\varepsilon$ together with
$\varphi(x_{0}-u)\leq\varphi(x_{0})+\varepsilon$ yields

$\varphi(x_{0}+u)+\varphi(x_{0}-u)<2\varphi(x_{0})$

contradicting (2.2). This means that $\varphi$ is continuous at $x_{0}$.

$(ii)\Leftrightarrow(iii)$ Assume that $\varphi$ is upper semi-continuous on
$L$ and let $G=\\{x\in L:\varphi(x)<1\\}$. It is easy to see that $G$ is non-
empty convex subset of $L$. Let $x_{0}\in G$ be an arbitrary vector and let
$\varepsilon=\frac{1-\varphi(x_{0})}{2}$. Then there exists
$U_{\varepsilon}\in\mathfrak{B}$ such that

$\begin{split}x\in
x_{0}+U_{\varepsilon}\Rightarrow\varphi(x)<\varphi(x_{0})+\varepsilon<1\Rightarrow
x\in G.\end{split}$

Thus, $x_{0}+U_{\varepsilon}\subset G$, and since $x_{0}$ was arbitrary it
follows that $G$ is an open subset of $L$.

Conversely, assume that $G$ is a non-empty, convex, open subset of $L$ on
which $\varphi$ is bounded above and let $x_{0}\in L$ be an arbitrary vector.
Then, there exists $M>0$ such that $\varphi(x)\leq M$ for all $x\in G$. If
$x_{0}\in G$, then by Lemma 2.1 $\varphi$ is upper semi-continuous at $x_{0}$.

It remains to consider $x_{0}\not\in G$. Fix a vector $y_{0}\in G$ and a real
positive number $\rho>1$ and define $w=y_{0}+\rho(x_{0}-y_{0})$. The function
$h:L\to L$ defined by

$h(y)=\frac{\rho-1}{\rho}y+\frac{1}{\rho}w,\text{ for all }y\in L$

is a homeomorphism, c.f. [12, 1.1, p.13]. The function $h$ transforms $y_{0}$
into $x_{0}$ and $G$ into an open and convex set $h(G)$ containing $x_{0}$.
For every $x\in h(G)$, we have

$x=h(h^{-1}(x))=\frac{\rho-1}{\rho}h^{-1}(x)+\frac{1}{\rho}w$

and

$\begin{split}\varphi(x)&=\varphi\left(\frac{\rho-1}{\rho}h^{-1}(x)+\frac{1}{\rho}w\right)\leq\frac{\rho-1}{\rho}\varphi(h^{-1}(x))+\frac{1}{\rho}\varphi(w)\\\
&\leq\frac{\rho-1}{\rho}M+\frac{1}{\rho}\varphi(w)=M^{\prime}<+\infty.\end{split}$

Thus, $\varphi$ is bounded over $h(G)$. Therefore, by Lemma 2.2, $\varphi$ is
upper semi-continuous at $x_{0}\in h(G)$. Since $x_{0}$ was arbitrary
$\varphi$ is upper semi-continuous at every $x\in G$ and at every $x\not\in
G$. Thus, $\varphi$ is upper semi-continuous on $L$, and this ends the proof.
∎

We are now ready to present the main result.

######  Theorem 2.3.

Let $M$ be a subspace of a vector space $L$ over $\mathbb{R}$ and let
$\varphi:L\to\mathbb{R}$ be a sub-additive convex function. If
$g:M\to\mathbb{R}$ is a linear function such that $g$ ”dominated” by
$\varphi$, i.e.,

$(\forall x\in M)[g(x)\leq\varphi(x)],$

then there exists a linear function $f:L\to\mathbb{R}$ such that

  * (i)

$f$ extends $g$ to $L$, i.e.,

$f|_{M}=g,$

  * (ii)

$f$ ”dominated” by $\varphi$, i.e.,

$(\forall x\in L)[f(x)\leq\varphi(x)].$

###### Proof.

Define

$U_{n}=\left\\{x\in L:\varphi(x)<\frac{1}{n}\right\\},n\in\mathbb{N}$

and

$\mathfrak{B}=\left\\{V_{n}\in
2^{L}:V_{n}=U_{n}\cap(-U_{n}),n\in\mathbb{N}\right\\},$

where

$-U_{n}=\left\\{-x\in L:\varphi(x)<\frac{1}{n}\right\\}=\left\\{x\in
L:\varphi(-x)<\frac{1}{n}\right\\}.$

The family $\mathfrak{B}$ is a filter base in $L$, since

$V_{n}\cap V_{m}\supset V_{m},\quad m>n$

and $\theta\in V_{n}$ for all $n\in\mathbb{N}$.

Note that $(-1)V_{n}=V_{n}$ and $0V_{n}=\\{\theta\\}\subset V_{n}$; given
$0<\lambda<1$ or $-1<\lambda<0$, we have also

$\lambda V_{n}\subset V_{n},$

since

$\begin{split}\varphi(x)&=\varphi\left(\frac{\lambda}{\lambda}x\right)=\varphi\left(\lambda\frac{1}{\lambda}x+(1-\lambda)\theta\right)\leq\lambda\varphi\left(\frac{1}{\lambda}x\right).\end{split}$

Thus, $V_{n}$ is circled. Let us show next that $V_{n}$ is radial. Since
$V_{n}$ is circled, it is enough to show that given $x_{0}\in L$ there exists
$\mu\in\mathbb{R}$ such that $x_{0}\in\mu V_{n}$. If $x_{0}\in V_{n}$ then
$\mu=1$. Otherwise $x_{0}\not\in U_{n}$ or $x_{0}\not\in-U_{n}$;
$\varphi(x_{0})\geq\frac{1}{n}$ or $\varphi(-x_{0})\geq\frac{1}{n}$. There
exists $\mu>1$ such that $\frac{1}{\mu}\varphi(x_{0})<\frac{1}{n}$ and
$\frac{1}{\mu}\varphi(-x_{0})<\frac{1}{n}$. Since

$\begin{split}\varphi\left(\frac{1}{\mu}x_{0}\right)=&\varphi\left(\frac{1}{\mu}x_{0}+\left(1-\frac{1}{\mu}\right)\theta\right)\leq\frac{1}{\mu}\varphi(x_{0})<\frac{1}{n},\\\
\varphi\left(-\frac{1}{\mu}x_{0}\right)=&\varphi\left(\frac{1}{\mu}(-x_{0})+\left(1-\frac{1}{\mu}\right)\theta\right)\leq\frac{1}{\mu}\varphi(-x_{0})<\frac{1}{n},\end{split}$

it follows that $x_{0}\in\mu U_{n}\cap\mu(-U_{n})=\mu V_{n}$. Thus,
$\mathfrak{B}$ is a filter base in $L$ such that for each $n\in\mathbb{N}$, we
have

  * (a)

$V_{2n}+V_{2n}\subset V_{n}$, since $\varphi$ is sub-additive function,

  * (b)

$V_{n}$ is radial and circled.

Therefore, by [12, I, 2.1, p.15] the family $\mathfrak{B}$ is a
$0$-neighborhood base for a unique topology $\mathfrak{T}$ under which $L$ is
a topological vector space.

Given any vector $x_{0}\in U_{1}$, there exists $n\in\mathbb{N}$ such that
$\frac{1}{n}<1-\varphi(x_{0})$. Then, $x_{0}+V_{n}\subset U_{1}$. This means
that $U_{1}$ is an open set, and since $U_{1}$ is also convex, by Lemma 2.2 it
follows that $\varphi$ is continuous on $L$.

We now consider $L\times\mathbb{R}$ as the product of two topological vector
spaces. Note that if $T:L\times\mathbb{R}\to\mathbb{R}$ is a linear function,
then

$T(x,t)=h(x)+\alpha t,\text{ for all }(x,t)\in L\times\mathbb{R},$

where $h(x)=T(x,0)$ for all $x\in L$ and, $\alpha=T(0,1)$. Hence,

$H=\\{(x,t)\in M\times\mathbb{R}:g(x)-t=1\\}$

is a linear manifold in the subspace $M\times\mathbb{R}$. Since $(x,t)\to
S(x,t)=g(x)-t$ is a non-zero linear functional, by [12, I, 4.1] $H$ is a
hyperplane in $M\times\mathbb{R}$ and a linear manifold in
$L\times\mathbb{R}$.

Since $\varphi$ is a continuous function on $L$ it follows that

$G=\\{(x,t)\in L\times\mathbb{R}:\varphi(x)-t<1\\}$

is an open set, and since $g(x)\leq\varphi(x)$ for $x\in M$ it follows that
$G\cap H=\emptyset$; $G$ is also a non-empty convex set. Therefore, by [12,
Theorem II, 3.1] there exists a closed hyperplane $H_{1}$ in
$L\times\mathbb{R}$ such that $H_{1}\supset H$ and $H_{1}\cap G=\emptyset$.
Since $H_{1}$ is a hyperplane, $H_{1}=(x_{0},t_{0})+H_{0}$ where $H_{0}$ is a
maximal subspace of $L\times\mathbb{R}$. If we choose $(x_{0},t_{0})\in
H\subset H_{1}\cap M\times\mathbb{R}$, then

$H_{1}\cap(M\times\mathbb{R})=(x_{0},t_{0})+(H_{0}\cap(M\times\mathbb{R})),$

and since

$H_{1}\cap(M\times\mathbb{R})\neq(M\times\mathbb{R}),\quad((0,0)\not\in
H_{1}\text{ and }(0,0)\in(M\times\mathbb{R}))$

it follows that $H_{1}\cap(M\times\mathbb{R})$ is a hyperplane in
$M\times\mathbb{R}$. Further, $H_{1}\cap(M\times\mathbb{R})\supset H$ implies

(2.3)
$H_{1}\cap(M\times\mathbb{R})=H=(x_{0},t_{0})+(H_{0}\cap(M\times\mathbb{R}));$

($M\times\mathbb{R}$ is the algebraic sum of subspaces
$H_{0}\cap(M\times\mathbb{R})$ and
$\\{\lambda(x_{0},t_{0}):\lambda\in\mathbb{R}\\}$). By [12, I,4.1], there
exits a non-zero linear functional $F:L\times\mathbb{R}\to\mathbb{R}$ such
that $H_{1}=\\{(x,t)\in L\times\mathbb{R}:F(x,t)=1\\}$. The least result
together with (2.3) yields

$F(x,0)=S(x,0)=g(x),\text{ for all }x\in M;$

that is, the linear function $f(\cdot)=F(\cdot,0)$ is an extension of $g$ to
$L$.

It remains to prove that $f$ ”dominated” by $\varphi$. Note that

$F(x,t)=f(x)-t,\text{ for all }(x,t)\in L\times\mathbb{R},$

since $(0,-1)\in H\subset H_{1}$. By $H_{1}\cap G=\emptyset$ it follows that

(2.4) $f(x)=1+t\leq\varphi(x),\text{ for all }(x,t)\in H_{1}.$

While $(x,t)\not\in H_{1}$ implies that there exists a real number $r\neq 1$
such that

$f(x)-t=r.$

Then, $(x,t+r-1)\in H_{1}$ and therefore (2.4) yields

$f(x)=1+(t+r-1)\leq\varphi(x).$

Thus, $f$ ”dominated” by $\varphi$ and this ends the proof. ∎

######  Corollary 2.4.

Let $M$ be a subspace of a vector space $L$ over $\mathbb{C}$ and let
$p:L\to[0,+\infty)$ be a semi-norm. If $g:M\to\mathbb{C}$ is a linear
functional such that $g$ ”dominated” by $p$, i.e.,

$(\forall x\in M)[|g(x)|\leq p(x)],$

then there exists a linear functional $f:L\to\mathbb{C}$ such that

  * •

$f$ extends $g$ to $L$, i.e.,

$f|_{M}=g,$

  * •

$f$ ”dominated” by $p$, i.e.,

$(\forall x\in L)[|f(x)|\leq p(x)].$

###### Proof.

There exists a linear functional $g_{0}:M\to\mathbb{R}$ such that

$g(x)=g_{0}(x)-ig_{0}(ix),\text{ for all }x\in M.$

We have also that $p$ is a convex functional such that $p(\lambda x)=\lambda
p(x)$ for every $\lambda\geq 0$ and $x\in L$. Thus, regarding $L$ as a real
linear space, and applying Theorem 2.3, a real linear function
$f_{0}:L\to\mathbb{R}$ is obtained for which $f_{0}|_{M}=g_{0}$ and
$f_{0}(x)\leq p(x)$ for all $x\in L$. Let the function $f:L\to\mathbb{C}$ be
defined by the equation

$f(x)=f_{0}(x)-if_{0}(ix).$

Clearly, $f$ is a linear functional and

$f(x)=f_{0}(x)-if_{0}(ix)=g_{0}(x)-ig_{0}(ix)=g(x),\text{ for all }x\in M.$

Thus $f$ is an extension of $g$. Finally, let $f(x)=|f(x)|e^{i\theta(x)}$,
where $\theta(x)$ is an argument of the complex number $f(x)$; then

$|f(x)|=f_{0}\left(\frac{x}{e^{i\theta(x)}}\right)\leq
p\left(\frac{x}{e^{i\theta(x)}}\right)=p(x),$

which proves that $|f(\cdot)|\leq p(\cdot)$ and this ends the proof. ∎

Since [12, Theorem II, 3.1] implies Theorem 2.3 and Theorem 2.3 implies [12,
Theorem II, 3.2], we conclude that Theorem 2.3, [12, Theorem II, 3.2] and [12,
Theorem II, 3.1] imply each other.

## References

  * [1] S. Banach, Sur les fonctionnelles linéaires II, Studia Math. (1929), 223 - 229.
  * [2] S. Çakan and Y. Yilmaz, A Generalization of the Hahn-Banach Theorem in Seminormed Quasilinear Spaces, Journal of Mathematics and Applications 42 (2019), 79-94.
  * [3] N. Dunford and J. T. Schwartz, Linear operators. Part I: General theory, Interscience New York and London (1958).
  * [4] L. Forg c, On the Hahn-Banach theorem, Mathematica Slovaca, vol. 26/1 (1976), 39-45.
  * [5] H. Hahn, Über lineare Gleichungssysteme in linearen Räumen, J. Reine Angew. Math. 157 (1926), 214-229.
  * [6] N. Hirano, H. Komiya, W. Takahashi, A Generalization of the Hahn-Banach Theorem, JMAA 88 (1982), 333-340.
  * [7] S. Kakutani, Two fixed-point theorems concerning bicompact convex sets, Proc. Imp. Acad. Tokyo 14 (1938), 242-245.
  * [8] H. Nakano, On an extension theorem, Proc. Japan. Acad., 35/3 (1959), 127.
  * [9] H. Nakano, On the Hahn Banach theorem, Bull. Acad. Polon. Sci., Ser. Sci. Math. Astronom. Phys., 19/8 (1971), 743-745.
  * [10] J. Plewnia, A generalization of the Hahn-Banach theorem, Annanes Polonici Math. LVIII.1 (1993).
  * [11] G. Rano Hahn-Banach extension theorem in quasi-normed linear spaces, Advances in Fuzzy Mathematics, 12/4 (2017), 825-833.
  * [12] H. H. Schaeffer, Topological Vector Spaces, (Third Printing, Graduate Texts in Mathematics), Springer-Verlag New York (1971).
  * [13] S. Simons, A new version of the Hahn-Banach theorem, Arch. Math. 80 (2003), 630-648.


# Obstacle Identification and Ellipsoidal Decomposition for Fast Motion
Planning in Unknown Dynamic Environments

Mehmetcan Kaymaz and Nazım Kemal Ure Mehmetcan Kaymaz is with the Faculty of
Aeronautics and Astronautics, Istanbul Technical University, Istanbul, Turkey
<EMAIL_ADDRESS>Kemal Ure is with the Faculty of Computer and
Informatics Engineering, Istanbul Technical University, Istanbul, Turkey
<EMAIL_ADDRESS>
###### Abstract

Collision avoidance in the presence of dynamic obstacles in unknown
environments is one of the most critical challenges for unmanned systems. In
this paper, we present a method that identifies obstacles in terms of
ellipsoids to estimate linear and angular obstacle velocities. Our proposed
method is based on the idea of any object can be approximately expressed by
ellipsoids. To achieve this, we propose a method based on variational Bayesian
estimation of Gaussian mixture model, the Kyachiyan algorithm, and a
refinement algorithm. Our proposed method does not require knowledge of the
number of clusters and can operate in real-time, unlike existing optimization-
based methods. In addition, we define an ellipsoid-based feature vector to
match obstacles given two timely close point frames. Our method can be applied
to any environment with static and dynamic obstacles, including ones with
rotating obstacles. We compare our algorithm with other clustering methods and
show that when coupled with a trajectory planner, the overall system can
efficiently traverse unknown environments in the presence of dynamic
obstacles.

##  I INTRODUCTION

Obstacle avoidance in unknown environments is a common problem for unmanned
systems. In general, the identification of obstacles and trajectory planning
are considered separate problems, and the computational complexity of both
algorithms is desired to be as low as possible[1] to compensate for the
uncertainty of the environment. In scenarios that involve dynamic obstacles,
the identification of obstacles has a critical effect on motion estimation, as
the estimation performance strongly depends on the accurate identification of
such obstacles.

Clustering is a popular method for separating obstacles based on point cloud
data provided by perception systems. The number of estimated clusters can make
a big difference in the overall performance. There are several approaches to
compute the number of clusters; such as assuming the constant number of
clusters[2][3], calculating the number of clusters from the number of
points[4], or estimating the number of clusters based on point data[5][6].
Most of the previous work assumes that each dynamic obstacle can be expressed
as a single circle[7] or an ellipsoid[8][9] and moves without rotation[1]. It
has been demonstrated that the ellipsoid shape is a better approximation to
approximate obstacles in such environments[10]. Although these assumptions are
true for a lot of real-world applications, these methods face severe
limitations as the environment structure gets more complex.

In this work, we model the obstacle identification problem as the computation
of minimum volume ellipsoids[11], which is defined as finding $n$ ellipsoids
that cover $m$ given points where the total volume of ellipsoids is minimized.
The optimal solution can be obtained via mixed-integer semidefinite
programming(MISDP). However, MISDP formulation requires providing the number
of ellipsoids beforehand, and the computational complexity of solving MISDP is
not feasible for online planning. We propose a method that provides an
approximate solution for the minimum volume ellipsoids problem without knowing
the number of ellipsoids a priori. In addition, the computational requirement
of our proposed algorithm is much lesser than MISDP, which enables running
trajectory planning algorithms concurrently with obstacle identification.

In summary, the main contributions of the paper are as follows:

  * •

To solve the minimum volume ellipsoids problem without the knowledge of the
number of ellipsoids, we combine variational Bayesian estimation of Gaussian
mixture model(VIGMM)[5] with the Khachiyan algorithm[12] and add custom
refinements to improve the performance.

  * •

To match obstacles given two timely close point frames, we define a feature
vector based on ellipsoid parameters and estimate motion by using matched
ellipsoids.

  * •

Overall, our proposed work can rapidly identify dynamic obstacles in cluttered
environments in real time. The simulation results show that the proposed
identification method can be coupled with a trajectory planning method to
traverse a variety of dynamic environments without too much computational
overhead.

##  II RELATED WORK

Clustering is a common method to identify obstacles in collision avoidance
approaches, especially when the obstacles are dynamic. K-means[13][14] is a
famous method in the literature for clustering due to its low computation
requirement. In [2], the authors used the K-means algorithm with a fixed
number of clusters to classify point cloud data from monocular SLAM in a
static and simple environment. Hierarchical clustering was used with the
single-linkage clustering method in [15], where point cloud data was estimated
from the stereo vision system. A similar idea is applied in [16] with lidar to
generate point cloud data. DBSCAN[6] based on Euclidean distance is a common
method for obstacle avoidance especially when the point data is noisy. In [9],
point cloud data obtained from the event camera was clustered using the DBSCAN
algorithm for dynamic obstacle avoidance. The same clustering method is also
used in [1] for both dynamic and static obstacles where the point cloud is
estimated from the depth image. The problem with the methods mentioned above
is that they fail when the clusters are Gaussian data points in which the
point distribution is ellipsoidal, so Euclidian distance is not a sufficient
metric to separate clusters.

The MISDP is proposed to cluster Gaussian data points in [11] but the
algorithm’s computational complexity is high and it requires the number of
clusters as input. A method based on C-means was proposed in [17], which
successfully clusters Gaussian data points but the computational cost is high,
similar to the MISDP. The possibilistic C-means algorithm[18][19], which is
the primitive version of [17] is more computationally efficient but it fails
when the clusters overlap. Also, C-means methods require the number of
clusters. Another method to cluster Gaussian data points is Gaussian mixture
models(GMM) which are based on fitting multivariate Gaussian distributions to
data. In [3], GMM was concatenated with the artificial potential field method
for static obstacle avoidance where the number of clusters is selected
constant. In [20], online mapping based on GMM with a constant number of Gauss
components is proposed for static collision avoidance. [4] used GMM to
generate the online map but the number of Gauss components is estimated via
$N/R$ where $N$ is the number of points and $R$ is a constant parameter.

Matching the obstacles/clusters from sequential point frame measurements is
important for estimating the motion of dynamic objects. In [21] and [22], the
centers of the clusters were used to match obstacles whereas the ones with a
minimum distance between center points are assumed as the same obstacle.
However, this method may fail when the centers of two different obstacles are
close to each other. To improve the matching robustness, the feature vector is
designed in [1] which is more robust than other alternatives but requires RGB
vision inputs.

##  III PROBLEM DEFINITION

In this section, we present a minimum volume ellipsoids problem to cover a
given point of data. The problem is defined in two dimensional (2D) plane.

We use three different ellipsoid definitions. The first one is called the
general form, which norm based definition:

$\mathcal{E}=\\{x\textbf{ }|\textbf{ }||Ax+b||_{2}\leq 1\textbf{ ,
}x\in\mathbb{R}^{2}\\},$ (1)

where $A\in\mathbb{S}^{2}_{++}$ and $b\in\mathbb{R}^{2}$. The second
definition is called the quadratic form, which is the open version of the
general form:

$\mathcal{E}=\\{x\textbf{ }|\textbf{ }x^{T}A^{q}x+2(b^{q})^{T}x+c^{q}\leq
0\textbf{ , }x\in\mathbb{R}^{2}\\},$ (2)

where $A^{q}\in\mathbb{S}^{2}_{++}$, $b^{q}\in\mathbb{R}^{2}$ and
$c^{q}\in\mathbb{R}$. The last definition is called the standard form, which
is a modified version of the quadratic form:

$\mathcal{E}=\\{x\textbf{ }|\textbf{
}(x-x_{c})^{T}R(\theta)^{T}H(r_{1},r_{2})R(\theta)(x-x_{c})\leq 1\textbf{ ,
}x\in\mathbb{R}^{2}\\},$ (3)

where $x_{c}\in\mathbb{R}^{2}$ is the center of the ellipsoid,
$R(\theta)\in\mathbb{R}^{2\times 2}$ is the rotation matrix,
$\theta\in[0,\pi]$ is the rotation angle,
$H(r_{1},r_{2})\in\mathbb{S}^{2}_{++}$ is the ellipsoidal zone,
$r_{1}\in\mathbb{R}^{+}$ and $r_{2}\in\mathbb{R}^{+}$ are the minor and major
axes of the ellipsoid.

Let the total number of ellipsoids be given as $n_{ell}\in\mathbb{N}^{+}$,
hence the space covered by all ellipsoids is defined as
$\mathcal{E}_{all}=\mathcal{E}_{1}\cup...\cup\mathcal{E}_{n_{ell}}\subset\mathbb{R}^{2}$.
Let the shape of the obstacle be represented by $n_{obs}\in\mathbb{N}^{+}$
number of points and let $K=\\{1,2,...,n_{obs}\\}$. The position of $k$th
point is denoted as $\nu_{k}\in\mathbb{R}^{2}$. The set of all points in given
shape is donated as
$O_{obs}=\\{\nu_{1},\nu_{2},...,\nu_{K}\\}\subset\mathbb{R}^{2}$. Given points
are covered if and only if all points in the shape are the elements of at
least one ellipsoid. In other words, the set $O_{obs}$ must be a subset of
$\mathcal{E}_{all}$.

##  IV METHOD

Our proposed framework is composed of two sub-modules: first, we present a
method for solving the minimum volume ellipsoids problem by using variational
Bayesian estimation of the Gaussian mixture model, the Khachiyan algorithm,
and our refinement algorithm. Then, we show how to match obstacles from given
two frames measured in a short time interval and how to estimate the motion of
the obstacles via ellipsoids.

###  IV-A Minimum Volume Ellipsoids

Our approximate method starts with the Gaussian mixture model(GMM), which is
used as the clustering method. The GMM is a probabilistic model with multiple
Gaussian components. Given a set of points $\nu$, the probability density
function of a GMM can be expressed as:

$p(\nu|\pi,\mu,\Sigma)=\sum_{j=1}^{J}\pi_{j}\mathcal{N}(\nu;\mu_{j},\Sigma_{j}),$
(4)

where $J$ is the number of clusters, $\mathcal{N}(\nu;\mu_{j},\Sigma_{j})$ is
multivariate Gaussian distribution, $\mu_{j}$ is mean vector, $\Sigma_{j}$ is
covariance matrix and $\pi_{j}$ is the weight coefficient of the $j$th
component, satisfying $\sum\pi_{j}=1$.

The latent variables $\lambda=\\{\lambda_{1},...,\lambda_{J}\\}$ can be
estimated efficiently via the Expectation-Maximization(EM)[23] method, but the
number of Gaussian components must be determined manually. The Maximum
Likelihood(ML)[24] method can be applied, but computing the global optimum is
difficult. The Monte-Carlo Markov Chain(MCMC)[25] method addresses both of the
issues in EM and ML, but the computational requirement is too big. On the
other hand, Variational Inference(VI)[26] estimates the parameters of GMM
through optimization, which is usually faster than MCMC. In VI, the maximum
number of Gaussian components must be determined, the method eliminates the
unnecessary ones.

The problem of estimating GMM parameters can be expressed as finding a
posterior distribution $p(\lambda|\nu)$. In VI, a variational distribution
$q(\lambda)$ is used to approximate $p(\lambda|\nu)$ by minimizing Kullback-
Leible(KL) divergence between them. It is not possible to calculate
$KL(q(\lambda)||p(\lambda|\nu))$ since $p(\lambda|x)$ is unknown. However, the
following equation can be written using the Bayes rule:

$KL(q(\lambda)||p(\lambda|\nu))=\mathcal{L}(q)-\int_{\lambda}q(\lambda)lnp(\nu)d\lambda,$
(5)

where $\mathcal{L}(q)$ is Evidence Lower Bound(ELBO). $p(\nu)$ is a definite
value, so minimizing KL diverge is equivalent to maximizing $L(q)$.

After clustering, the minimum volume ellipsoids problem is converted to
minimum volume ellipsoid for $J$ clusters which is much easier to calculate.
The minimum volume ellipsoid problem is defined as finding the minimum volume
ellipsoid that covers given $n_{obs}$ points which can be written in semi-
definite programming form[27]:

[2] A,b(detA)^1n ——Aν+b——_2≤1,

where $A$ and $b$ represent ellipsoid in the general form. The Khachiyan
algorithm[12] is used for solving (IV-A), which provides a solution with a
predefined error. The Kyachiyan algorithm is faster than semi-definite
programming solvers even with a very small predefined error[28]. The issue in
the variational Bayesian estimation of the Gaussian mixture model is observed
in Fig. 1. Although the VI method eliminates unnecessary Gaussian components,
the result depends on the number of maximum components. The refinement
algorithm uses the Khachiyan algorithm outputs, which are $n_{ell}$ ellipsoids
as input and combines a subset of them based on their volumes.

Figure 1: Variational Bayesian estimation of Gaussian mixture model clustering
based on the number of maximum components. Each color shows a different
cluster and n is the maximum number of clusters.

Assume we have $n_{ell}$ ellipsoids and $J=\\{1,2,...,n_{ell}\\}$. We define
ellipsoids volume ratio $r_{ell_{i,j}}$ as:

$r_{ell_{i,j}}=\frac{Vol(\mathcal{E}_{i})+Vol(\mathcal{E}_{j})}{Vol(OBB(\mathcal{E}_{i},\mathcal{E}_{j}))},$
(6)

where $i,j\in J$, $Vol$ is the volume and
$OBB(\mathcal{E}_{i},\mathcal{E}_{j})$ is oriented bounding box of
$\mathcal{E}_{i}$ and $\mathcal{E}_{j}$. Calculation of $r_{ell_{i,j}}$ is
computationally efficient as the volume of $\mathcal{E}_{i}$ and
$\mathcal{E}_{j}$ are already known, and calculating OBB is a simple
operation[29]. In order to combine two ellipsoids, we use the minimum volume
ellipsoid covering the union of ellipsoids[27], which is a convex optimization
problem that can be solved efficiently for two ellipsoids in two or three
dimensions. The semi-definite programming for given two ellipsoids in the
quadratic form is defined as: [3] A^2,~b,τ_1,τ_2(detA)^1n τ_i ≥0
[A2-τiAqi~b-τibqi0 (~b-τibqi)T-1-τicqi~bT0 ~b-A2]⪯0 i=1,2,

where $\tilde{b}=Ab$. The Mosek software[30] is used for solving the semi-
definite problem defined in (6). The higher values of $r_{ell_{i,j}}$ means
combining $i$th and $j$th ellipsoids is a acceptable approximation. We define
a threshold $r_{ell_{threshold}}$ to decide which ellipsoids should be
combined. The refinement algorithm is presented in Alg. 1

Algorithm 1 The Refinement Algorithm

0: $\\{\mathcal{E}_{1},...,\mathcal{E}_{n_{ell_{v}}}\\}$ list of ellipsoids,
$r_{ell_{threshold}}$ threshold to combine ellipsoids.

0: $\\{\mathcal{E}_{1},...,\mathcal{E}_{n_{ell_{r}}}\\}$ list of ellipsoids

1: for $i=1$ to $n_{ell_{v}}$ do

2: for $j=i+1$ to $n_{ell_{v}}$ do

3: $r_{ell_{i,j}}=f(\mathcal{E}_{i},\mathcal{E}_{j})$ (6)

4: if $r_{ell_{i,j}}\geq r_{ell_{threshold}}$ then

5: Delete $i$th and $j$th ellipsoids

6: $A_{i,j},b_{i,j}=f(\mathcal{E}_{i},\mathcal{E}_{j})$(6)

7: Add combined ellipsoid $A_{i,j},b_{i,j}$

8: end if

9: end for

10: end for

11: return $\\{\mathcal{E}_{1},...,\mathcal{E}_{n_{ell_{r}}}\\}$

The refinement algorithm can be used until the number of ellipsoids does not
change. Yet in our simulations, we did not encounter any scenario that
requires using the refinement algorithm more than once.

Figure 2: Overall algorithm pipeline. Each color shows the different clusters
obtained from VIGMM.

###  IV-B Obstacle Matching and Motion Estimation

We define the obstacle matching method and motion estimation method similar to
[1]. The feature vector is defined as ellipsoid parameters in the standard
form.

$fte(\mathcal{E})=[x_{c},r_{1},r_{2},\theta]$ (7)

The idea is that if the change in future vectors given two timely close frames
is small, they are considered the same object. The Euclidean distance is used
to determine the change in future vectors. Assume $i$th ellipsoid generated at
time $t_{1}$ and $j$th ellipsoid is generated at time $t_{2}$, then the
distance between them is defined as:

$d_{i,j}=||fte(\mathcal{E}_{i})-fte(\mathcal{E}_{j})||_{2}.$ (8)

The $j$th ellipsoid is matched with the one from the previous frame with the
smallest distance. For motion estimation, we assume that the obstacles are
moving with constant linear and angular velocity. When the two ellipsoids are
matched, the linear and angular velocity of $j$th ellipsoid is calculated by a
simple numerical derivative.

$v_{j}=\frac{x_{c_{j}}-x_{c_{i}}}{t_{2}-t_{1}}\quad\omega_{j}=\frac{\theta_{j}-\theta_{i}}{t_{2}-t_{1}}$
(9)

This calculation is not directly accurate due to the fact that ellipsoids are
generated with an approximation algorithm so, a Kalman filter is integrated to
reduce error with more samples.

##  V PLANNING & CONTROL

We use a single model predictive control algorithm for both planning and
control. We use the simple point mass dynamics as the vehicle dynamics. The
states $\xi=\\{p_{x},p_{y},v_{x},v_{y}\\}$ are positions and velocities,
respectively. The inputs $u=\\{F_{x},F_{y}\\}$ are forces in the $x$ and $y$
directions. The state space model $\dot{\xi}=A\xi+Bu$ is defined as follows:

$\begin{bmatrix}\dot{p}_{x}\\\ \dot{p}_{y}\\\ \dot{v}_{x}\\\
\dot{v}_{y}\end{bmatrix}=\begin{bmatrix}0&0&1&0\\\ 0&0&0&1\\\ 0&0&0&0\\\
0&0&0&0\end{bmatrix}\begin{bmatrix}p_{x}\\\ p_{y}\\\ v_{x}\\\
v_{y}\end{bmatrix}+\begin{bmatrix}0&0\\\ 0&0\\\ 1/m&0\\\
0&1/m\end{bmatrix}\begin{bmatrix}F_{x}\\\ F_{y}\end{bmatrix},$ (10)

where $m$ is the mass of the point. The collision avoidance constraint is
defined similarly to [8]:

$c_{i}^{j}=(\Delta\xi_{i,j})^{T}R(\theta_{i,j})^{T}H(\alpha,\beta)R(\theta_{i,j})(\Delta\xi_{i,j})>1,$
(11)

where $\Delta\xi_{i,j}$ is the distance between the center of $j$th ellipsoid
and the vehicle at time $i$, $\alpha$, and $\beta$ are the semi-axes of an
enlarged ellipse defined in [8]. We use the model predictive controller
defined in [31] where soft constraints are used for collision avoidance.

[2] ξ, UJ^t+J^u+J^c (ξ_i+1-ξ_i)/dt=Aξ_i+Bu_i c_i^j+ψs_i¿1 ξ_min≤ξ_i ≤ξ_max
u_min ≤u_i ≤u_max ξ_0=ξ_init i={1,2,…,N},

where $J^{t}$ is tracing cost, $J^{u}$ control cost, $s_{k}$ is slack
variable, $J^{c}$ is the cost of the slack variable, $\psi$ is the sensitivity
of relaxed collision constraint, $N$ is the prediction horizon,
$\xi_{max},\xi_{min}$ are the state bounds, $u_{max},u_{min}$ are the control
bounds and $\xi_{init}$ is the initial state. Interior Point
Optimizer(IPOPT)[32] was used to solve the Nonlinear Programming(NLP) problem
defined in (V).

##  VI RESULTS

We test the proposed algorithm across five different scenarios. The parameters
used in these scenarios are given in Table I. In the first scenario shown in
Fig. 6 where there are only static obstacles, our algorithm generated six
ellipsoids that closely represent the ground truth map. Affinity
propagation[33] method also provided a suitable result, yet it generated 21
clusters and its computation time was more than five seconds, whereas our
algorithm generates six ellipsoids in less than 100 milliseconds for the same
map.

Parameter |  Value
---|---
$n_{ell_{max}}$ |  30
$e_{kyachiyan}$ |  0.05
$r_{ell_{threshold}}$ |  0.6
$\psi$ |  0.15
$N$ |  20
$Q$ |  eye(4)
$P$ |  diag(.1,.1)
$S$ |  eye($n_{ell}$)
$u_{max},u_{min}$ |  20,-20 N
$\xi_{max}$ |  $\infty,\infty,20,20$ m,m/s
$\xi_{min}$ |  $-\xi_{max}$
$m$ |  1 kg
TABLE I: Parameters for the test. $S$ soft constraint penalty matrice, $P$ is
the control penalty matrice and $Q$ is the tracking penalty matrice.

The second scenario presented in Fig. 3, is designed to test motion estimation
performance. A plus-shaped obstacle, with linear and angular velocities, was
considered, and its motion was estimated via ellipsoids. The error in motion
estimation is seen in Fig. 4 where the error becomes close to zero after a few
samples.

In the third scenario presented in Fig. 7, both obstacles are dynamic and move
with constant velocity. The fourth scenario provided in Fig. 8 involves two
static and one dynamic obstacle, while the final scenario illustrated in Fig.
9 involves four static and four dynamic obstacles. We combined some known
clustering methods with the Khachiyan algorithm and test them in these
scenarios. The GMM was not involved in this comparison due to the fact that it
gives the result of our algorithm with the right number of clusters. For the
algorithms that require the number of clusters, the optimal values which are:
$6$ for Map-1, $4$ for Map-3, $4$ for Map-4, and $12$ for Map-5 were used.
Other parameters were tested and the most suitable ones were selected. The
results are given in Table II.

In the first scenario, The Spectral[34] method provided the best performance;
however, the results were close to each other. Although the Spectral method
provides convenient results when the number of clusters is chosen correctly,
it beings to fail when the clusters overlap with each other. In scenario 3,
our algorithm and the Spectral algorithm obtained the same result, whereas
generated ellipsoids were also almost identical. In scenarios four and five,
the other methods failed to reach the target location as a consequence of the
wrong clustering. The blue plots in the figures show the motion generated by
using our algorithm. All computations were run on a desktop computer with
Intel core 7th generation i7 processor. The computation time of our algorithm
in test maps is given in Fig. 5.

Figure 3: Map-2. This scenario which includes 214 points is generated to test motion estimation performance where the obstacle moves with the velocity of 5 m/s in the x direction, 2 m/s in the y direction, and rotates $\pi$/2 rad/s in the counter-clockwise direction. The result can be seen in Fig. 4. Figure 4: This figures show motion estimation errors in Map-2 seen in Fig. 3. The left graph shows the error in linear velocity estimation and the right one is the error in angular velocity estimation. The error is the absolute value of the difference between the real value and the estimated value. Method | Time to reach final location(s)
---|---
Map-1 |  Map-3 |  Map-4 |  Map-5
Ours |  0.96 |  1.24 |  1.32 |  3.9
Spectral∗ |  0.9 |  1.24 |  - |  -
DBSCAN |  1.08 |  1.79 |  - |  -
Affinity |  0.92 |  1.48 |  - |  -
K-means∗ |  0.98 |  - |  - |  -
Mean-Shift |  0.97 |  - |  - |  -
Hierarchical∗ |  0.97 |  - |  - |  -
TABLE II: Planning comparison in custom maps. ∗ means the algorithm requires
the number of clusters and - means the vehicle did not reach the target
location. Figure 5: Computation time plot for each map. Map-1 includes 520
points, Map-2 includes 214 points, Map-3 includes 248 points, Map-4 includes
908 points and Map-5 includes 2418 points. Figure 6: Map-1. This scenario was
created using 520 points where all obstacles are static.  Figure 7: Map-3. In
this scenario, obstacles(248 points) are moving with a velocity of -5 m/s in
the x(horizontal) direction and the vehicle is on the left side of the
obstacle while the target location is on the right side of the obstacle.
Figure 8: Map-4. This scenario includes both static and dynamic obstacles
where the total map is generated by using 908 points. While the walls on the
upper and lower side do not move, the obstacle between the walls rotates with
$\pi$/2 rad/s in the counter-clockwise direction.  Figure 9: Map-5. This
scenario is generated by using 2418 points which include both static and
dynamic obstacles. The four walls are static and the four obstacles that are
inside the walls rotate with $\pi$/2 rad/s in the counter-clockwise direction.

##  VII CONCLUSION

In this work, we propose an ellipsoid-based obstacle identification method to
represent obstacles and estimate their motions without requiring the number of
clusters beforehand. We validate our method in different scenarios that
involve both static and dynamic obstacles including the rotating obstacles and
compare our approach with other methods. In addition, we showed the motion
estimation performances and the computational requirements of our algorithm.
Our algorithm was able to produce accurate outputs even in environments with
complex obstacles, outperforming the compared approaches.

## References

  * [1] H. Chen and P. Lu, “Real-time identification and avoidance of simultaneous static and dynamic obstacles on point cloud for uavs navigation,” _Robotics and Autonomous Systems_ , vol. 154, p. 104124, 2022. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0921889022000665
  * [2] O. Esrafilian and H. D. Taghirad, “Autonomous flight and obstacle avoidance of a quadrotor by monocular slam,” in _2016 4th International Conference on Robotics and Mechatronics (ICROM)_ , 2016, pp. 240–245.
  * [3] J. Mok, Y. Lee, S. Ko, I. Choi, and H. S. Choi, “Gaussian-mixture based potential field approach for uav collision avoidance,” in _2017 56th Annual Conference of the Society of Instrument and Control Engineers of Japan (SICE)_ , 2017, pp. 1316–1319.
  * [4] M. Corah, C. O’Meadhra, K. Goel, and N. Michael, “Communication-efficient planning and mapping for multi-robot exploration in large environments,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 2, pp. 1715–1721, 2019\.
  * [5] D. M. Blei and M. I. Jordan, “Variational inference for Dirichlet process mixtures,” _Bayesian Analysis_ , vol. 1, no. 1, pp. 121 – 143, 2006. [Online]. Available: https://doi.org/10.1214/06-BA104
  * [6] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, “A density-based algorithm for discovering clusters in large spatial databases with noise,” in _Proceedings of the Second International Conference on Knowledge Discovery and Data Mining_ , ser. KDD’96. AAAI Press, 1996, p. 226–231.
  * [7] J. Lin, H. Zhu, and J. Alonso-Mora, “Robust vision-based obstacle avoidance for micro aerial vehicles in dynamic environments,” 02 2020.
  * [8] B. Brito, B. Floor, L. Ferranti, and J. Alonso-Mora, “Model predictive contouring control for collision avoidance in unstructured dynamic environments,” _IEEE Robotics and Automation Letters_ , vol. PP, pp. 1–1, 07 2019.
  * [9] D. Falanga, K. Kleber, and D. Scaramuzza, “Dynamic obstacle avoidance for quadrotors with event cameras,” _Science Robotics_ , vol. 5, no. 40, p. eaaz9712, 2020. [Online]. Available: https://www.science.org/doi/abs/10.1126/scirobotics.aaz9712
  * [10] A. Chakravarthy and D. Ghose, “Collision cones for quadric surfaces,” _IEEE Transactions on Robotics_ , vol. 27, no. 6, pp. 1159–1166, 2011.
  * [11] R. Shioda and L. Tunçel, “Clustering via minimum volume ellipsoids,” _Computational Optimization and Applications_ , vol. 37, no. 3, pp. 247–295, Jul 2007. [Online]. Available: https://doi.org/10.1007/s10589-007-9024-1
  * [12] L. G. Khachiyan, “Rounding of polytopes in the real number model of computation,” _Mathematics of Operations Research_ , vol. 21, no. 2, pp. 2:307–2:320, 1996. [Online]. Available: https://doi.org/10.1287/moor.21.2.307
  * [13] J. MacQueen, “Some methods for classification and analysis of multivariate observations,” 1967.
  * [14] N. S. Altman, “An introduction to kernel and nearest-neighbor nonparametric regression,” _The American Statistician_ , vol. 46, no. 3, pp. 175–185, 1992\. [Online]. Available: http://www.jstor.org/stable/2685209
  * [15] J. Park and H. Baek, “Stereo vision based obstacle collision avoidance for a quadrotor using ellipsoidal bounding box and hierarchical clustering,” _Aerospace Science and Technology_ , vol. 103, p. 105882, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1270963820305642
  * [16] D. D. Pasqual, S. J. Sooriyaarachchi, and C. Gamage, “Dynamic obstacle detection and environmental modelling for agricultural drones,” in _2022 7th Asia-Pacific Conference on Intelligent Robot Systems (ACIRS)_ , 2022, pp. 70–75.
  * [17] S. Xenaki, K. Koutroumbas, and A. Rontogiannis, “Generalized adaptive possibilistic c-means clustering algorithm,” in _Proceedings of the 10th Hellenic Conference on Artificial Intelligence_ , ser. SETN ’18. New York, NY, USA: Association for Computing Machinery, 2018. [Online]. Available: https://doi.org/10.1145/3200947.3201012
  * [18] R. Krishnapuram and J. Keller, “A possibilistic approach to clustering,” _IEEE Transactions on Fuzzy Systems_ , vol. 1, no. 2, pp. 98–110, 1993.
  * [19] ——, “The possibilistic c-means algorithm: insights and recommendations,” _IEEE Transactions on Fuzzy Systems_ , vol. 4, no. 3, pp. 385–393, 1996.
  * [20] A. Dhawale, X. Yang, and N. Michael, “Reactive collision avoidance using real-time local gaussian mixture model maps,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2018, pp. 3545–3550.
  * [21] T. Eppenberger, G. Cesari, M. Dymczyk, R. Siegwart, and R. Dubé, “Leveraging stereo-camera data for real-time dynamic obstacle detection and tracking,” 2020\. [Online]. Available: https://arxiv.org/abs/2007.10743
  * [22] Y. Wang, J. Ji, Q. Wang, C. Xu, and F. Gao, “Autonomous flights in dynamic environments with onboard vision,” 2021. [Online]. Available: https://arxiv.org/abs/2103.05870
  * [23] R. D. Hjelm, A. Fedorov, S. Lavoie-Marchildon, K. Grewal, A. Trischler, and Y. Bengio, “Learning deep representations by mutual information estimation and maximization,” 08 2018.
  * [24] P. Bekker and J. van Essen, “Ml and gmm with concentrated instruments in the static panel data model,” _Econometric Reviews_ , vol. 39, no. 2, pp. 181–195, 2020. [Online]. Available: https://doi.org/10.1080/07474938.2019.1580946
  * [25] X. Li, Y. Hu, M. Li, and J. Zheng, “Fault diagnostics between different type of components: A transfer learning approach,” _Applied Soft Computing_ , vol. 86, p. 105950, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1568494619307318
  * [26] G. Wu, “Fast and scalable variational bayes estimation of spatial econometric models for gaussian data,” _Spatial Statistics_ , vol. 24, pp. 32–53, 2018\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2211675318300162
  * [27] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [28] M. J. Todd and E. A. Yıldırım, “On khachiyan’s algorithm for the computation of minimum-volume enclosing ellipsoids,” _Discrete Applied Mathematics_ , vol. 155, no. 13, pp. 1731–1744, 2007. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0166218X07000716
  * [29] S. Gottschalk, M. Lin, and D. Manocha, “Obbtree: A hierarchical structure for rapid interference detection,” _Computer Graphics_ , vol. 30, 10 1997.
  * [30] M. ApS, _MOSEK Fusion API for Python 10.0.20_ , 2022. [Online]. Available: https://docs.mosek.com/10.0/pythonfusion/index.html
  * [31] M. Castillo-Lopez, S. A. Sajadi-Alamdari, J. L. Sanchez-Lopez, M. A. Olivares-Mendez, and H. Voos, “Model predictive control for aerial collision avoidance in dynamic environments,” pp. 1–6, 2018.
  * [32] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical Programming_ , vol. 106, no. 1, pp. 25–57, Mar 2006. [Online]. Available: https://doi.org/10.1007/s10107-004-0559-y
  * [33] B. J. Frey and D. Dueck, “Clustering by passing messages between data points,” _Science_ , vol. 315, no. 5814, pp. 972–976, 2007. [Online]. Available: https://www.science.org/doi/abs/10.1126/science.1136800
  * [34] A. Ng, M. Jordan, and Y. Weiss, “On spectral clustering: Analysis and an algorithm,” in _Advances in Neural Information Processing Systems_ , T. Dietterich, S. Becker, and Z. Ghahramani, Eds., vol. 14. MIT Press, 2001. [Online]. Available: https://proceedings.neurips.cc/paper/2001/file/801272ee79cfde7fa5960571fee36b9b-Paper.pdf


# A polynomial algorithm for the most degree-central shortest path problem

Johnson Phosavanh  The University of Sydney  Dmytro Matsypura  The University
of Sydney

###### Abstract

The degree centrality of a node, defined to be the number of nodes adjacent to
it, can be used as a measure of importance of a node to the structure of a
network. This metric can be extended to paths in a network, where the degree
centrality of a path is defined to be the number of nodes adjacent to it. In
this paper, we reconsider the problem of finding the most degree-central
shortest path in an unweighted network. We propose a polynomial algorithm with
the worst-case running time of $O(|V|^{3}(\Delta(G))^{2})$, where $|V|$ is the
number of vertices in the network and $\Delta(G)$ is the maximum degree of the
graph. We conduct a numerical study of our algorithm on synthetic and real-
world networks and compare our results to the existing literature. In
addition, we show that the same problem is NP-hard when a weighted graph is
considered.

##  1 Introduction

###  1.1 Motivation

We extend the work of Matsypura et al., (2023), where the authors consider
problems of finding the most degree-central walk, trail, path, induced path,
and shortest path in a graph. The authors show that the problems are nested in
the sense that a shortest path is an instance of an induced path, an induced
path is an instance of a path, which is an instance of a trail, which is also
an instance of a walk. Further, all the problems considered in this study were
shown to be NP-hard, except for the problem of finding the most degree-central
shortest path, which can be solved in polynomial time using an algorithm based
on the breadth-first search technique. We will refer to this algorithm as the
MVP algorithm. For all the remaining problems, the authors developed mixed-
integer programs and proposed heuristic solutions as a warm start to improve
performance.

In this paper, we focus on the most degree-central shortest path problem and
propose a more efficient algorithm than the one in Matsypura et al., (2023).
The problem under consideration has many applications, including the design of
distribution channels and social network analysis with security implications.
In the context of distribution network design, a planner may be interested in
constructing a hub line network by connecting cities in a single path whilst
maximising the network’s reach (de Sá et al.,, 2015). During natural
disasters, authorities may also use this model to consider prioritising
locations to set up evacuation centres based on maximising the number of
affected communities that can be reached most efficiently. In social networks,
by considering each user as a node, we can use this model to find the most
influential users that can be reached most efficiently (Peng et al.,, 2018),
which, from an adversarial player point of view, can be used with great effect
when information and software with malicious intent are distributed.

###  1.2 Problem definition

Let $G=(V,E)$ be an unweighted graph with a set of vertices $V$ and edges $E$.
We define $\mathcal{N}(\langle i_{1},i_{2},\dots,i_{n}\rangle)$ to be the set
of neighbours of the path $\langle i_{1},i_{2},\dots,i_{n}\rangle$, and
$C(P)=|\mathcal{N}(P)|$, to be the centrality of path $P$. Following the
definitions provided by Matsypura et al., (2023), we define the neighbourhood
$\mathcal{N}(P)$ of a path $P$ to be the set of vertices adjacent to the path,
excluding the vertices on the path itself, i.e.,

$\mathcal{N}(P)=\\{v:(u,v)\in E,u\in P\\}\setminus\\{u:u\in P\\}.$

Note that when the path consists of a single node, this reduces to the
traditional definition of the degree centrality of a node.

Given a graph $G$, we aim to solve the following problem:

$\displaystyle\max\\{C(P):P\text{ is a shortest path between two nodes in
}G\\},$ (1)

i.e., find a path $P$ with the largest degree centrality under the condition
that $P$ is a shortest path between a pair of nodes in $G$.

###  1.3 Structure of the paper

The remaining paper is structured as follows: in Section 2, we present a
polynomial time algorithm that can be used to solve Problem (1). This is
followed by Section 4, where we conduct an extensive numerical study of the
proposed algorithm on synthetic and real-world graph instances and provide
insights into the results. Finally, in Section 5, we provide concluding
remarks.

##  2 Algorithm

Before formally introducing the algorithm, we prove the following lemma that
underpins many efficient algorithms for solving problems involving shortest
paths.

######  Lemma 1.

If $\langle s,p_{1},\ldots,p_{k},t\rangle$ is a shortest path from $s$ to $t$,
then $\langle s,p_{1},\ldots,p_{k}\rangle$ is a shortest path from $s$ to
$p_{k}$.

###### Proof.

Suppose there is a shorter path from $s$ to $p_{k}$, and let this be $\langle
s,\tilde{p}_{1},\ldots,\tilde{p}_{j},p_{k}\rangle$. This implies that $\langle
s,\tilde{p}_{1},\ldots,\tilde{p}_{j},p_{k},t\rangle$ would be a shorter path
from $s$ to $t$. However, this cannot be true as $\langle
s,p_{1},\ldots,p_{k},t\rangle$ is a shortest path from $s$ to $t$. ∎

We now extend the lemma above to consider the degree centrality of the path.

######  Lemma 2.

If $\langle s,p_{1},\ldots,p_{k-1},p_{k},t\rangle$ is a most degree-central
shortest path from $s$ to $t$, then $\langle s,p_{1},\ldots,p_{k-1}\rangle$ is
a most degree-central shortest path from $s$ to $p_{k-1}$.

###### Proof.

First, observe that vertex $t$ cannot share any neighbours with $\langle
s,p_{1},\ldots,p_{k-2}\rangle$. To see this, assume that vertex $p_{i}$ on the
path, where $i\leq k-2$, shares a neighbour $v$ with vertex $t$. Then $\langle
s,p_{1},\ldots,p_{k-1},p_{k},t\rangle$ cannot possibly be the most degree-
central shortest path between $s$ and $t$ as $\langle
s,p_{1},\ldots,p_{i},v,t\rangle$ would be a shorter path.

Now consider the graph in Figure 1 where vertex $t$ does not share any
neighbours with any node on path $\langle s,p_{1},\ldots,p_{k-2}\rangle$, but
with $p_{k-1}$ and $p_{k}$.

$s$$p_{1}$$p_{k-2}$$p_{k-1}$$p_{k}$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$$t$
Figure 1: Why $\langle s,p_{1},\ldots,p_{k-1}\rangle$ must be a most degree-
central shortest path between $s$ and $p_{k-1}$ if $\mathcal{N}(\langle
s,p_{1},\ldots,p_{k},t\rangle)$ is a most degree-central shortest path between
$s$ and $t$

From the diagram, it is clear that if $\langle
s,p_{1},\ldots,p_{k-1},p_{k},t\rangle$ is a most degree-central shortest path
from $s$ to $t$, then $\langle s,p_{1},\ldots,p_{k-1}\rangle$ is a most
degree-central shortest path from $s$ to $p_{k-1}$. Now, suppose that there is
an alternate shortest path $\langle
s,\tilde{p}_{1},\ldots,\tilde{p}_{k-2},p_{k-1}\rangle$ to vertex $p_{k-1}$
that shares no neighbours with vertex $t$ and has a higher centrality. In this
case, path $\langle
s,\tilde{p}_{1},\ldots,\tilde{p}_{k-2},p_{k-1},p_{k},t\rangle$, would have a
centrality greater than that of path $\langle
s,p_{1},\ldots,p_{k-1},p_{k},t\rangle$ as vertices $p_{k}$ and $t$ contribute
the same vertices to the neighbourhoods of both paths. ∎

We now introduce the following notation to simplify the exposition. First,
consider a source node $s$. Let $d_{u}$ and $P_{u}$ be the shortest distance
and most degree-central shortest path between nodes $s$ and $u$, respectively,
and let $\mathcal{P}_{u}$ be the set of vertices that can directly precede
vertex $u$ on any shortest path from $s$ to $u$. Additionally, we introduce
the operation $\langle P,u,v\rangle$, which we will use to denote extending
the path $P$ to vertices $u$, then $v$.

To solve Problem (1), we modify Dijkstra’s algorithm (Dijkstra,, 1959).
Dijkstra’s algorithm can be used to identify the shortest path from a source
node $s$ to any other node $t$ in the network by considering nodes already
connected to node $s$ from previous iterations. At each stage of the
algorithm, we consider all vertices $S\subseteq V$ such that $u\in S$ if $u$
has been connected to $s$ from a previous iteration. We then consider each
vertex $u\in S$ in increasing order of their path length from the source node
$s$ and each edge $\\{(u,v):(u,v)\in E\\}$, i.e., all the edges that can be
connected to $u$. Then, we compare the length of the previously calculated
shortest path from $s$ to $v$ and a new path from $s$ to $u$ and then to $v$
and take the minimum (this is the traditional Dijkstra’s algorithm). To obtain
the most degree-central shortest path, we modify the algorithm in the
following manner: for each vertex $v$, we place all nodes that can directly
precede vertex $v$ on the shortest path from $s$ to $v$ in the set
$\mathcal{P}_{v}$. Each time we consider extending the shortest path $\langle
s,\dots,u\rangle$ to vertex $v$, we compare the neighbourhood of the existing
most degree-central shortest path from $s$ to vertex $v$ against the
neighbourhoods of paths obtained by extending the most degree-central shortest
path to $w$, $P_{w}$, to $u$ then $v$ for all $w\in\mathcal{P}_{u}$.
Furthermore, instead of only considering the case where a path is shorter in
length than a previous path in the traditional Dijkstra’s algorithm, we also
consider tied shortest paths as we need to compute and compare the centrality
to ensure the most degree-central path is found.

This is summarised in Algorithm 1.

Algorithm 1 Finding the most degree-central shortest path

1:Graph $G=(V,E)$, starting vertex $s$.

2:Most degree-central shortest paths $P_{v}$ from node $s$ to $v$ for $v\in
V$.

3:for $v\in V$ do

4: $d_{v}\leftarrow+\infty$, $P_{v}\leftarrow\text{undefined}$,
$\mathcal{P}_{v}=\emptyset$.

5:end for

6:$P_{s}\leftarrow\langle s\rangle$.

7:Insert $(0,s)$ to the priority queue $Q$.

8:while $Q$ not empty do

9: $(d_{u},u)\leftarrow$ pop minimum element of $Q$.

10: for $v\in\\{v^{\prime}\in V:(u,v^{\prime})\in E$} do

11: $d_{\text{new}}\leftarrow d_{u}+1$.

12: if $d_{\text{new}}=1$ then

13: Insert $(d_{\text{new}},v)$ to the priority queue $Q$.

14: $d_{v}\leftarrow d_{\text{new}}$, $P_{v}=\langle s,v\rangle$,
$\mathcal{P}_{v}\leftarrow\\{s\\}$.

15: else if $d_{\text{new}}<d_{v}$ then

16: Insert $(d_{\text{new}},v)$ to the priority queue $Q$.

17: $d_{v}\leftarrow d_{\text{new}}$,
$\mathcal{P}_{v}\leftarrow\mathcal{P}_{v}\cup\\{u\\}$.

18: for $w\in\mathcal{P}_{u}$ do

19: if $C(\langle P_{w},u,v\rangle)>C(P_{v})$ then

20: $P_{v}\leftarrow\langle P_{w},u,v\rangle$.

21: end if

22: end for

23: else if $d_{\text{new}}=d_{v}$ then

24: $\mathcal{P}_{v}\leftarrow\mathcal{P}_{v}\cup\\{u\\}$.

25: for $w\in\mathcal{P}_{u}$ do

26: if $C(\langle P_{w},u,v\rangle)>C(P_{v})$ then

27: $P_{v}\leftarrow\langle P_{w},u,v\rangle$.

28: end if

29: end for

30: end if

31: end for

32:end while

######  Proposition 1.

Given a graph $G=(V,E)$ and a source vertex $s\in V$, Algorithm 1 returns the
most degree-central shortest paths from $s$ to all other vertices in $V$.

###### Proof.

We show this by induction.

Base case: The most degree-central shortest path to the source vertex $s$ is
just $\langle s\rangle$. Also, it is easy to see that any path connected in
the first iteration of the while loop must be the shortest path from $s$ as it
is directly connected to the source node. Furthermore, this will be the most
degree-central path as there cannot be another path of the same or shorter
length to this set of vertices.

Inductive step: Assume that in the $k$th iteration of the while loop, the most
degree-central shortest paths to all nodes in $S\subseteq V$ have been found.

In the $(k+1)$th iteration, we consider node $u$ where $d_{u}$ is the smallest
distance between the source vertex $s$ and all other previously connected
nodes that have not yet been processed. The for-loop in line 8 allows us to
branch to all adjacent nodes, while lines 10, 13 and 21 guarantee that only
the shortest path will be recorded. The for-loops and if-statements in lines
16-20 and 23-27 then allow us to extend paths $P_{w}$ to $u$ then $v$ for all
$w\in\mathcal{P}_{u}$, which, by the inductive hypothesis and Lemma 2, will
return the most degree-central shortest path from $s$ to $v$ that passes
through vertex $u$. After considering extending all shortest paths of length
$d_{u}$ in subsequent iterations of the while-loop, and therefore considering
all possible $u\in\mathcal{P}_{v}$, we are guaranteed the most degree-central
shortest path from $s$ to $v$. ∎

######  Proposition 2.

Given a graph $G=(V,E)$, the worst-case running time of Algorithm 1 is
$O(|V|^{4})$.

###### Proof.

Algorithm 1 is a modification of Dijkstra’s algorithm. Dijkstra’s algorithm
(with a priority queue) has a runtime of $O(|V|+|E|\log|V|)\sim O(|E|\log|V|)$
as initialising $d$ and $P$ take $O(|V|)$ time and at most, there are $|E|$
iterations of the for-loop, each taking $O(\log|V|)$ time to add to the
priority queue $Q$ (Johnson,, 1977). However, for an unweighted graph, this
can be reduced to $O(|V|\log(|V|))$ time as we only add to the priority queue
at most $|V|$ times.

Assuming that $\mathcal{N}(P_{w})$ and $\mathcal{N}(P_{v})$ have been stored
and can be easily accessed at each iteration of the while loop, computing
$\mathcal{N}(\langle P_{w},u,v\rangle)$ in the inner for-loop for all possible
$w$’s takes $O(|V|^{2})$ time in the worst case as we need to consider at most
$|V|$ $P_{w}$’s, and computing $C(\langle P_{w},u,v\rangle)$ takes $O(|V|)$
time, which dominates the running time of the inner for-loop. In total, this
operation is completed at most $O(|V|^{2})$ times as there are $|V|$
iterations of the while-loop and $|V|$ iterations of the outer for-loop as
each vertex can have at most $|V|$ connected edges. As a result, Algorithm 1
has runtime of $O(|V|^{4})$. ∎

This bound can be tightened. Let $\Delta(G)$ be the maximum degree of graph
$G$. Then, the following is true.

######  Proposition 3.

Given a graph $G=(V,E)$, the worst-case running time of Algorithm 1 is
$O(|V|^{2}(\Delta(G))^{2})$.

###### Proof.

This tighter upper bound of the running time can be obtained for the following
reasons. For each iteration of the while-loop, we can only have at most
$\Delta(G)$ iterations of the outer for-loop as we can only consider extending
the path to node $u$ to one of its neighbours and not all the vertices in the
graph. Additionally, there are only $|\mathcal{P}_{u}|\leq\Delta(G)$
iterations of the inner for-loops in lines 16-20 and 23-27 as there are at
most $\Delta(G)$ nodes that may directly precede vertex $u$ on a shortest path
to vertex $u$. As a result, two orders of $|V|$ can be replaced with
$\Delta(G)$ from the runtime quoted in Proposition 2. ∎

######  Theorem 1.

Problem (1) can be solved in $O(|V|^{3}(\Delta(G))^{2})$ time.

###### Proof.

By Proposition 3, Algorithm 1 can be used to solve the problem of finding all
most degree-central shortest paths from a given source node in
$O(|V|^{2}(\Delta(G))^{2})$ time. To solve Problem (1), we can apply Algorithm
1 over all vertices in the network, resulting in an overall running time of
$O(|V|^{3}(\Delta(G))^{2})$. ∎

Comparing this to the running time of Algorithm MVP, which has a runtime of
$O(k|V|^{6})\sim O(|V|^{7})$, with $k$ being the diameter of the graph,
Algorithm 1 has a runtime of $O(|V|^{3}(\Delta(G))^{2})\sim O(|V|^{5})$, which
is an improvement of $O(|V|^{2})$, and in addition, we only require
$O(|V|^{2})$ space, rather than $O(|V|^{3})$.

##  3 The case of weighted graphs

We first show that the Problem (1) is NP-hard even with just two distinct
weights by showing it is at least as hard as the Maximum 2-Satisfiability (Max
SAT2) problem, which is known to be NP-complete (Garey and Johnson,, 1979).

Max SAT2
---
Instance: |  Set of $U$ variables, collection $C$ of clauses over $U$ such that each clause $c\in C$ has $|c|=2$ and positive integer $k\leq|C|$.
Question: |  Is there a truth assignment for $U$ that simultaneously satisfies at least $k$ of the clauses in $C$?

Before we show the reduction, we first formally define the decision variant of
Problem (1) for a weighted graph.

MDCWSP-D
---
Instance: |  A weighted graph $G$ and a positive integer $\alpha$.
Question: |  Is there a shortest path between two vertices in $G$ such that its centrality is greater than $\alpha$?

######  Theorem 2.

MDCWSP-D is NP-complete.

###### Proof.

First, observe that the centrality of a path can be computed in polynomial
time; hence MDCWSP-D is in NP.

Given an instance of Max SAT2, we can construct a directed weighted graph as
follows. First, consider an unweighted graph $H$ with $2|U|$ nodes and
$4(|U|-1)$ edges with the following structure: we start with two unlinked
paths $\langle x_{1},\dots,x_{|U|}\rangle$ and
$\langle\bar{x}_{1},\dots,\bar{x}_{|U|}\rangle$ both of length $|U|-1$. For
$i=1,\dots,|U|-1$, we add additional edges $(x_{i},\bar{x}_{i+1})$ and
$(\bar{x}_{i},x_{i+1})$ to $H$. This means that the shortest path between any
nodes $x_{i}$ or $\bar{x}_{i}$ and $x_{j}$ or $\bar{x}_{j}$ is of length $j-i$
for $i<j$. Finally, to complete $H$, we add four additional nodes $s$,
$\bar{s}$ that are both connected to $x_{1}$ and $\bar{x}_{1}$, and $t$ and
$\bar{t}$ that are connected to $x_{|U|}$ and $\bar{x}_{|U|}$. Note that we
connect the nodes such that the edges point towards $s$, $\bar{s}$ $t$ and
$\bar{t}$. To obtain the weighted graph $G$, we add vertices and edges to $H$
depending on the clauses. For each clause $c\in C$, we add a vertex $y_{c}$ to
$G$ and undirected edges $(z,y_{c})$ for each variable $z$ in clause $c$, each
with edge weight $|U|$. Given these edge weights, any shortest path in the $G$
cannot contain any vertex $y_{c}$ in the middle of the path because an
alternate shorter path will exist. Note: The additional edges can be directed;
however, we set them to be undirected to amplify the difference between the
unweighted and weighted settings.

An example of the structure of $G$ is given in Figure 2.

$s$$\bar{s}$$\bar{x}_{1}$$x_{1}$$\bar{x}_{2}$$x_{2}$$\bar{x}_{3}$$x_{3}$$\bar{x}_{4}$$x_{4}$$\bar{x}_{|U|-1}$$x_{|U|-1}$$\bar{x}_{|U|}$$x_{|U|}$$t$$\bar{t}$$y_{1}$$y_{2}$$y_{3}$$y_{4}$$y_{|C|-1}$$y_{|C|}$………
Figure 2: Example of a weighted graph $G$ used in Theorem 2, with the blue
vertices and edges indicating the initial unweighted graph $H$ and the red
vertices and edges indicating the additional vertices and edges with edge
weight $|U|$.

We now show that there exists a truth assignment for $U$ with more than $k$
satisfiable clauses if and only if there exists a shortest path in $G$ with a
degree of at least $|U|+k+4$.

$\Longrightarrow$ If there is some truth assignment for $U$ with at least $k$
satisfiable clauses, then the corresponding path in $G$ will have at least
$|U|+k+4$ neighbours by construction as we have $|U|$ neighbours from the
False assignments and $k$ neighbours from the satisfied clauses, with an
additional $4$ neighbours coming from the source and sink vertices.

$\Longleftarrow$ Suppose we are given a shortest path in $G$ with at least
$|U|+k+4$ neighbours. We can transform this path into the form
$\langle(x_{1}\oplus\bar{x}_{1}),\dots,(x_{U}\oplus\bar{x}_{|U|})\rangle\>$
and maintain a degree centrality of at least $|U|+k+4$. Consider the following
cases:

  * •

If the given path ends at $s$, $\bar{s}$, $t$ or $\bar{t}$, we remove this
node, and the centrality will increase by one.

  * •

If the given path has starting and end nodes that are both in $H$, then any
arbitrary extension of the path so that it starts at
$(x_{1}\oplus\bar{x}_{1})$ and ends at $(x_{|U|}\oplus\bar{x}_{|U|})$ will
maintain or increase the centrality as the neighbourhood will contain the
unselected path nodes.

  * •

If the path starts or ends at one of the additional nodes $y_{c}$, for
example, $\langle y_{c},x_{i},\dots,x_{j}\rangle$, we can remove $y_{c}$ from
the path and extend the path using the result from the above case, resulting
in the path
$\langle(x_{1}\oplus\bar{x}_{1}),\dots,(x_{i-1}\oplus\bar{x}_{i-1}),x_{i},\dots,x_{j},(x_{j+1}\oplus\bar{x}_{j+1}),\dots,(x_{|U|}\oplus\bar{x}_{|U|})\rangle$.
This new path will not have a centrality smaller than the original path as we
will gain $y_{c}$ as a neighbour.

In all the above cases, the new paths can be computed in polynomial time, and
the path returned will be a solution to Max SAT2 because after the removal of
the two start and end nodes, as well as the unvisited $|U|$ nodes in $H$, the
remaining $k$ nodes in the neighbourhood correspond to the satisfied clauses
and the path corresponds to the truth assignment. ∎

######  Corollary 1.

Problem (1) with a weighted graph is NP-hard.

Although Dijkstra’s algorithm can be used to solve the problem of finding the
shortest path in a graph with positive weighted edges, we note that Algorithm
1 does not guarantee an optimal solution as Lemma 2 does not hold. However, if
the weights are positive integers, we can augment the network with auxiliary
vertices and edges to represent the weighted edges. Then Algorithm 1 can be
applied with the worst-case running time of
$O(w_{\text{sum}}|V|^{2}(\Delta(G))^{2})$, where $w_{\text{sum}}$ is the sum
of all edge weights. Note that the running time is not
$O(w_{\text{sum}}^{3}(\Delta(G))^{2})$ as we only need to change the number of
times the while-loop is executed. The remaining runtime remains the same as we
only need to consider setting the starting vertex to be one that is in the
original graph, and we only need to compute and compare the centrality of the
paths when we consider a vertex in the original graph. By applying Algorithm 1
on the augmented graph, the runtime is now pseudo-polynomial, as it depends on
the sum of the edge weights. However, if all the edge weights $w_{uv}$ are
generated from some positive continuous distribution $W_{uv}$ for all
$(u,v)\in E$, the shortest path will be unique with probability 1 as
$P(W_{uv}=w_{uv})=0$ for all $(u,v)\in E$. In this case, the problem can be
solved by first solving the all pairs shortest path problem and then
evaluating and comparing the centralities of the paths. This can be done in
$O(|V|^{3})$ time with the Floyd-Warshall algorithm.

##  4 Computational experiments

We tested Algorithm 1 on a large number of synthetic and real-world graphs.
The experiments were conducted on a Mac Studio (2023) with Apple M2 Ultra and
128 GB of RAM running macOS Ventura 13.6. The algorithm was implemented in
Python 3.9.

###  4.1 Data and results

The synthetic graphs were randomly generated using two procedures outlined
below (each type of graph was generated 30 times for each parameter setting):

  * •

Watts–Strogatz graphs. We considered the ‘small-world’ network proposed by
Watts and Strogatz, (1998), which represents networks that are “highly
clustered, like regular lattices, yet have small characteristic path lengths”.
We considered Watts–Strogatz graphs with 100, 500, 1000, 5000, and 10000
nodes. The number of initial neighbours was set to 4, and the rewiring
probability was set to one of 0.1 or 0.2.

  * •

Barabási-Albert graphs. The ‘scale-free preferential attachment’ model
proposed by Barabási and Albert, (1999) generates graphs with edges whose
degree follows the power rule. We considered graphs with the same number of
nodes as the Watts–Strogatz instances and set the number of edges attached to
any new node to 2 to mirror the edge density of the Watts-Strogatz graphs.

In total, we generated and analysed 450 synthetic instances. We averaged the
results for each parameter instance and summarised them in Tables 1 and 2.

To compare our results to Algorithm MVP, we also considered the following
real-world instances:

  * •

Krebs: The terrorist network of the 9/11 hijackers (Krebs,, 2002).

  * •

Dolphins: A network representing the social interactions of bottlenose
dolphins (Lusseau et al.,, 2003).

  * •

Sandi Auths: An academic collaboration network (Rossi and Ahmed,, 2015).

  * •

IEEE Bus: A subnetwork of the US Electric Power System in 1962 (Davis and Hu,,
2011).

  * •

Santa Fe: A collaboration network for the Santa Fe Institute (Girvan and
Newman,, 2002).

  * •

US Air 97: A transportation network of US Air in 1997 (Batagelj and Mrvar,,
2006; Davis and Hu,, 2011).

  * •

Bus: Bus power system network (Davis and Hu,, 2011).

  * •

Email: Network of e-mail interchanges (Guimera et al.,, 2003; Davis and Hu,,
2011).

  * •

Cerevisiae: Biological network of yeast protein-protein interactions (Davis
and Hu,, 2011).

The results from the real-world graph instances are presented in Table 3.

| Watts–Strogatz (4, 0.1) | Watts–Strogatz (4, 0.2)
---|---|---
$|V|$ | 100 | 500 | 1000 | 5000 | 10000 | 100 | 500 | 1000 | 5000 | 10000
$|E|$ | 200 | 1000 | 2000 | 10000 | 20000 | 200 | 1000 | 2000 | 10000 | 20000
$\Delta(G)$ | 5.87 | 6.40 | 6.70 | 7.27 | 7.33 | 6.60 | 7.40 | 7.43 | 8.03 | 8.20
diam | 10.43 | 15.40 | 17.73 | 22.33 | 24.50 | 8.43 | 11.50 | 12.97 | 16.13 | 17.47
diam centrality | 22.23 | 31.20 | 35.37 | 42.53 | 45.20 | 21.67 | 28.80 | 30.60 | 37.73 | 40.17
path length | 9.03 | 12.97 | 14.73 | 18.63 | 19.37 | 7.03 | 9.33 | 10.17 | 12.53 | 13.27
path centrality | 23.17 | 33.87 | 37.60 | 46.90 | 50.40 | 23.33 | 31.73 | 35.17 | 43.77 | 47.07
MVP runtime | 0.41 | 59.73 | 513.84 | - | - | 0.41 | 55.73 | 482.51 | - | -
Alg 1 runtime | 0.07 | 2.08 | 10.26 | 366.88 | 1712.52 | 0.06 | 2.08 | 9.75 | 359.95 | 1704.02
Table 1: Results for Watts–Strogatz graphs. | Barabási-Albert
---|---
$|V|$ | 100 | 500 | 1000 | 5000 | 10000
$|E|$ | 196 | 996 | 1996 | 9996 | 19996
$\Delta(G)$ | 25.27 | 53.83 | 76.53 | 188.20 | 257.40
diam | 5.57 | 7.03 | 7.27 | 8.53 | 9.00
diam centrality | 46.17 | 120.10 | 179.30 | 382.57 | 553.27
path length | 3.87 | 4.67 | 5.07 | 5.80 | 6.03
path centrality | 51.87 | 138.33 | 199.83 | 483.93 | 676.47
MVP runtime | 0.36 | 45.92 | 381.45 | - | -
Alg 1 runtime | 0.07 | 2.31 | 13.46 | 713.06 | 3792.77
Table 2: Results for Barabási-Albert graphs. | Kreb | Dolphins | Sandi Auths | IEEE Bus | Santa Fe | US Air 97 | Bus | Email | Cerevisiae
---|---|---|---|---|---|---|---|---|---
$|V|$ | 62 | 62 | 86 | 118 | 118 | 332 | 662 | 1133 | 1458
$|E|$ | 153 | 159 | 124 | 179 | 200 | 2126 | 906 | 5451 | 1948
$\Delta(G)$ | 22 | 12 | 12 | 9 | 29 | 139 | 9 | 71 | 56
diam | 5 | 8 | 11 | 14 | 12 | 6 | 25 | 8 | 19
diam centrality | 38 | 30 | 32 | 32 | 90 | 167 | 45 | 159 | 57
path length | 3 | 6 | 7 | 8 | 10 | 3 | 20 | 4 | 7
path centrality | 40 | 36 | 34 | 33 | 92 | 206 | 50 | 187 | 156
MVP runtime | 0.09 | 0.11 | 0.25 | 0.73 | 0.66 | 21.96 | 147.03 | 807.92 | 1317.34
Alg 1 runtime | 0.03 | 0.03 | 0.04 | 0.09 | 0.09 | 2.17 | 3.88 | 26.53 | 22.89
Table 3: Results for real-world instances

###  4.2 Discussion and insights

We observe that by using Algorithm 1, we can now solve problems much more
efficiently than Algorithm MVP. For instance, the Watts-Strogatz instances
with 1000 nodes, which took around 8 minutes to solve, can now be solved in
around 10 seconds, and problems involving 5000 and 10000 nodes can now be
solved within a fifth and half of an hour, respectively. A similar observation
can be made for the Barabási-Albert graphs. We did not consider solving larger
instances with Algorithm MVP, as one 5000-node instance took over 13 hours to
solve. This is consistent with the expected runtime based on the difference
between solving 100 and 500-node instances. For real-world instances, the
improvement in runtime is also significant. For example, the Cerevisiae
instance no longer requires over twenty minutes to solve; instead, it can be
solved in just over twenty seconds. We also note that in practice, Algorithm 1
scales much more efficiently than stated in Theorem 1 because, for most
vertices $u$ and paths $P_{u}$, the size of $\mathcal{P}_{u}$ and
$\mathcal{N}(P_{u})$ is small as the number of immediate neighbours is smaller
than the maximum degree of the graph/vertices in the graph. We also observe
the difference in running time over the two types of synthetic instances.
Although the number of vertices and edges in both sets of graphs are nearly
identical, the differences in running times for both algorithms on both sets
of graphs are pronounced. Algorithm 1 returns a solution faster in Watts-
Strogatz graphs than Barabási-Albert graphs of similar sizes, while the
converse is true for Algorithm MVP. This is consistent with the theoretical
running times of Algorithm MVP and Theorem 1 as Watts-Strogatz graphs have
smaller maximum degrees but larger diameters than Barabási-Albert graphs.
Lastly, we note that although in our experiments, both Algorithm MVP and
Algorithm 1 were run sequentially over all starting vertices, both could be
run in parallel for an improved runtime.

##  5 Concluding remarks

In this paper, we consider the problem of finding the most degree-central
shortest path in a graph. We show that the problem with an unweighted graph
can be solved in polynomial time by modifying Dijkstra’s algorithm and that
this new algorithm is more efficient than the one proposed by Matsypura et
al., (2023). We then test our algorithm numerically on synthetic and real-
world graph instances to examine the improvement. In addition to this, we show
that the problem with a weighted graph is NP-hard, even with two distinct
weights, and we provide a pseudo-polynomial time extension of our algorithm
for the integer weighted problem, as well as a polynomial time solution to the
case with continuous weights.

## References

  * Barabási and Albert, (1999) Barabási, A.-L. and Albert, R. (1999).  Emergence of scaling in random networks.  Science (American Association for the Advancement of Science), 286(5439):509–512.
  * Batagelj and Mrvar, (2006) Batagelj, V. and Mrvar, A. (2006).  Pajek datasets.
  * Davis and Hu, (2011) Davis, T. A. and Hu, Y. (2011).  The University of Florida Sparse Matrix Collection.  ACM Transactions on Mathematical Software (TOMS), 38(1).
  * de Sá et al., (2015) de Sá, E. M., Contreras, I., Cordeau, J.-F., de Camargo, R. S., and de Miranda, G. (2015).  The hub line location problem.  Transportation Science, 49(3):500–518.
  * Dijkstra, (1959) Dijkstra, E. W. (1959).  A note on two problems in connexion with graphs.  Numerische Mathematik, 1(1):269–271.
  * Garey and Johnson, (1979) Garey, M. R. and Johnson, D. S. (1979).  Computers and intractability: a guide to the theory of NP-completeness.  A Series of books in the mathematical sciences. W. H. Freeman, San Francisco.
  * Girvan and Newman, (2002) Girvan, M. and Newman, M. E. (2002).  Community structure in social and biological networks.  Proceedings of the National Academy of Sciences, 99(12):7821–7826.
  * Guimera et al., (2003) Guimera, R., Danon, L., Diaz-Guilera, A., Giralt, F., and Arenas, A. (2003).  Self-similar community structure in a network of human interactions.  Physical Review E, 68(6):065103.
  * Johnson, (1977) Johnson, D. B. (1977).  Efficient algorithms for shortest paths in sparse networks.  Journal of the ACM, 24(1):1–13.
  * Krebs, (2002) Krebs, V. (2002).  Uncloaking terrorist networks.  First Monday, 7(4).
  * Lusseau et al., (2003) Lusseau, D., Schneider, K., Boisseau, O. J., Haase, P., Slooten, E., and Dawson, S. M. (2003).  The bottlenose dolphin community of doubtful sound features a large proportion of long-lasting associations.  Behavioral Ecology and Sociobiology, 54(4):396–405.
  * Matsypura et al., (2023) Matsypura, D., Veremyev, A., Pasiliao, E. L., and Prokopyev, O. A. (2023).  Finding the most degree-central walks and paths in a graph: Exact and heuristic approaches.  European Journal of Operational Research, 308(3):1021–1036.
  * Peng et al., (2018) Peng, S., Zhou, Y., Cao, L., Yu, S., Niu, J., and Jia, W. (2018).  Influence analysis in social networks: A survey.  Journal of Network and Computer Applications, 106:17–32.
  * Rossi and Ahmed, (2015) Rossi, R. and Ahmed, N. (2015).  The network data repository with interactive graph analytics and visualization.  Proceedings of the AAAI Conference on Artificial Intelligence, 29(1).
  * Watts and Strogatz, (1998) Watts, D. J. and Strogatz, S. H. (1998).  Collective dynamics of ‘small-world’ networks.  Nature (London), 393(6684):440–442.