2.8cm2.8cm2.8cm2.8cm

# Constructing entire functions with non-locally connected Julia set by
quasiconformal surgery

Yanhua Zhang & Gaofei Zhang  School of Mathematical Sciences, Qufu Normal
University, Qufu 273165, P. R. China  g<EMAIL_ADDRESS>Department of
Mathematics, Nanjing University, Nanjing 210093, P. R. China
<EMAIL_ADDRESS>
###### Abstract.

We give an alternative way to construct an entire function with quasiconformal
surgery so that all its Fatou components are quasi-circles but the Julia set
is non-locally connected.

###### Key words and phrases:

entire functions, non-locally connected Julia sets

###### 2010 Mathematics Subject Classification:

Primary 37F45; Secondary 37F20, 37F10

##  1\. Statements of the main results

In a recent manuscript [1] the authors there construct an interesting example
of entire function for which all the Fatou components are quasi-circles but
the Julia set is non-locally connected. The idea used in their construction
follows Eremenco and Sodin [4] and is called Maclane-Vinberg method. The
purpose of this article is to provide an alternative way to construct such
functions by quasiconformal surgery. Our construction is very much inspired by
[2]. Compared with the method used in [1], which is purely function theoretic,
ours has more geometric feature.

The starting point of our construction is a quasi-regular map $F$ which is
described as follows. The building blocks of $F$ are the unit disk
$\mathbb{D}$ and a sequence of squares $\\{P_{n}\\}_{n\geq 1}$ with diameters
being even integers. We assemble these objects from the right to the left so
that $\mathbb{D}$ is the rightmost one. See Figure 1 for an illustration.
Recall that a map $f:\mathbb{C}\to\mathbb{C}$ is called $K$-quasi-regular if
$f=g\circ\phi$ for some entire function $g$ and $K$-quasiconformal map $\phi$
of the plane.

Figure 1. Constructing the quasi-regular map $F$

######  Main Lemma.

There is a $K>1$ such that for any $\\{P_{n}\\}_{n\geq 1}$ there is a
$K$-quasi-regular branched covering $F:\mathbb{C}\to\mathbb{C}$ satisfying
$F^{-1}(\mathbb{D})=\mathbb{D}\cup\bigcup_{n\geq 1}P_{n}$, and moreover,

  * •

$F:\mathbb{D}\to\mathbb{D}$ is the square map, and $F$ is symmetric about the
real axis, that is, $F(\bar{z})=\overline{F(z)}$ for all $z\in\mathbb{C}$,

  * •

$F:P_{n}\to\mathbb{D}$ is a branched covering of degree $d_{n}$ with $d_{n}$
being the diameter of $P_{n}$, and moreover, each $P_{n}$ contains exactly one
critical point in its interior which is the center and is mapped to the
origin,

  * •

the set of all the other critical points consists of the vertices of the
squares which belong to the negative real axis, which are all mapped to the
point $1$,

  * •

$F$ has no asymptotic values and has exactly two critical values $0$ and $1$,

  * •

For any compact subset $X\subset\mathbb{C}$, there is a $B>0$ depending only
on $X$ such that $|F(z)|<B$ for all $z\in X$.

As a consequence we have

######  Theorem 1.1.

There is an entire function $f$ with no asymptotic values and with exactly two
critical values $0$ and $1$ such that $f(0)=0$ and $f(1)=1$, and moreover,

  * •

$f$ has a critical point at $0$, and all the other critical points of $f$ are
negative reals and are mapped either to $0$ or to $1$,

  * •

all the Fatou components of $f$ are quasi-disks and are eventually mapped to
the super-attracting basin at the origin,

  * •

the Julia set of $f$ is non-locally connected.

Just after the submission of the paper, L. Rempe sent us the following
comments on this work. The entire function constructed here is subhyperbolic,
has one super-attracting cycle and is unbounded on the real axis, while the
one constructed in [1] is hyperbolic, has two super-attracting cycles and is
bounded on the real axis. The difference is not very significant however. This
is because the function in [1] can be used to obtain the function here through
a square root transformation, and on the other hand, the method here can also
be used to obtain the function in [1] through an easy adaptation.

##  2\. Proof of Theorem 1.1 assuming the main lemma

Let $\mathbb{C}$, $\mathbb{C}^{*}$, $\mathbb{T}$ and $\mathbb{D}$ denote the
complex plane, the punctured plane, the unit circle and the unit disk
respectively. Let $F$ be the quasi-regular branched covering map of the plane
guaranteed by the main lemma. Let $\mu$ be the pull back of the standard
complex structure by $F$. By the symmetry of $F$ it follows that $\mu$ is
symmetric about the real axis. Let $\psi:\mathbb{C}\to\mathbb{C}$ be the
quasiconformal homeomorphism which solves the Beltrami equation
$\psi_{\bar{z}}=\mu\psi_{z}$ and fixes $-1$ and $1$. Then $F\circ\psi^{-1}$ is
an entire function. By the symmetry of $\mu$ we have
$\psi(\bar{z})=\overline{\psi(z)}$. Then $\psi(\mathbb{T})$ is a quasi-circle
which is symmetric about the real axis. Let $\eta$ be the Riemann isomorphism
which maps the exterior of $\mathbb{T}$ to the exterior of $\psi(\mathbb{T})$
and fixes $\infty$ and $1$. By the symmetry of $\psi(\mathbb{T})$,
$\eta(\bar{z})=\overline{\eta(z)}$. Note that this implies that $\eta(-1)=-1$.
Let

$G=F\circ\psi^{-1}\circ\eta.$

Then $G$ is holomorphic in the exterior of $\mathbb{T}$ and extends
continuously to $\mathbb{T}$. Note that $G$ maps $\mathbb{T}$ to $\mathbb{T}$,
by Schwarz reflection principle, $G$ can be continued to be a meromorphic
function in $\mathbb{C}^{*}$. Let us still denote the map by $G$. By the
construction, it is not difficult to see that $G(\bar{z})=\overline{G(z)}$,
$G(-1)=1$ and $G(1)=1$. This implies that $G$ is post-critically finite on
$\mathbb{T}$. Since $\psi$ is normalized and with the quasiconformal constant
being uniformly bounded by the main lemma, $G$ is holomorphic in an annular
neighborhood of $\mathbb{T}$ with modulus bounded away from zero. This implies
that $G|\mathbb{T}$ belongs to a compact family of analytic circle mappings.
By Theorem 6.2 of Chapter 3 in [7], there is a quasi-symmetric circle
homeomorphism $h:\mathbb{T}\to\mathbb{T}$ such that
$G|\mathbb{T}=h^{-1}\circ\Lambda\circ h$ with $\Lambda$ being the square map.
It is not difficult to see that $h$ is also symmetric about the real axis in
the sense $\overline{h(e^{it})}=h(e^{-it})$. In particular $h(1)=1$ and
$h(-1)=-1$. Since $G|\mathbb{T}$ belongs to a compact family of analytic
circle mappings, the quasi-symmetric constant of $h$ is bounded by some
uniform constant $1<M<\infty$.

######  Lemma 2.1.

The quasi-symmetric circle mapping $h$ can be extended to a quasiconformal
homeomorphism $H:\mathbb{D}\to\mathbb{D}$ so that

  * •

$H(0)=0$ and $H$ is symmetric about the real axis, that is,
$H(\bar{z})=\overline{H(z)}$,

  * •

the quasiconformal constant of $H$ is bounded by some constant depending only
on $M$.

###### Proof.

We first define a homeomorphism $h:[-1,1]\to[-1,1]$. To do that, note that $G$
is symmetric about the real axis and thus $G^{\prime}(1)=\lambda>1$ is a real
number. For $k\geq 0$, let $x_{k}=1-\frac{1}{{\lambda}^{k}}$ and
$y_{k}=1-\frac{1}{{2}^{k}}$. Define $h(x_{k})=y_{k}$ for all $k\geq 0$, and
define $h$ on each $[x_{k},x_{k+1}]$ so that $h$ maps $[x_{k},x_{k+1}]$
linearly onto $[y_{k},y_{k+1}]$. For $-1\leq x\leq 0$, define $h(x)=-h(-x)$.
In this way we have defined a homeomorphism $h:\partial U\to\partial U$ so
that $h(0)=0$ where $U$ is the upper half unit disk. Let $\xi:U\to\mathbb{D}$
be the Riemann isomorphism which fixes $-1$ and $1$ and maps $0$ to $-i$. Then
$\xi$ behaves like the square map near $-1$ and $1$ because it opens the right
angles to straight angles at the two points. From the definition of $h$ on
$[-1,1]$, especially its behavior near the two end points $1$ and $-1$, it is
not difficult to check that $\xi\circ h\circ\xi^{-1}:\mathbb{T}\to\mathbb{T}$
is quasi-symmetric with the quasi-symmetric constant depending only on that of
$h$. Hence $\xi\circ h\circ\xi^{-1}$ can be extended to a quasiconformal
homeomorphism $\Pi:\mathbb{D}\to\mathbb{D}$ with the quasiconformal constant
bounded by some uniform constant. Define $H:U\to U$ by setting
$H(z)=\xi^{-1}\circ\Pi\circ\xi$. On the lower half unit disk, we define $H$ by
setting $H(\bar{z})=\overline{H(z)}$. This completes the proof of the lemma. ∎

Define

(2.1) $\widehat{G}(z)=\begin{cases}G(z)&\text{for
$z\in{\mathbb{C}}-\mathbb{D}$},\\\ H^{-1}\circ\Lambda\circ H(z)&\text{ for
$z\in\mathbb{D}$}.\end{cases}$

Let $\nu_{0}$ be the complex structure in $\mathbb{D}$ given by the pull back
of the standard one by $H$. Since $H$ is symmetric about the real axis,
$\nu_{0}$ is symmetric about the real axis. We then pull back $\nu_{0}$ to the
complex plane by the iteration of $G$ and get a $\widehat{G}$-invariant
complex structure $\nu$. Thus $\nu$ is symmetric about the real axis. Since
all the maps involved are uniformly quasiconformal, we have
$\|\nu\|_{\infty}<\kappa$ for some uniform $0<\kappa<1$. Let $\Phi$ be the
quasiconformal homeomorphism of the plane which fixes $0$ and $1$ and solves
the Beltrami equation given by $\nu$. Since $\nu$ is symmetric about the real
axis, so is $\Phi$. Let

$f=\Phi\circ\widehat{G}\circ\Phi^{-1}.$

Then $f$ is an entire function.

######  Lemma 2.2.

The following assertions hold.

  * •

$f$ has no asymptotic values,

  * •

$f$ is symmetric about the real axis, and fixes $0$ and $1$,

  * •

$f$ has a critical point at $0$ and all the other critical points are negative
reals and are mapped either to $0$ or to $1$,

  * •

all the Fatou components of $f$ are quasi-disks and are eventually mapped to
the super-attracting component at the origin.

###### Proof.

Suppose $f$ has an asymptotic value. From the construction of $f$,
$\widehat{G}$ and thus $F$ would have an asymptotic value. This contradicts
the main lemma. This proves the first assertion.

Since $\widehat{G}$ and $\Phi$ are both symmetric about the real axis, so is
$f$. Since $\Phi$ and $\widehat{G}$ fix $0$ and $1$, $f$ fix $0$ and $1$ also.
This proves the second assertion.

By the symmetry of $\Phi$, $\Phi$ maps the real axis to the real axis. Since
$\Phi$ fixes $0$ and $1$, $\Phi$ maps the negative reals to the negative
reals. From the construction of $f$, the critical set of $f$ are the
$\Phi$-image of that of $\widehat{G}$. The critical points of $\widehat{G}$,
except the one at the origin, are exactly those of $G$ which belong to
$(-\infty,-1]$. By the construction of $G$, these are the
$(\eta^{-1}\circ\psi)$-image of the critical points of $F$ which belongs to
$(-\infty,-1]$. By the symmetry and the normalization conditions of $\psi$ and
$\eta$, $\eta^{-1}\circ\psi$ maps $(-\infty,-1]$ homeomorphically to
$(-\infty,-1]$. This, together with the fact that $\Phi$ maps $(-\infty,0]$
homeomorphically to $(-\infty,0]$, implies that all the other critical points
of $f$ are negative reals. Since $F$ maps any of its critical points either to
$1$ or to $0$, so does $\widehat{G}$ by the definition of $\widehat{G}$. Since
$\Phi$ fixes $0$ and $1$, it follows that $f$ maps any of its critical points
either to $0$ or to $1$. This proves the third assertion.

As we have seen, $f$ has two singularly values, $0$ and $1$, both of which are
critical values. By [3] and [5], $f$ has neither Baker domains nor wandering
domains. Since entire functions have no Herman rings, every Fatou component of
$f$ must be eventually mapped into some periodic cycle of either Siegel disks,
or parabolic components, or attracting components, or supper-attracting
components. Note that both the two singular values of $f$ are fixed. Since
each point of the boundary of a Siegel disk must be accumulated by the forward
orbit of the singular values, $f$ has no Siegel disks. Since any parabolic
basin or attracting basin must contain an infinite forward orbit of some
singular value, $f$ have not periodic parabolic or attracting cycles. Thus all
the periodic cycles of Fatou components of $f$ must be super-attracting. Let
$\Omega_{0}$ be the Fatou component containing the origin. Since the origin is
the only periodic (fixed) critical point of $f$, all the Fatou components must
be eventually mapped to $\Omega_{0}$. Since $\widehat{G}$ is conjugate to the
square map on $\overline{\mathbb{D}}$ and since the iteration of the square
map is normal in $\mathbb{D}$ but not normal at any point in
$\partial\mathbb{D}$, we must have $\Omega_{0}=\Phi(\mathbb{D})$. In
particular, $\Omega_{0}$ is a quasi-disk. By the construction of $f$, any
component of $f^{-1}(\Omega_{0})$ other than $\Omega_{0}$ is the
$(\Phi\circ\eta^{-1}\circ\psi)$-image of some $P_{n}$ with $n\geq 1$, and is
thus a quasi-disk. Since the closure of each
$\Phi\circ\eta^{-1}\circ\psi(P_{n}),n\geq 1$, does not intersect the forward
orbit of the singular values, its pre-images under the iteration of $f$ must
be quasi-disks also. So all the Fatou components of $f$ are quasi-disks. This
proves the last assertion of the lemma. ∎

Now let us prove a uniform contraction property of $G$. Let ${\rm
diam}(\cdot)$ denote the diameter with respect to the Euclidean metric.

######  Lemma 2.3.

For any $L>1$ there is an $n\geq 1$ depending only on $L$ and a curve segment
$\Gamma$ attached to $\mathbb{T}$ from the outside such that

  * •

${\rm diam}(\Gamma)<1$ and

  * •

$G^{n}:\Gamma\to[-L,-1]$ is a homeomorphism.

###### Proof.

Let $h$ be the quasi-symmetric circle homeomorphism which conjugates
$G|\mathbb{T}$ to the square map. For $k\geq 1$ let
$z_{k}=h^{-1}(e^{i\pi/2^{k}})$. Since the quasi-symmetric constant of $h$ is
bounded by some uniform $1<M<\infty$, it follows that $z_{k}\to 1$ uniformly
in the following sense: for any $\epsilon>0$, there is an $N\geq 1$
independent of $G$ so that $|z_{k}-1|<\epsilon$ for all $k\geq N$. Let
$\Gamma_{k}$ be the curve segment attached to $z_{k}$ from the outside of
$\mathbb{T}$ so that $G^{k}:\Gamma_{k}\to[-L,-1]$ is a homeomorphism. Let
$X=\mathbb{C}\setminus\\{0,1\\}$ and $Y=G^{-1}(X)$. Then $G:Y\to X$ is a
holomorphic covering map. From $G(-1)=1$ it follows that $Y\subset
X\setminus\\{-1\\}$. Thus the pull back by $G$ decreases the hyperbolic metric
in $X$. Let $l_{X}(\cdot)$ denote the length with respect to the hyperbolic
metric in $X$. Then $l_{X}(\Gamma_{k})\leq l_{X}([-L,-1])$ for all $k\geq 1$.
This, together with that $z_{k}\to 1$ uniformly, implies the existence of the
integer $n$ such that ${\rm diam}(\Gamma_{n})<1$. ∎

To construct the entire function $f$ with non-locally connected Julia sets, we
need to take appropriate squares $\\{P_{n}\\}_{n\geq 1}$ in the construction
of the quasi-regular map $F$. As we have seen above, $\mathbb{D}$ corresponds
to a quasi-disk $\Omega_{0}$ which is fixed by $f$, and each $P_{n}$
corresponds to a quasi-disk $\Omega_{n}$ which is mapped to $\Omega_{0}$ by
$f$. More precisely,

$\partial\Omega_{n}=\Phi\circ\eta^{-1}\circ\psi(\partial P_{n})\>\>\>\>\hbox{
and }\partial\Omega_{0}=\Phi(\mathbb{T}).$

By symmetry, all $\Omega_{n}$ are symmetric about the real axis. From the
construction, it also follows that $\\{-1,1\\}\subset\partial\Omega_{0}$ and
$f(-1)=f(1)=1$. Let $\partial\Omega_{l}\cap\mathbb{R}=\\{a_{l},b_{l}\\}$ with
$a_{l}>b_{l}$. Since $\Phi$ is uniformly quasiconformal and normalized, it
follows that there exists some $R>0$ independent of the choice of
$\\{P_{n}\\}_{n\geq 1}$ such that

$\overline{\Omega_{0}}\subset B_{R}(0)$

where $B_{R}(0)$ is the disk centered at the origin and with radius $R$. From
the fact that $\psi(-1)=\eta(-1)=\Phi(-1)=-1$ it follows that $a_{1}=-1$.
Since $\Phi$, $\eta$ and $\psi$ are uniformly quasiconformal and normalized,
we can take $P_{1}$ large enough so that $|b_{1}|>R+1$ holds for any choice of
the remaining squares $\\{P_{n}\\}_{n\geq 2}$. Let $n_{1}=0$. Now for $k\geq
1$, suppose the diameters of all $P_{l}$, $1\leq l\leq k$, have been
determined such that for each $1\leq l\leq k$, there is an integer $n_{l}$ and
a branch of $f^{-n_{1}}$, say $\chi_{l}$, such that $|\chi_{l}(a_{l})|<R$ and
$|\chi_{l}(b_{l})|>R+1$ hold for any choice of $\\{P_{n}\\}_{n\geq l+1}$. Let
us determine the diameter of $P_{k+1}$ and $n_{k+1}$ as follows. Since $P_{k}$
has determined and since $\Phi$, $\eta$ and $\psi$ are uniformly
quasiconformal and normalized, it follows that $a_{k+1}=b_{k}$ is bounded and
the bound does not depend on the choice of $\\{P_{n}\\}_{n\geq k+1}$. By Lemma
2.3 and the fact that $\Phi$ is uniformly quasiconformal and normalized, we
can choose $n_{k+1}$ independent of the choice of $\\{P_{n}\\}_{n\geq k+1}$,
such that the corresponding branch of $f^{-n_{k+1}}$ , say $\chi_{k+1}$,
satisfying $|\chi_{k+1}(a_{k+1})|<R$. Since $F$ is uniformly bounded on
compact sets by the main lemma, and since $\Phi$, $\eta$ and $\psi$ are
uniformly quasiconformal and normalized, it follows that $f$ is uniformly
bounded on compact sets. Thus for the $n_{k+1}$, we can choose $P_{k+1}$ large
enough to ensure that $|\chi_{k+1}(b_{k+1})|>R+1$ holds for any choice of
$\\{P_{n}\\}_{n\geq k+2}$.

In this way we get a sequence of squares $\\{P_{n}\\}$. For such
$\\{P_{n}\\}$, let $F$ be the quasi-regular map in the main lemma and $f$ be
the entire function obtained as above.

######  Lemma 2.4.

There exist an $\epsilon>0$ and an infinite sequence of Fatou components of
$f$ such that the spherical diameter of each Fatou component in this sequence
is greater than $\epsilon$.

###### Proof.

From the construction of $P_{n}$, it follows that there exist a constant $R>0$
and an infinite sequence of Fatou components of $f$ which intersect both of
the two Euclidean circles $\\{z:\>|z|=R\\}$ and $\\{z:\>|z|=R+1\\}$. It is
clear that the spherical diameter of these Fatou components have a uniform
positive lower bound depending only on $R$. ∎

Before we show that $J(f)$ is not locally connected let us recall a criterion
of the local connectedness.

######  Theorem 2.5 (Whyburn, [10]).

A closed set in $\widehat{\mathbb{C}}$ is locally connected if and only if it
satisfies the following two conditions:

  1. (a)

The boundary of each complementary component of it is locally connected;

  2. (b)

For an arbitrary $\epsilon>0$, the number of the complementary components
whose diameters with respect to the spherical metric exceed $\epsilon$ is
finite.

######  Lemma 2.6.

$J(f)$ is not locally connected.

###### Proof.

By Lemma 2.2 all the Fatou components of $f$ are bounded and simply connected.
By Theorem 2 in [6] $J(f)$ is connected. Thus $J(f)\cup\\{\infty\\}$ is a
continuum in $\widehat{\mathbb{C}}$. By Corollary 5.13 in [8] a continuum in
$\widehat{\mathbb{C}}$ can not fail to be locally connected at a single point.
It follows that $J(f)$ is locally connected if and only if
$J(f)\cup\\{\infty\\}$ is locally connected. By Lemma 2.4 and Theorem 2.5 it
follows that $J(f)\cup\\{\infty\\}$ is not locally connected. Thus $J(f)$ is
not locally connected. This proves Lemma 2.6. ∎

Theorem 1.1 now follows from Lemmas 2.2 and 2.6.

##  3\. Proof of the main lemma

Suppose $\\{P_{n}\\}$ is a sequence of squares in the Main Lemma. See Figure 2
for an illustration. Here the diameters of $P_{1}$, $P_{2}$ , $P_{3}$ and
$P_{4}$ are $2$, $4$, $2$ and $8$ respectiely. We divide the plane into
disjoint union of “half strips” and squares. Up to symmetry, there are five
types of half strips. Except the three special ones whose shapes are
independent of the choice of $\\{P_{n}\\}_{n\geq 1}$ (see Figures 3-5), each
of the other ones has the following common properties.

  * •

The two bottom angles are $\pi/4$ and $3\pi/4$,

  * •

the width is equal to one,

  * •

the length of the bottom is equal to $\sqrt{2}$,

  * •

it is bent to a right angle at some place.

Figure 2. Constructing the quasi-regular map $F$

We will construct $F$ by defining it to be $K$-quasiconformal on each of these
half strips and squares with $K>1$ being some constant independent of the
choice of $\\{P_{n}\\}$. We then glue these maps together along the boundary
of their domains so that the global map satisfies the properties in the lemma.
First we will define a quasiconformal map $g$ on each of the half strip so
that it maps each half strip onto the standard half strip
$S=\\{x+iy\>|\>-1<x<1,y>0\\}$. In this process the situations for types II and
III are similar with the type I and the situation for type V is similar with
the type IV. So we will only present the details of the construction for the
half strips of types I and IV. It should be easy for the reader to supply the
details for the other situations.

Let $H$ be the half strip of the type I. Let $AB$ denote the $1/4$-unit circle
and $A^{\prime}B^{\prime}$ denote the bottom side of $S$. See Figure 3 for an
illustration. Define $g$ on $AB$ so that it maps $AB$ proportionally onto
$A^{\prime}B^{\prime}$. Let $L$ and $R$ be the two components of $\partial
H\setminus AB$ so that $L$ is the one with $B$ being the end point and $R$ is
the one with $A$ being the end point. Similarly, let $L^{\prime}$ and
$R^{\prime}$ be the two corresponding components of $\partial S\setminus
A^{\prime}B^{\prime}$. For any $X\in L$, define $g(X)=X^{\prime}\in
L^{\prime}$ such that the subarc $BX\subset L$ has the same length as the
subarc $B^{\prime}X^{\prime}\subset L^{\prime}$, and for $Y\in R$, define
$g(Y)=Y^{\prime}\in R^{\prime}$ such that the subarc $AY\subset R$ has the
same length as the subarc $A^{\prime}Y^{\prime}\subset R^{\prime}$. In this
way we have defined $g:\partial H\to\partial S$ so that it preserves the
length on the two boundary sides $L$ and $R$. To define $g$ in the interior of
$H$, let $C$ be the vertex of $H$ at the upper left corner and take $D\in L$
and $E\in R$ such that $|CD|=2$ and $|AE|=1$. Let $C^{\prime}=g(C)$,
$D^{\prime}=g(D)$ and $E^{\prime}=g(E)$. Define $g$ on $DE$ so that it maps
$DE$ linearly onto $D^{\prime}E^{\prime}$. Now it is not difficult to extend
$g$ to be a quasiconformal homeomorphism from the domain $ABCDE$ to the
pentagon $A^{\prime}B^{\prime}C^{\prime}D^{\prime}E^{\prime}$. For instance,
one can divide the two domains into the same number of star domains and then
apply Theorem 6 of [9]. Let $H_{0}$ and $S_{0}$ denote the remaining part of
$H$ and $S$ respectively. Note that $g$ has been defined on $\partial H_{0}$
and is linear on each straight boundary piece of $H_{0}$. We can thus extend
$g$ to be a linear homeomorphism from $H_{0}$ to $S_{0}$. It is clear the map
$g:H\to S$ obtained in this way is a quasiconformal homeomorphism, and most
importantly, it preserve the length on the two boundary sides of $H$. The
similar construction can be applied to the half strips of types II and III.
See Figures 4 and 5 for an illustration. We leave the details to the reader.

Figure 3. Define $g$ in the half strip of type I Figure 4. Defining $g$ in the
half strip of type II Figure 5. Define $g$ in the half strip of type III

Now let us define $g$ in a half strip $H$ of type IV. See Figure 6 for an
illustration. As before we first define $g$ on the bottom side $AB$ of $H$ to
be linear on $AB$ so that $g$ maps $AB$ to $A^{\prime}B^{\prime}$. Then let
$L$ and $R$ denote respectively the two components of $\partial H\setminus
AB$, one with $B$ being the end point and the other one with $A$ being the end
point. Similarly, Let $L^{\prime}$ and $R^{\prime}$ denote the corresponding
components of $\partial S\setminus A^{\prime}B^{\prime}$. For $X\in L$, define
$g(X)=X^{\prime}\in L^{\prime}$ so that the subarc $BX\subset L$ has the same
length as the subarc $B^{\prime}X^{\prime}\subset L^{\prime}$. Similarly, for
$Y\in R$, define $g(Y)=Y^{\prime}\in R^{\prime}$ so that the subarc $AY\subset
R$ has the same length as the subarc $A^{\prime}Y^{\prime}\subset R^{\prime}$.
In this way we have defined $g:\partial H\to\partial S$.

Let $D$ and $E$ denote the two corner points in $\partial H$ as indicated in
Figure 6. To define $g$ in the interior of $H$, we take a points $C\in L$ such
that $|BC|=|AE|$. Let $C^{\prime},D^{\prime},E^{\prime}\in\partial S$ so that
$C^{\prime}=g(C)$, $D^{\prime}=g(D)$ and $E^{\prime}=g(E)$. Then
$|CD|=|C^{\prime}D^{\prime}|$, $|BC|=|B^{\prime}C^{\prime}|$,
$|AE|=|A^{\prime}E^{\prime}|$ and $g$ maps $CD$, $BC$ and $AE$ linearly onto
$C^{\prime}D^{\prime}$, $B^{\prime}C^{\prime}$ and $A^{\prime}E^{\prime}$,
respectively. Define $g$ on $CE$ and $DE$ so that $g$ maps them linearly onto
$C^{\prime}E^{\prime}$ and $D^{\prime}E^{\prime}$ respectively. Then $g$ can
be extended linearly to the interiors of the parallelogram $ABCE$ and the
triangle $CED$ so that $g$ maps them onto the interiors of the rectangle
$A^{\prime}B^{\prime}C^{\prime}E^{\prime}$ and the triangle
$C^{\prime}E^{\prime}D^{\prime}$ respectively. The remaining part of $H$, say
$H_{0}$, is a half strip with $DE$ being the bottom side. Correspondingly, the
remaining part of $S$, say $S_{0}$, is also a half strip with
$D^{\prime}E^{\prime}$ being the bottom side. Note that $g$ has been defined
on $\partial H_{0}$ so that $g:\partial H_{0}\to\partial S_{0}$ is linear on
each straight boundary piece. Thus $g$ can be extended to the interior of
$H_{0}$ so that $g$ maps $H_{0}$ linearly onto $S_{0}$. Now we define $g:H\to
S$ by gluing the maps on the three pieces along the boundary. It is clear that
$g$ is a quasiconformal map. Since the quasiconformal constant of $g$ on each
piece is some uniform constant, thus the quasiconformal constant of $g$ on $H$
is a uniform constant. Moreover, from the construction we see that $g$
preserve the length on the two boundary side of $\partial H$.

Figure 6. Define $g$ in the half strip of type IV Figure 7. Define $g$ in the
half strip of type V

For a half strip $H$ of type V, the method above can be applied in a similar
way to get a $K$-quasiconformal homeomorphism $g:H\to S$ such that (1) $K>1$
is some uniform constant, and (2) $g$ is linear on the bottom side and
preserves the length on the two boundary sides of $H$. See Figure 7 for an
illustration. We leave the details to the reader.

Now let us define $F$ on the half strips. Let $H_{+}$ and $H_{-}$ denote
respectively the upper and lower components of
$\mathbb{C}\setminus(-\infty,-1]\cup\overline{\mathbb{D}}\cup[1,+\infty)$.
Define $\sigma:S\to H_{+}$ by $\sigma(z)=e^{i\pi(1-z)/2}$ and $\tau:S\to
H_{-}$ by $\tau(z)=-e^{i\pi(1-z)/2}$. It remains to verify that whenever a
point belong to the intersection of the boundaries of two half strips, the two
maps defined on the two half strips must coincide on this point. To see this,
we may assume that $z$ belong to the $L$ part of the boundary of some red half
strip (that is, marked with red rectangle) and belongs to the $R$ part of the
boundary of some blue half strip (that is, marked with blue rectangle). Let
$g_{r}$, $F_{r}$, and $g_{b}$, $F_{b}$ denote respectively the maps $g$ and
$F$ defined on the red and blue half strips. Since both $g_{r}$ and $g_{b}$
preserve the length on the two boundary side of the half strips, $g_{r}(z)\in
L^{\prime}$ and $g_{b}(z)\in R^{\prime}$ have the same imaginary part. By the
definition of $\sigma$, it follows that

$F_{r}(z)=\sigma(g_{r}(z))=-\sigma(g_{b}(z))=\tau(g_{b}(z))=F_{b}(z).$

Figure 8. $\theta\to\xi(\theta)$,
$r\to\frac{r^{d}}{R(\theta)^{d}}e^{\xi(\theta)}$

It remains to define $F$ in the interior of the squares. Suppose that the
diameter of $P_{n}$ is $d=2m$ for some integer $m\geq 1$. Note that we have
defined $F$ on the boundary of $P_{n}$ in the above. From the construction,
$F$ maps $\partial P_{n}$ proportionally onto $\mathbb{T}$ and $F:\partial
P_{n}\to\mathbb{T}$ is a covering of degree $d$. Suppose $\omega$ is the
center of $P_{n}$. For $0\leq\theta\leq 2\pi$, let $R(\theta)>0$ and
$0\leq\xi(\theta)\leq 2\pi$ be the numbers such that
$\omega+R(\theta)e^{i\theta}\in\partial P_{n}$ and

$F(\omega+R(\theta)e^{i\theta})=e^{i\xi(\theta)}.$

Since $F$ maps $\partial P_{n}$ proportionally onto $\mathbb{T}$, it follows
that $\xi^{\prime}(\theta)$ exists except the fours angles corresponding to
the four vertices of $P_{n}$, and moreover,
$1/C<\xi^{\prime}(\theta_{1})/\xi^{\prime}(\theta_{2})<C$ for some uniform
$C>1$. Since $F:\partial P_{n}\to\mathbb{T}$ is a covering of degree $d$ we
get $\xi^{\prime}(\theta)/d>0$ and is uniformly bounded away from zero and
$1$. For $0\leq r\leq R(\theta)$, define

$F(\omega+re^{i\theta})=\frac{r^{d}}{R(\theta)^{d}}e^{i\xi(\theta)}.$

We claim that $F$ is $K$-quasi-regular for some universal $K>1$ in $P_{n}$.
For the convenience of the computation we adopt the log coordinate. Consider
the map

$\zeta\to\log F(\omega+e^{\zeta})=u+iv$

where $\zeta=x+iy$. Then $x=\log r$ and $y=\theta$. We thus get

$u(x,y)=dx-d\log R(y)\hbox{ and }v(x,y)=\xi(y).$

Let $h=u+iv$. Let us compute $h_{\bar{z}}$ and $h_{z}$. For $A,B>0$ we denote
$A=O^{*}(B)$ if $A/B$ is uniformly bounded away from zero and $1$. We use
$O(1)$ to denote a number whose absolute value is bounded by some universal
constant. By a simple computation we have

$h_{x}=d\hbox{ and }h_{y}=-d\frac{R^{\prime}(y)}{R(y)}+i\xi^{\prime}(y).$

Note that $R^{\prime}(y)=O(1)$, $R(y)=O^{*}(d)$ and
$\xi^{\prime}(y)=O^{*}(d)$. From $h_{z}=\frac{1}{2}(h_{x}-ih_{y})$ and
$h_{\bar{z}}=\frac{1}{2}(h_{x}+ih_{y})$ we get

$h_{\bar{z}}=\frac{1}{2}(d-O^{*}(d))+i\lambda(z)$

and

$h_{z}=\frac{1}{2}(d+O^{*}(d))-i\lambda(z)$

where $\lambda(z)=O(1)$. This implies that $|h_{\bar{z}}/h_{z}|<k<1$ for some
uniform $0<k<1$ independent of $d$.

Now from the construction the first four assertions of the main lemma hold. To
prove the last one, we first take a smallest closed disk centered at the
origin which contains the compact set. Note that the bound of $F$ on the
closed disk is essentially determined by the length of the longest half strip
contained in it. But this is then bounded above by some constant depending
only on the diameter of the disk. Since the diameter of the disk depends on
the compact set, the last assertion follows. The proof of the Main Lemma is
completed.

$\mathbb{Acknowledgements.}$ We would like to thank D. Sixsmith and L. Rempe
for their detailed comments on an early version of this work. Further thanks
are due to the anonymous referees for their suggestions which greatly improve
the paper. The second author is partially supported by NSFC(11171144,
11325104).

## References

  * [1] Bergweiler W, Fagella N, Rempe-Gillen L. Hyperbolic entire functions with bounded Fatou components. Comment Math Helv(4), 2015, 90: 799–829
  * [2] Bishop C. Constructing entire functions by quasiconformal folding. Acta Math(1), 2015, 214: 1–60
  * [3] Eremenko A, Lyubich M. Dynamical properties of some classes of entire functions. Ann Inst Fourier(4), 1992, 42: 989–1020
  * [4] Eremenko A, Sodin M. Parametrization of entire functions of sinetype by their critical values, Entire and subharmonic functions. Adv Soviet Math, vol. 11, Amer Math Soc, Providence, RI, 1992, 237–242
  * [5] Goldberg L, Keen L. A finiteness theorem for a dynamical class of entire functions. Ergd Th & Dyna Sys(2),1986, 6:183–192
  * [6] Kisaka M. On the connectivity of Julia sets of transcendental entire functions. Ergd Th & Dyna Sys(1), 1998,18:189–205
  * [7] De Melo W, Strien S. One-Dimensional Dynamics. Springer-Verlag Berlin Heidelberg, 1993
  * [8] Sam B. Nadler, Jr. Continuum theory. An introduction. Monograph and Textbooks in Pure and Applied Mathematics, vol.158, Marcel Dekker Inc., New York, 1992
  * [9] Nicks D A, Sixsmith D J. Periodic domains of quasiregular maps. Mathematics(4), 2017, 141: 1385–1392
  * [10] Whyburn G T. Analytic Topology. American Mathematical Society Colloquium Publications, 28, American Mathematical Society, New York, 1942


# Orthogonal Time Frequency Space Modulation – Part III: ISAC and Potential
Applications

Weijie Yuan, , Zhiqiang Wei, , Shuangyang Li, ,
Robert Schober, , and Giuseppe Caire,
(_Invited Paper_) W. Yuan is with the Department of Electrical and Electronic
Engineering, Southern University of Science and Technology, Shenzhen 518055,
China (e-mail: yuanwj@sustech.edu.cn). Z. Wei is with the School of
Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, China
(e-mail: zhiqiang.wei@xjtu.edu.au). S. Li was with the School of Electrical
Engineering and Telecommunications, University of New South Wales, Sydney, NSW
2052, Australia, when this letter was submitted. He is now with the Department
of Electrical, Electronic, and Computer Engineering, University of Western
Australia, Perth, WA 6009, Australia (e-mail: shuangyang.li@uwa.edu.au). R.
Schober is with the Institute for Digital Communications (IDC), the Friedrich-
Alexander University Erlangen-Nuremberg, Erlangen 91054, Germany (e-mail:
robert.schober@fau.de). G. Caire is with the Department of Electrical
Engineering and Computer Science, Technical University of Berlin, Berlin
10587, Germany (caire@tu-berlin.de).

###### Abstract

The first two parts of this tutorial on orthogonal time frequency space (OTFS)
modulation have discussed the fundamentals of delay-Doppler (DD) domain
communications as well as some advanced technologies for transceiver design.
In this letter, we will present an OTFS-based integrated sensing and
communications (ISAC) system, which is regarded as an enabling technology in
next generation wireless communications. In particular, we illustrate the
sensing as well as the communication models for OTFS-ISAC systems. Next, we
show that benefiting from time-invariant DD channels, the sensing parameters
can be used for inferring the communication channels, leading to an efficient
transmission scheme. As both functionalities are realized in the same DD
domain, we briefly discuss several promising benefits of OTFS-based ISAC
systems, which have not been completely unveiled yet. Finally, a range of
potential applications of OTFS for the future wireless networks will be
highlighted.

###### Index Terms:

OTFS, integrated sensing and communications (ISAC), delay-Doppler domain

##  I Introduction

The orthogonal time frequency space (OTFS) modulation has gained considerable
attentions due to its capability of providing more reliable communications
compared to the existing orthogonal frequency division multiplexing (OFDM)
modulation, especially in high-mobility environments [1]. As overviewed in
Part I of this tutorial, OTFS utilizes the delay-Doppler (DD) domain for data
multiplexing, which mirrors the geometry of the scatterings comprising the
wireless channel. In history, the delay and Doppler parameters are typically
used for sensing-related applications, representing the range and velocity
characteristics corresponding to the moving targets. In future wireless
networks, both sensing and communications functionalities are highly desired,
leading to the integrated sensing and communications (ISAC) technology by
using the same signaling waveform and hardware architecture [2, 3]. In
contrast to the OFDM-based ISAC system, OTFS-ISAC system provides the direct
interaction between transmitted signals and the channel responses in a unified
DD domain for both functionalities, on top of the resilience to delay and
Doppler spreads. To this end, OTFS has been regarded as a promising signal
waveform to fully unleash the potentials of ISAC system.

The pioneering work [4] compared the radar sensing performance based on OTFS
and OFDM schemes, showing that OTFS is more robust for velocity estimation.
Some recent works discussed the ISAC performance limits, the beamforming
designs, the resource allocation schemes, and the receiver structures for
OTFS-based ISAC systems [5, 6, 7, 8, 9]. All these contributions have
demonstrated the effectiveness of OTFS for ISAC systems. However, considering
the inherent relationship between DD domain communications and sensing, there
are several promising benefits of adopting OTFS as the waveform for realizing
ISAC, which have been not well discussed in the literature.

In this letter, we will first study the OTFS signaling model for both
communications and sensing in the DD domain. It is noted that the delay and
Doppler associated with the DD domian communication channel are half of the
round-trip delay and Doppler parameters extracted from the sensing echoes.
This motivates us that the sensing parameters can be used for predicting the
communication channels given the slow time-varying property of the DD domain
channel. Afterwards, we briefly discuss several benefits of OTFS-ISAC system,
i.e., unified signaling design framework, sensing channel exploitation, and
communications-assisted sensing. Finally, we discuss potential applications
related to OTFS modulation, including underwater communications, optical
transmission, and space-air-ground integrated networks.

_Notations:_ $\left(\cdot\right)^{\rm H}$ and $\left(\cdot\right)^{\rm T}$
denotes the Hermittian operation and the transpose operation, respectively;
$\left(\cdot\right)_{N}$ denotes the modulus operation w.r.t. $N$;
$\mathbf{I}_{MN}$ denotes an identity matrix of dimension $MN$;
$\delta\left(\cdot\right)$ denotes the Dirac delta function;
$\mathbb{E}[\cdot]$, $\textrm{Tr}[\cdot]$, and $\textrm{det}[\cdot]$ denote
the expectation operation, trace operation, and the determinant of a matrix,
respectively.

##  II OTFS-ISAC System Model

Let us consider a general system model with a mono-static base station (BS)
and $P$ targets of interests, as shown in Fig. 1. For ease of exposition, we
assume that all targets of interest are communications user equipments (UEs)
as well, which aligns well with the use case such as vehicular networks. The
BS is equipped with a transmit uniform linear array (ULA) of $N_{t}$ antennas
and a separate receive ULA of $N_{r}$ antennas. Under the assumption of
sufficient isolation between the transmit and receive arrays, the sensing
echoes would not interfere with the transmitted signals. For receiving the
communication information, each target is equipped with $N_{u}$ antennas.

Figure 1: OTFS-ISAC system model.

###  II-A Communication Model

Without loss of generality, we consider a DD domain grid of $M\times N$, where
$l\in\\{0,\ldots,M-1\\}$ and $k\in\\{0,\ldots,N-1\\}$ denotes the delay and
Doppler indices, respectively. The information symbols sent to the $i$-th
target are defined as $X_{\rm DD}^{i}\left[l,k\right]$, which is transformed
to the time-frequency (TF) domain via inverse symplectic finite Fourier
transform (ISFFT), yielding $X^{i}_{\rm TF}\left[m,n\right]$, where $m$ and
$n$ denote the indices of subcarriers and time slots. By performing
multicarrier modulation, the transmitted signal to the $i$-th target is given
by

$\displaystyle s_{i}(t)=\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}X^{i}_{\rm
TF}\left[m,n\right]g_{\rm tx}(t-nT)e^{j2\pi m\Delta f(t-nT)},$ (1)

where $g_{\rm tx}(t)$ is the transmit shaping pulse, $T$ and $\Delta f$ denote
the symbol duration and frequency spacing, respectively, satisfying $T\Delta
f=1$ to maintain the orthogonality.

Initially at the target detection stage, the BS could use the antenna array to
formulate a wide beam or an omni-directional signal to detect all targets of
interests. Afterwards, in the tracking mode, the BS formulates multiple beams
for information transmission and target tracking. Through a beamforming matrix
$\mathbf{F}\in\mathbb{C}^{N_{t}\times P}$, the multi-beam signal is expressed
as

$\displaystyle\tilde{\mathbf{s}}(t)=\mathbf{F}\mathbf{s}(t),$ (2)

where $\mathbf{s}(t)=\left[s_{1}(t),...,s_{P}(t)\right]^{\rm T}$ is the
transmitted signal to all $P$ targets. The $i$-th column of $\mathbf{F}$,
denoted by
$\mathbf{f}_{i}=\sqrt{\frac{p_{i}}{N_{t}}}\mathbf{a}_{N_{t}}(\tilde{\theta}_{i})$
is used to allocate power $p_{i}$ and to steer the transmitted signal towards
the intended direction $\tilde{\theta}_{i}$, where the steering vector is
given by

$\displaystyle\mathbf{a}_{N_{t}}(\theta_{i})=\left[1,e^{j\pi\sin\theta_{i}},...,e^{j(N_{t}-1)\pi\sin\theta_{i}}\right]^{\rm
T}.$ (3)

Due to the asymptotic orthogonality of the massive antenna array, after
transmit beamforming, the communication channel to the $i$-th target is
asymptotically line-of-sight (LoS) dominated, and given by

$\displaystyle\mathbf{C}^{i}_{\rm
DD}(\tau,\nu)=h_{i}\mathbf{a}_{N_{u}}(\theta_{i})\mathbf{a}_{N_{t}}^{\rm
H}(\theta_{i})\delta(\tau-\tau_{i})\delta(\nu-\nu_{i}),$ (4)

with channel gain $h_{i}=\sqrt{\frac{c}{4\pi f_{c}d^{2}_{i}}}$, signal
propagation speed $c$, carrier frequency $f_{c}$, range distance $d_{i}$,
delay $\tau_{i}$, and Doppler shift $\nu_{i}$. After passing through the
channel [10], the received signal of the $i$-th target can be modelled as

$\displaystyle{y}_{i}(t)=h_{i}\mathbf{u}^{\rm
H}_{i}\mathbf{a}_{N_{u}}(\theta_{i})\mathbf{a}_{N_{t}}^{\rm
H}(\theta_{i})\mathbf{f}_{i}s_{i}(t-\tau_{i})e^{j2\pi\nu_{i}(t-\tau_{i})}+z_{i}(t),$
(5)

where $\mathbf{u}_{i}\in\mathbb{C}^{N_{u}\times 1}$ is the receive beamformer,
$z_{i}(t)$ is the time domain noise signal. We adopt a receive filter $g_{\rm
rx}(t)$ and after OTFS demodulation, the received OTFS signal is expressed
as111For simplicity, we assume integer delay and Doppler shifts in this
letter. For the fractional case, off-grid estimation schemes can be adopted,
e.g., [11].

$\displaystyle Y^{i}_{\rm DD}\left[{l,k}\right]=$ $\displaystyle
h_{i}\mathbf{u}^{\rm
H}_{i}\mathbf{a}_{N_{u}}(\theta_{i})\mathbf{a}_{N_{t}}^{\rm
H}(\theta_{i})\mathbf{f}_{i}X^{i}_{\rm
DD}\left[{(l-l_{i})_{M},(k-k_{i})_{N}}\right]$ $\displaystyle+Z^{i}_{\rm
DD}\left[{l,k}\right],$ (6)

where the integers $k_{i}={\nu_{i}}{NT}$, $l_{i}={\tau_{i}}{M\Delta f}$,
$Z^{i}_{\rm DD}\left[{l,k}\right]$ is the additional white Gaussian noise
sample with a power spectral density of $N_{0}$.

###  II-B Sensing Model

The sensing channel is given by

$\displaystyle\mathbf{H}(t,\tau)=\sum_{i=1}^{K}\gamma_{i}\mathbf{a}_{N_{r}}(\theta_{i})\mathbf{a}_{N_{t}}^{\rm
H}(\theta_{i})\delta(\tau-\eta_{i})e^{j2\pi\upsilon_{i}t},$ (7)

where $\gamma_{i}$, $\eta_{i}$, and $\upsilon_{i}$ are the reflected
coefficient, the delay, and the Doppler associated with the $i$-th target.
Then, the acquired sensing echoes at the BS can be expressed as

$\displaystyle\mathbf{r}(t)=\sum_{i=1}^{K}\gamma_{i}\mathbf{a}_{N_{r}}(\theta_{i})\mathbf{a}_{N_{t}}^{\rm
H}(\theta_{i})\tilde{\mathbf{s}}(t-\eta_{i})e^{j2\pi\upsilon_{i}t}+\mathbf{w}(t),$
(8)

where $\mathbf{w}(t)$ is the measurement noise. For massive MIMO receive
antenna arrays, steering vectors with different angular values are
asymptotically orthogonal [12], i.e., $\mathbf{a}_{N_{r}}^{\rm
H}(\theta_{i})\mathbf{a}_{N_{r}}(\theta_{i^{\prime}})\approx
0,~{}\forall\theta_{i}\neq\theta_{i^{\prime}}$. Thus, the interference between
different targets in the sensing echoes can be neglected and the BS can
distinguish different targets in terms of their angles-of-arrival (AoAs). As a
result, the sensing echo of the $i$-th target can be extracted from
$\mathbf{r}(t)$ with a receive beamformer
$\mathbf{b}_{i}=\mathbf{a}_{N_{r}}(\tilde{\theta}_{i})$, expressed as

$\displaystyle{r}_{i}(t)=\gamma_{i}\mathbf{b}^{\rm
H}_{i}\mathbf{a}_{N_{r}}(\theta_{i})\mathbf{a}_{N_{t}}^{\rm
H}(\theta_{i})\mathbf{f}_{i}s_{i}(t-\eta_{i})e^{j2\pi\upsilon_{i}t}+{w}(t).$
(9)

We note that angular parameter $\theta_{i}$ can also be inferred from the
receive beamforming process by comparing the gains obtained for different beam
directions. Hence, (9) is simplified as

$\displaystyle{r}_{i}(t)=G_{a}s_{i}(t-\eta_{i})e^{j2\pi\upsilon_{i}t}+{w}_{i}(t),$
(10)

where $G_{a}$ is the composite antenna array gain. By adopting an ideal
receive filter and the OTFS demodulation, we finally arrive at the DD domain
input-output relationship,

$\displaystyle R^{i}_{\rm DD}\left[{l,k}\right]=$ $\displaystyle
G_{a}\sum_{k^{\prime}=0}^{N-1}\sum_{l^{\prime}=0}^{M-1}H^{i}_{\rm
DD}[l^{\prime},k^{\prime}]$ (11) $\displaystyle\cdot X^{i}_{\rm
DD}\left[{(l-l^{\prime})_{M},(k-k^{\prime})_{N}}\right]+W^{i}_{\rm
DD}\left[{l,k}\right],$

where $H^{i}_{\rm DD}[l^{\prime},k^{\prime}]$ denotes the gain of the $i$-th
target at the DD grid (bin) with indices $l^{\prime},~{}k^{\prime}$
corresponding to the delay of $\frac{l^{\prime}}{M\Delta f}$ and Doppler of
$\frac{k^{\prime}}{NT}$. By stacking all received samples $R^{i}_{\rm DD}$ in
a vector $\mathbf{r}^{i}_{\rm DD}\left[{l,k}\right]$, (11) is rewritten in a
matrix form as

$\displaystyle\mathbf{r}^{i}_{\rm DD}=G_{a}\mathbf{X}^{i}_{\rm
DD}\mathbf{h}^{i}_{\rm DD}+\mathbf{w}^{i}_{\rm DD},$ (12)

where the $(Nk^{\prime}+l^{\prime})$-th element in $\mathbf{h}^{i}_{\rm DD}$
is $H^{i}_{\rm DD}[l^{\prime},k^{\prime}]$. Based on the estimate of sensing
channel $\mathbf{h}^{i}_{\rm DD}$, it is able to extract the delay and Doppler
associated to the $i$-th target.

##  III Preliminary Results and Discussions

###  III-A Sensing Parameters Estimation

The classic matched filtering scheme in radar community can be used for
obtaining the estimate of $\mathbf{h}^{i}_{\rm DD}$, given by

$\displaystyle\tilde{\mathbf{r}}^{i}_{\rm DD}=\mathbf{X}^{i,{\rm H}}_{\rm
DD}\mathbf{r}^{i}_{\rm DD}=G_{a}\mathbf{X}^{i,{\rm H}}_{\rm
DD}\mathbf{X}^{i}_{\rm DD}\mathbf{h}^{i}_{\rm DD}+\mathbf{X}^{i,{\rm H}}_{\rm
DD}\mathbf{w}^{i}_{\rm DD}.$ (13)

Obviously, (13) is a linear model and it is straightforward to find an
estimator, e.g., maximum likelihood (ML) estimator or linear minimum mean
square error (LMMSE) estimator, to obtain the estimate
$\hat{\mathbf{h}}^{i}_{\rm DD}$. Such estimators are Cramér-Rao bound (CRB)
achieving if the elements of $\mathbf{h}^{i}_{\rm DD}$ are also Gaussian
distributed. However, $\mathbf{h}^{i}_{\rm DD}$ corresponds to the round-trip
delay and Doppler shifts of the $i$-th target, indicating that
$\mathbf{h}^{i}_{\rm DD}$ is a sparse vector. The delay and Doppler indices
can be determined by finding element in $\hat{\mathbf{h}}^{i}_{\rm DD}$ that
has the largest response, i.e., $[\hat{l}_{i},\hat{k}_{i}]$, yielding the
delay and Doppler estimates $\hat{\eta}_{i}$ and $\hat{\phi}_{i}$,
respectively222In addition to the adopted estimator, the sensing accuracy also
depends on the resolutions of the delays and Doppler shifts, which are
associated with the bandwidth and signal duration. Therefore, the OTFS frame
size, i.e., $M$ and $N$ will also affect the estimation performance.. Due to
the existence of noise, the above estimation is biased in the low signal-to-
noise ratio (SNR) regime, which has been reported in [5].

###  III-B Sensing-assisted OTFS Communications

From the estimates of delay $\hat{\eta}_{i}$, Doppler $\hat{\phi}_{i}$, and
angle $\hat{\theta}_{i}$, we are now able to determine the location as well as
the speed of target $i$. Assuming that the BS is located at the origin
$[0,0]^{\rm T}$, the location of the $i$-th target
$\mathbf{p}_{i}=[p_{x,i},p_{y,i}]^{\rm T}$ is written as

$\displaystyle\hat{p}_{x,i}=\frac{c\hat{\eta}_{i}\sin\hat{\theta}_{i}}{2},$
(14)
$\displaystyle{\hat{p}_{y,i}=\frac{c\hat{\eta}_{i}\cos\hat{\theta}_{i}}{2}.}$
(15)

Similarly, the velocity of the $i$-th target
$\hat{\mathbf{v}}_{i}=[\hat{v}_{x,i},\hat{v}_{y,i}]^{\rm T}$ can be obtained
based on $\hat{\phi}_{i}$ and $\hat{\theta}_{i}$.

According to the definition of OTFS channel model, the Doppler $\nu_{i}$ is in
theory half of the Doppler associated with the sensing echo $\phi_{i}$. The
slow time-varying property of the OTFS communication channel motivates us to
reuse the sensing parameters for communications purpose. Relying on the
estimates of target locations as well as the speeds, the locations of the
targets in the following time instant are predictable. In particular, the
predicted location of the $i$-th target can be expressed as

$\displaystyle\tilde{\mathbf{p}}_{i}=\hat{\mathbf{p}}_{i}+\hat{\mathbf{v}}_{i}\Delta
T,$ (16)

where $\Delta T$ is the duration of one transmission block. For the speeds of
the targets, it is naturally to set the predicted velocity
$\tilde{\mathbf{v}}_{i}$ identical to the estimated velocity, since the target
speed will not have a burst change in a relatively short duration. The
predicted locations of targets enable us to further obtain the prediction of
the angular parameter $\theta_{i}$ in the following time instant, based on the
geometric relationship as

$\displaystyle\tilde{\theta}_{i}=\arctan\frac{\tilde{p}_{y,i}}{\tilde{p}_{x,i}}.$
(17)

The predicted angle $\tilde{\theta}_{i}$ can be used for steering the transmit
and receive antennas before transmitting the multi-beam OTFS-ISAC signal and
receiving the echoes. Compared to the conventional beam pairing and alignment
schemes, using the predicted angle avoids the adoption of dedicated pilots for
estimating the angular parameters, which can improve the resource efficiency.
The predicted angle can also be contained in the signal to targets, which is
used for designing the beamformer $\mathbf{u}_{i}$ [13, 14].

In addition, the prediction of the range distance is obtained as
$\tilde{d}_{i}=\|\tilde{\mathbf{p}}_{i}\|$. Consequently, the downlink channel
gain $h_{i}$ can be predicted by replacing ${d}_{i}$ by $\tilde{d}_{i}$.
Moreover, the delay and Doppler indices in the downlink OTFS communications
signal model are determined as well by exploiting the predicted range and
speed. Motivated by the fact that all communication channel parameters can be
predicted at the BS side, it is capable of compensating the channel effects
before downlink transmission. Provided sufficiently accurate prediction of the
communication channel and compensation at the BS, the received OTFS signal is
simplified as

$\displaystyle Y^{i}_{\rm DD}\left[{l,k}\right]=G_{t}X^{i}_{\rm
DD}\left[{l,k}\right]+Z^{i}_{\rm DD}\left[{l,k}\right],$ (18)

with $G_{t}$ being the composite array gain. It is observed from (18) that the
communications receiver can perform data detection while bypassing the channel
estimation.

To verify the effectiveness of the proposed scheme, we show the bit-error-rate
(BER) performance versus the SNR based on the proposed scheme in Fig. 2. The
size of OTFS frame is $M=128$ and $N=20$. The BS is operating at a carrier
frequency of $3$ GHz and the subcarrier spacing is $6$ kHz. Therefore, the
occupied bandwidth is $768$ kHz. The radar cross section (RCS) of the target
is assumed to be $25$, given the reflection area of the target being
$1\textrm{m}^{2}$. Finally, the transmit power is set to $40$ dBm. Moreover,
we consider the binary phase shift keying (BPSK) symbol mapping. The speed of
the target is randomly generated from the uniform distribution $[10,15]$ m/s.
For comparison purpose, the BER performance corresponding to the classic
pilot-based beam alignment and channel estimation scheme, and the ideal case
with perfect knowledge of channel are illustrated. We can see that our
proposed algorithm can approach the performance of the ideal case, which
implies the highly accurate channel prediction. Nevertheless, due to the
random noise, the predicted channel may deviate from the actual channel,
resulting in a small BER performance gap. Compared to the classic pilot-based
scheme, the proposed algorithm yields a significant BER performance gain. This
is because the transmitted signal is known to the BS, providing a high SNR
gain to improve the estimation performance of the sensing parameters. Despite
the performance improvement, the proposed algorithm also reduces the
communication overhead as well as latency of the classic scheme, showing the
advantage of adopting OTFS-ISAC signal.

Figure 2: The BER performance of different algorithms for uplink transmission.

###  III-C Discussions of OTFS-ISAC system

So far, we have proposed an example of sensing-assisted communications
protocol for OTFS-ISAC systems. We show that benefiting from the slow-varying
OTFS channel, the parameters extracted from sensing echoes can be exploited
for predicting the downlink DD domain communication channel. As a result, no
pilots are required either for beam alignment or for channel estimation,
leading to a reduced communication overhead. Moreover, since the angle
feedback and channel estimation are bypassed, the communications receiver is
capable of directly detecting the data symbols, which is attractive for
supporting ultra low latency communications in the future networks. Relying on
the same DD domain for both communications and sensing, OTFS offers more
promising benefits for ISAC systems, which will be briefly discussed in the
following.

####  III-C1 Unified Signaling Design Framework

Let us consider a general single antenna case. When transmitting a vectorized
OTFS signal $\mathbf{x}_{\rm DD}$, after OTFS demodulation, the acquired
sensing echo and the received OTFS signal can be written in concise forms as

$\displaystyle\mathbf{r}_{\rm DD}$ $\displaystyle=\mathbf{H}_{s,{\rm
DD}}\mathbf{x}_{\rm DD}+\mathbf{w}_{\rm DD},$ (19)
$\displaystyle\mathbf{y}_{\rm DD}$ $\displaystyle=\mathbf{H}_{c,{\rm
DD}}\mathbf{x}_{\rm DD}+\mathbf{z}_{\rm DD},$ (20)

where $\mathbf{H}_{s,{\rm DD}}$ and $\mathbf{H}_{c,{\rm DD}}$ denotes the
sensing and communication channels in DD domain, respectively. For sensing,
the knowledge of $\mathbf{x}_{\rm DD}$ at the BS is exploited for inferring
the sensing parameters $\bm{\eta}$ from channel $\mathbf{H}_{s,{\rm DD}}$. In
contrast, for communications, the transmitted data symbols $\mathbf{x}_{\rm
DD}$ are detected. The sensing performance is characterized by ambiguity
function and the CRB while the communication performance can be evaluated
using channel capacity.  In particular, the CRB of the sensing parameters
$\bm{\eta}$ is derived as [15],

$\displaystyle\textrm{CRB}_{\bm{\eta}}=\textrm{Tr}\left[\mathbb{E}\left[-\frac{\partial^{2}\ln
p(\mathbf{r}_{\rm DD}|\bm{\eta})}{\partial\bm{\eta}^{2}}\right]^{-1}\right],$
(21)

where $p(\mathbf{r}_{\rm DD}|\bm{\eta})$ is the likelihood function, and the
capacity of the communication channel $C_{\mathbf{H}_{c}}$ can be expressed as

$\displaystyle
C_{\mathbf{H}_{c}}=\log_{2}\textrm{det}\left(\mathbf{I}_{MN}+\frac{1}{N_{0}}\mathbf{H}_{c,\rm{DD}}\mathbf{R}_{x}\mathbf{H}_{c,\rm{DD}}^{\rm
H}\right),$ (22)

where $\mathbf{R}_{x}$ is the covariance matrix of the transmitted signal. The
OTFS waveform enables us to formulate the signaling design problem in the same
DD domain for achieving satisfactory sensing and communications performance.
For example, in a communication-centric scenario, we aim for maximizing the
channel capacity $C_{H_{c}}$ while maintaining the sensing performance. Then,
it is straightforward to formulate the following optimization problem

$\displaystyle\max_{\mathbf{x}_{\rm DD}}C_{\mathbf{H}_{c}}$ (23)
$\displaystyle\mbox{s.t.}\;\;\mbox{{C1}}:\;$
$\displaystyle\textrm{CRB}_{\bm{\eta}}\leq T_{\rm CRB}$
$\displaystyle\mbox{{C2}}:\;$ $\displaystyle\frac{1}{MN}\|\mathbf{x}_{\rm
DD}\|=P_{T},$

where C1 indicates the acceptable sensing performance and C2 constraints the
power consumption. A possible solution of (23) requires a specific waveform
design, which can achieve the optimal tradeoff between sensing and
communication via taking into account their respectively optimal waveforms in
the DD domain. Moreover, the objective of (23) can be extended to include both
communication and sensing metrics for optimizing OTFS signal that can strike a
good trade-off between both functionalities. Instead of directly designing
$\mathbf{x}_{\rm DD}$, we can also design the DD domain precoding matrix
$\mathbf{P}_{\rm DD}$ through optimization given certain constraints such that
the transmitted OTFS signal becomes $\mathbf{P}_{\rm DD}\mathbf{x}_{\rm DD}$

####  III-C2 Sensing Channel Exploitation

As discussed in Section III, the delay and Doppler parameters associated with
the sensing channel are reused for OTFS communication channel prediction. In
this case, we assume that all targets of interest are communication users as
well. However, this is not always the case in other scenarios, e.g., only part
of the targets in the area of coverage are UEs. Therefore, a fundamental
question is how much information we can glean from the sensing channel
$\mathbf{H}_{s,{\rm DD}}$ for inferring $\mathbf{H}_{c,{\rm DD}}$. As the
first step, we need to characterize the similarity between OTFS communication
and sensing channels. A nature choice is to use the cross-correlation matrix
$\mathbf{R}_{cs}=\mathbb{E}\left[\mathbf{H}_{c,{\rm DD}}\mathbf{H}_{s,{\rm
DD}}^{\rm H}\right]$ and the entries of $\mathbf{R}_{cs}$ will provide the
correlation between sensing targets and communication reflectors corresponding
to different paths. In particular, when sensing and communication channels are
highly correlated, the matrix $\mathbf{R}_{cs}$ will be diagonally dominant.
Although OTFS modulation offers the opportunity to connect two channels in the
same DD domain, new analysis tools and framework are still to be developed to
further study the exploitation of sensing channel in OTFS-ISAC system.

####  III-C3 Communications-assisted Sensing

The OTFS channel model reveals the underlying geometry of the wireless
propagation environments. After channel estimation, the estimated delay and
Doppler indices can be converted to range distance and speed. Given a temporal
sequence of estimated delays corresponding to the $i$-th path, i.e.,
$\bm{\tau}_{i}=[{\hat{\tau}_{i,1},...\hat{\tau}_{i,t}}]^{\rm T}$, the location
of the $i$-th reflector can be obtained via an ML estimator, i.e.,

$\displaystyle\hat{\mathbf{p}}_{i}=\arg\max_{{\mathbf{p}}_{i}}p(\bm{\tau}_{i}|\mathbf{p}_{i}),$
(24)

where $p(\bm{\tau}_{i}|\mathbf{p}_{i})$ is the likelihood function of measured
delays conditioned on the location of the $i$-th reflector. If the reflector
is also a moving target, Bayesian estimator such as extended Kalman filtering
(EKF) can be used to exploit the state transition information. By combining
the sensing information gleaned from the reflected echoes at BS and from the
OTFS channel estimation, we can obtain better knowledge of the environment.

##  IV Potential Applications of OTFS

In addition to unifying communications and sensing in the same DD domain, OTFS
is also attractive for several other emerging applications, due to its
advantages of low complexity implementation, sparse channel representation,
and high resilience to large delay and Doppler spread.

###  IV-A Underwater Communications

Due to the propagation medium in underwater communications, mechanical waves,
e.g., acoustic wave and seismic wave instead of the radio signal, are usually
used for information transmission. The relatively low signal propagation speed
of mechanical waves results in large delay and Doppler spreads for underwater
communication channels. In general, the delay could spread up to hundreds
milliseconds, which is much higher than typical wireless radio channels having
the delay level of $\mu$s. The Doppler spread is also severe due to internal
waves and the movements of communication transceivers. Moreover, multiple
reflectors from the sea-surface and sea floor lead to complex multi-path
propagation environments. On the other hand, the inherent wideband propagation
in underwater communications will introduce Doppler-scale effect instead of a
simple Doppler shift. To characterize different delays and Doppler scales, the
multiscale multilag (MSML) channel model is usually used, expressed as

$\displaystyle
h(t,\tau)=\sum_{i=1}^{P}h_{i}\delta(\tau-\tau_{i}-\alpha_{i}t)e^{-j\alpha_{i}f_{c}t},$
(25)

where $P$ is the number of paths, $\alpha_{p}$ is the time scaling parameter.
In fact, in underwater environments, the received signal is a time-delayed,
time-scaled, and phase-rotated version of the transmitted signal caused by
moving scatterings. Although OFDM has been applied in underwater
communications because of its robustness against long delay spread, the inter-
carrier interference (ICI) induced by multiple Dopplers will deeply degrade
the performance. Compared with OFDM modulation, OTFS modulation could benefit
from the possible sparsity and stability of the DD domain channel response,
enabling efficient channel estimation and equalization.

###  IV-B Optical Transmission

The optical spectrum offers huge and unregulated bandwidth without
electromagnetic interference. Information symbols are transmitted via
modulating the luminous intensity of light-emitting diodes (LED) while the
received signal is converted to electrical signal using photodiodes. In
optical transmission, the channel is generally modeled as static multi-path
one. While OTFS is recognized for its superiority in high Doppler scenarios,
its capability of exploiting full time-frequency diversity will guarantee its
performance in static channels as well. Recent results have shown that OTFS-
based optical transmission system achieves a better performance than the
conventional OFDM-based counterpart with a single LED [16]. It is worth to
extend the research to multi-LED and multi-photodiodes systems based on OTFS
modulation.

###  IV-C Space-air-ground Integrated Networks (SAGIN)

To provide ubiquitous connectivity and global coverage, the future
communications will be SAGIN. The low-earth-orbit (LEO) satellite
communications, the aerial communications, and the vehicular communications
will extend the coverage range of the current deployed terrestrial networks.
On the one hand, extremely high-mobility is the principal characteristic for
the satellite and space communications. Such a high speed will inevitably
introduce a very high Doppler shift, which is challenging for reliable
information transmission. On the other hand, the networking with different
communication devices will result in complex signal transmission scenarios.
For example, in high-speed train communications, the environments will change
frequently, such as viaducts, cuttings, and tunnels, which imposes great
difficulty for robust communications. Relying on the quasi-periodic property,
OTFS modulation can support mobile users having a very wide range of speeds.
Moreover, it is highly symmetrical to any channel distortions caused by delay
and Doppler spreads, making OTFS a good solution for realizing the SAGIN. Some
recent contributions have applied OTFS in LEO system and have verified its
effectiveness [17]. However, how to design OTFS parameters, such as delay and
Doppler resolutions, requires extensive channel measurements and
investigations.

##  V Conclusions

In the last part of this tutorial on OTFS modulation, we aim for unveiling its
great potentials of supporting future networks. We first studied the OTFS-
based ISAC system, which shows the superiority of sensing-assisted OTFS
communications. Then several possible benefits of adopting OTFS as the ISAC
signaling waveform were discussed. As a further step, we summarized several
OTFS modulation-related applications including underwater communications,
optical transmission, and SAIGN, which are interesting and of great
importance, but have not been completely investigated.

## References

  * [1] Z. Wei, W. Yuan, S. Li, J. Yuan, G. Bharatula, R. Hadani, and L. Hanzo, “Orthogonal time-frequency space modulation: A promising next-generation waveform,” _IEEE Wireless Commun._ , vol. 28, no. 4, pp. 136–144, Aug. 2021\.
  * [2] G. N. Saddik, R. S. Singh, and E. R. Brown, “Ultra-wideband multifunctional communications/radar system,” _IEEE Trans. Microw. Theory Techn._ , vol. 55, no. 7, pp. 1431–1437, Jul. 2007.
  * [3] F. Liu, C. Masouros, A. Petropulu, H. Griffiths, and L. Hanzo, “Joint radar and communication design: Applications, state-of-the-art, and the road ahead,” _IEEE Trans. Commun._ , vol. 68, no. 6, pp. 3834–3862, Jun. 2020\.
  * [4] P. Raviteja, K. T. Phan, Y. Hong, and E. Viterbo, “Orthogonal time frequency space (OTFS) modulation based radar system,” in _IEEE Radar Conf._ , 2019, pp. 1–6.
  * [5] L. Gaudio, M. Kobayashi, G. Caire, and G. Colavolpe, “On the effectiveness of OTFS for joint radar parameter estimation and communication,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 9, pp. 5951–5965, Sep. 2020.
  * [6] L. Gaudio, M. Kobayashi, G. Caire, and G. Colavope, “Hybrid digital-analog beamforming and MIMO radar with OTFS modulation,” _arXiv preprint arXiv:2009.08785_ , 2020.
  * [7] W. Yuan, Z. Wei, S. Li, J. Yuan, and D. W. K. Ng, “Integrated sensing and communication-assisted orthogonal time frequency space transmission for vehicular networks,” _IEEE J. Sel. Topics Signal Process._ , vol. 15, no. 6, pp. 1515–1528, Dec. 2021.
  * [8] S. Li, W. Yuan, C. Liu, Z. Wei, J. Yuan, B. Bai, and D. W. K. Ng, “A novel ISAC transmission framework based on spatially-spread orthogonal time frequency space modulation,” __to appear in_ IEEE J. Sel. Areas Commun._, 2022.
  * [9] K. Wu, J. A. Zhang, X. Huang, and Y. J. Guo, “OTFS-based joint communication and sensing for future industrial IoT,” __to appear in_ IEEE Internet Things J._, 2021.
  * [10] W. Yuan, Z. Wei, J. Yuan, and D. W. K. Ng, “A simple variational Bayes detector for orthogonal time frequency space (OTFS) modulation,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 7, pp. 7976–7980, Jul. 2020.
  * [11] Z. Wei, W. Yuan, S. Li, J. Yuan, and D. W. K. Ng, “Off-grid channel estimation with sparse bayesian learning for OTFS systems,” _IEEE Trans. Wireless Commun._ , in press, 2022.
  * [12] T. L. Marzetta, _Fundamentals of massive MIMO_. Cambridge University Press, 2016.
  * [13] F. Liu, W. Yuan, C. Masouros, and J. Yuan, “Radar-assisted predictive beamforming for vehicular links: Communication served by sensing,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 11, pp. 7704–7719, Oct. 2020\.
  * [14] W. Yuan, F. Liu, C. Masouros, J. Yuan, D. W. K. Ng, and N. González-Prelcic, “Bayesian predictive beamforming for vehicular networks: A low-overhead joint radar-communication approach,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 3, pp. 1442–1456, Mar. 2021.
  * [15] S. M. Kay, _Fundamentals of statistical signal processing: estimation theory_. Prentice-Hall, Inc., 1993.
  * [16] A. Sharma, S. Jain, R. Mitra, and V. Bhatia, “Performance analysis of OTFS over mobile multipath channels for visible light communication,” in _IEEE Region 10 Conf._ , 2020, pp. 490–495.
  * [17] J. Shi, J. Hu, Y. Yue, X. Xue, W. Liang, and Z. Li, “Outage probability for OTFS based downlink LEO satellite communication,” __to appear in_ IEEE Trans. Veh. Technol._, 2022.


# Performance Evaluation and Comparison of a New Regression Algorithm

Sabina Gooljar, Kris Manohar and Patrick Hosein
Department of Computer Science, The University of the West Indies, St.
Augustine, Trinidad
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

In recent years, Machine Learning algorithms, in particular supervised
learning techniques, have been shown to be very effective in solving
regression problems. We compare the performance of a newly proposed regression
algorithm against four conventional machine learning algorithms namely,
Decision Trees, Random Forest, $k$-Nearest Neighbours and XG Boost. The
proposed algorithm was presented in detail in a previous paper but detailed
comparisons were not included. We do an in-depth comparison, using the Mean
Absolute Error (MAE) as the performance metric, on a diverse set of datasets
to illustrate the great potential and robustness of the proposed approach. The
reader is free to replicate our results since we have provided the source code
in a GitHub repository while the datasets are publicly available.

##  1 INTRODUCTION

Machine Learning algorithms are regularly used to solve a plethora of
regression problems. The demand for these algorithms has increased
significantly due to the push of digitalisation, automation and analytics.
Traditional techniques such as Random Forest, Decision Trees and XG Boost have
been integral in various fields such as banking, finance, healthcare and
engineering. Technology is always evolving and technological advancements are
driven by factors such as human curiosity, problem-solving and the desire for
increased efficacy and reliability. Researchers are constantly working on
improving these existing methods as well as exploring new improved strategies
as can be seen in [Hosein, 2022]. This approach uses a distance metric
(Euclidean distance) and a weighted average of the target values of all
training data points to predict the target value of a test sample. The weight
is inversely proportional to the distance between the test point and the
training point, raised to the power of a parameter $\kappa$. In our paper, we
investigate the performance of this novel approach and several well-
established machine learning algorithms namely XG Boost, Random Forest,
Decision Tree and $k$-NN using the Mean Absolute Error (MAE) as the
performance metric. We intend to showcase the potential of this new algorithm
to solve complex regression tasks across diverse datsets. In the next section,
we describe related work and then the theory of the proposed approach. After,
we present and discuss the findings such as any issues encountered. Finally,
we advocate that the proposed approach may be robust and efficient making it
extremely beneficial to the field.

##  2 RELATED WORK AND CONTRIBUTIONS

In this section, we summarize the various regression techniques we considered
and then discuss differences with the proposed approach. Our contribution,
which is a detailed comparison, is then outlined.

###  2.1 Decision Tree

Decision Tree is a supervised machine learning algorithm that uses a set of
rules to make decisions. Decision Tree algorithm starts at the root node where
it evaluates the input features and selects the best feature to split the data
[Quinlan, 1986]. The data is split in such a way so that it minimises some
metric that quantifies the difference between the actual values and the
predicted values such as Mean Squared Error and Sum of Squared error. Then a
feature and a threshold value is chosen that best divides the data into two
groups. The data is split recursively into two subsets until a stopping
condition is met such as having too few samples in a node. When the decision
tree is constructed, predictions are made by traversing from the root node to
a leaf node. The predicted value is calculated as the mean of the target
values in the training samples which is associated with that leaf node.

###  2.2 Random Forest

Random Forest builds decision trees on different samples and then averages the
outputs for regression tasks [Breiman, 2001]. It works on the principle of an
ensemble method called Bagging. Ensemble is combining multiple models and then
the set of models is used to make predictions instead of using an individual
model. Bagging which is also known as Bootstrap Aggregation selects random
samples from the original dataset. Each model is created from the samples that
are given by the original data with replacement. Individual decision trees are
constructed for each sample and each tree then produces its own output. These
outputs are numerical values. The final output is then calculated from the
average of these values which is known as aggregation.

###  2.3  k-Nearest Neighbours

k-Nearest Neighbours (kNN) is a supervised machine learning algorithm that is
used to solve both classification and regression tasks. Firstly, choose the
number of neighbours (k) which is used when making predictions. Then it
calculates the distances (Euclidean) between a new query and the existing data
and selects the specified number of neighbours (k) that is closest to the
query and finds the average of these values. The average is the predicted
value.

###  2.4 XGBoost

XGBoost (eXtreme Gradient Boosting) is an ensemble learning algorithm that
combines the output of weak learners (usually decision trees) to make more
accurate predictions [Chen and Guestrin, 2016]. New weak learners are added to
the model iteratively with each tree aiming to correct the errors made by the
previous learners. The training process is stopped when there is no
significant improvement in a predefined number of iterations.

###  2.5 New Algorithm

Regression models enable decision-making in a wide range of applications such
as finance, healthcare, education and engineering. It is imperative that these
regression models are precise and robust so that better decisions can be made
to enhance and improve these fields. While there are various popular machine
learning algorithms for solving regression tasks, we introduce a new
regression model that shows high accuracy and robustness, ensuring that real-
world applications are optimised. The core of the approach is similar to
$k$-NN but instead of using samples in a neighbourhood, all samples are used
and closer samples are weighted more heavily than those further away. In this
case, there is no parameter $k$ to specify but we do introduce a parameter
$\kappa$ that dictates the rate of decay of the weighting.

##  3 PROPOSED APPROACH

The proposed approach was originally designed to determine a suitable
insurance policy premium [Hosein, 2022]. Specifically, [Hosein, 2022] noted as
personlisation increases (i.e., more features), predictions become less robust
due to the reduction in the number of samples per feature, especially in
smaller datasets. His main goal was to achieve an optimal balance between
personlisation and robustness. Instead of using the samples available for each
feature, his algorithm computes the weighted average of the target variable
using all samples in the dataset. This algorithm uses the Euclidean distance
metric and a hyperparameter $\kappa$, which controls the influence of the
distance (the weights) between points in the data. Another aspect is that the
same unit of distance is used for each feature which allows one, for example,
to compare distance between a gender feature with the distance between an age
feature.

The $\kappa$ parameter introduced in [Hosein, 2022] is used as an exponent in
the weighting formula, where the weights are inversely proportional to the
distance between data points raised to the power of kappa. When $\kappa$ is
large, the influence of points further away from the test point decreases
quickly since their distance raised to a large power becomes very large which
in turn makes the weight very small. However, when $\kappa$ is small, the
influence of points further away decreases slowly since the distance raised to
a small power results in a relatively smaller value which then makes the
inverse weight larger.

The algorithm firstly normalises the ordinal features. Then the prediction is
done in two parts. For example, say we have a single categorical feature
Gender with two values, Male and Female. In the first stage, we compute the
mean for each feature value over all the training samples. That is, the
average target value for all females $\mu_{Gender,Female}$ and the average
target value for all males $\mu_{Gender,male}$. With these means, the distance
d between a Female and a Male is the distance between $\mu_{Gender,Female}$
and $\mu_{Gender,male}$. The second stage computes the prediction for a test
sample. For a given test sample i, its prediction is the weighted average of
the target value over all training samples. These weights are computed as
$\frac{1}{(1+d[i,j])^{\kappa_{2}}}$, where $d[i,j]$ is the Euclidean distance
between the test sample i and the training sample j.

Numerical features pose an interesting challenge because they can have a wider
(potentially infinite) range of values. For example, consider the feature of
Age. Compared to Gender, this feature can easily span 40 values (e.g., 18 to
58) instead of two. Thus, for the same training set, the number of samples per
age will be low which implies the means (i.e., $\mu_{Age,20}$, $\mu_{Age,21}$,
$\mu_{Age,30}$ etc.,) used to compute distances (and hence the weights) may
not be robust. Additionally, we may encounter some unique values in the
testing data set that are not in the training data set or vice versa.
Generally, categorical variables do not encounter these issues because the
number of samples per category value is sufficient. In order to solve this
problem for numerical features, we impute a value for means for each unique
value in both the training and testing data sets.

These imputed means are calculated based on the distances between the
attribute value and all training samples in the feature space. This distance
assists the algorithm in determining which training samples are most relevant
to the test point. Then the target values for these training samples are
combined, using a weighted average similar to before. The weights are computed
as $\frac{1}{(1+d_{f}[i,j])^{\kappa_{1}}}$, where $d_{f}[i,j]$ is the distance
between the numerical feature values in test sample i and training sample j.
For example, say f = Age, and test sample i’s age is 30 and training sample
j’s age is 40 then $d_{Age}[i,j]=|30-40|=10$.

This modification means our proposed approach uses two kappa values: one for
the pre-processing (i.e., $\kappa_{1}$) and one for predicting (i.e.,
$\kappa_{2}$). We determine the optimal combination of kappa values
($\kappa_{1},\kappa_{2})$ that minimises the error. According to [Hosein,
2022], as $\kappa_{2}$ increases, the error decreases up to a certain point
and then the error increases after this point. Therefore, an optimal
$\kappa_{2}$ can be found that minimizes the error.

We define a range of values for both the pre-processing and predicting parts.
The initial range is determined through a trial and error process. We observe
the MAE and adjust these values if needed. However, note that while the
initial range of kappa values involves some trial and error, the process of
finding the optimal combination within the range is essentially a grid search
which is systematic and data-driven and ensures that the model is robust.
Since the algorithm uses all the training data points in its prediction, it
will be robust for small data sets or where there are not enough samples per
category of a feature. The Pseudo code in Figure 1 summarizes the steps of the
proposed algorithm.

$\displaystyle 1:\quad C\equiv\text{set of categorical features}$
$\displaystyle 2:\quad O\equiv\text{set of ordinal features}$ $\displaystyle
3:\quad X\equiv\text{set of training samples}$ $\displaystyle 4:\quad
Y\equiv\text{set of testing samples}$ $\displaystyle
5:\quad\kappa_{1},\kappa_{2}>0\quad\text{tuning parameters}$ $\displaystyle
6:\quad\text{for each }f\in O\text{ do}$ $\displaystyle 7:\quad\quad\text{for
each sample }j\text{ in }X\text{ do}$ $\displaystyle 8:\quad\quad\quad
x_{\text{train},j,f}\leftarrow\frac{x_{\text{train},j,f}}{\max(f)-\min(f)}\quad\text{(normalize
feature values in the train set)}$ $\displaystyle 9:\quad\quad\text{end for}$
$\displaystyle 10:\quad\quad\text{for each sample }i\text{ in }Y\text{ do}$
$\displaystyle 11:\quad\quad\quad
x_{\text{test},i,f}\leftarrow\frac{x_{\text{test},i,f}}{\max(f)-\min(f)}\quad\text{(normalize
feature values in the test set)}$ $\displaystyle 12:\quad\quad\text{end for}$
$\displaystyle 13:\quad\text{end for}$ $\displaystyle 14:\quad
vf\equiv\text{set of categories for feature }f\in C$ $\displaystyle
15:\quad\text{for each }f\in C\text{ do}$ $\displaystyle
16:\quad\quad\text{for each }v\in vf\text{ do}$ $\displaystyle
17:\quad\quad\quad z\equiv\\{x\in X|x_{f}=v\\}$ $\displaystyle
18:\quad\quad\quad\mu_{f,v}\leftarrow\frac{1}{|z|}\sum_{x\in
z}y_{x}\quad(\text{mean target value over training samples where feature
}f\text{ has value }v)$ $\displaystyle 19:\quad\quad\text{end for}$
$\displaystyle 20:\quad\quad\text{Replace category values with their mean
}\mu_{f,v}\text{ in both }X\text{ and }Y$ $\displaystyle 21:\quad\text{end
for}$ $\displaystyle 22:\quad\text{for each}f\in O\text{ do}$ $\displaystyle
23:\quad\quad\text{for sample }i\text{ with unique feature value }v\text{ in
}X\text{ and }Y\text{ do}$ $\displaystyle
24:\quad\quad\quad\mu_{f,v}\leftarrow\frac{\sum_{j\in
X}\frac{y_{j}}{(1+d_{f}[i,j])^{\kappa_{1}}}}{\sum_{j\in
X}\frac{1}{(1+d_{f}[i,j])^{\kappa_{1}}}}\quad(\text{imputed mean target value
over training samples when feature }f\text{ has value }v)$ $\displaystyle
25:\quad\quad\text{end for}$ $\displaystyle 26:\quad\quad\text{Replace feature
values with the imputed mean values }\mu_{f,v}\text{ in both }X\text{ and }Y$
$\displaystyle 27:\quad\text{end for}$ $\displaystyle 28:\quad\text{for each
test sample }i\text{ in }Y\text{ do}$ $\displaystyle 29:\quad\quad\text{for
each training sample }j\text{ in }X\text{ do}$ $\displaystyle
30:\quad\quad\quad d[i,j]\leftarrow\left(\sum_{f\in
F}(x_{\text{test},i,f}-x_{\text{train},j,f})^{2}\right)^{\frac{1}{2}}\quad\text{(calculate
Euclidean distance) }$ $\displaystyle 31:\quad\quad\text{end for}$
$\displaystyle 32:\quad\text{end for}$ $\displaystyle 33:\quad\text{for each
test sample }i\text{ in }Y\text{ do}$ $\displaystyle 34:\quad\quad
c[i]\leftarrow\frac{\sum_{j\in
X}\frac{y_{j}}{(1+d[i,j])^{\kappa_{2}}}}{\sum_{j\in
X}\frac{1}{(1+d[i,j])^{\kappa_{2}}}}$ $\displaystyle 35:\quad\text{end for}$
$\displaystyle 36:\quad\text{Return the output values }c$ Figure 1: Pseudo
code for the algorithm.

##  4 NUMERICAL RESULTS

In this section, we describe the datasets that were used and apply the various
techniques in order to compare their performances. A GitHub repository,
[Gooljar, 2023] containing the code used in this assessment has been created
to facilitate replication and validation of the results by readers.

###  4.1 Data Set Description

The data sets used were sourced from the University of California at Irvine
(UCI) Machine Learning Repository. We used a wide variety of data sets to
illustrate the robustness of our approach. We removed samples with any missing
values and encoded the categorical variables. No further pre-processing was
done so that the results can be easily replicated. Table 1 shows a summary of
the data sets used.

Table 1: Summary of Data sets. Dataset | No. of Samples | No. of Attributes | Target Value | Citation
---|---|---|---|---
Student Performance | 394 | 32 | G3 | [Cortez, 2014]
Auto | 392 | 6 | mpg | [Quinlan, 1993]
Energy Y2 | 768 | 8 | Y2 | [Tsanas, 2012]
Energy Y1 | 768 | 8 | Y1 | [Tsanas, 2012]
Iris | 150 | 4 | Sepal Length | [Fisher, 1988]
Concrete | 1030 | 8 | Concrete Compressive Strength | [Yeh, 2007]
Wine Quality | 1599 | 11 | Residual Sugar | [Cortez et al., 2009]

###  4.2 Feature Selection

There are various ways to perform feature selection [Banerjee, 2020] but the
best subset of features can only be found by exhaustive search. However, this
method is computationally expensive so we select the optimal subset of
features for the Random Forest model using Recursive Feature Elimination with
Cross-Validation(RFECV) [Brownlee, 2021] and use these features for all other
models. Note that for each model, there may be a different optimal subset of
features and, in particular, this subset of features may not be optimal for
the proposed approach so it is not provided with any advantage. Table 2 shows
the selected attributes for each dataset. The optimal subset of features was
the full set of features for Auto, Energy Y2 and Iris datasets. The columns
are indexed just as they appear in the datasets from UCI.

Table 2: Summary of Selected Features (categorical features are color coded in red) Dataset |  Selected Attributes
---|---
Student Performance |  [2,23,28,29,31]
Auto |  [1,2,3,4,5,6]
Energy Y2 |  [0,1,2,3,4,5,6,7]
Energy Y1 |  [0,1,2,3,4,6]
Iris |  [1,2,3, 4]
Concrete |  [0,1,2,3,4,6]
Wine Quality |  [0,4,5,6,9]

###  4.3 Performance Results

We show the performances of the different algorithms Random Forest, Decision
Tree, k-Nearest Neighbors, XG Boost, and the proposed method. The models are
evaluated on seven datasets (Auto, Student Performance, Energy Y2, Energy Y1,
Iris, Concrete, and Wine Quality). We used Mean Absolute Error (MAE) to
measure the performance since it is robust and easy to interpret. It measures
the average magnitude of errors between the predicted values and the actual
values. Lower MAE values indicate better performance. We used 10 Fold cross-
validation for each model to ensure that we achieve a more reliable estimation
of the model’s performance.

The performance results (MAE) for each algorithm with each dataset are
provided in Table 3. The average MAE for the proposed algorithm is 1.380,
while the next best performer is XG Boost, with an average MAE of 1.653.
Random Forest, Decision Tree, and k-NN have higher average MAEs of 1.659,
2.105 and 2.537, respectively. The proposed algorithm consistently performs
better than the popular algorithms. We also provide a bar chart showing a
comparison between the algorithms for the various datasets in Figure 2.

Table 3: Comparison of algorithm performance on various datasets Dataset |  Algorithm |  MAE
---|---|---
Auto |  Random Forest |  2.048
Decision Tree |  2.684
k-NN |  3.271
XG Boost |  2.141
Proposed with $\kappa_{1}=11,\kappa_{2}=6$ |  1.912
Student Performance |  Random Forest |  1.529
Decision Tree |  1.717
k-NN |  1.868
XG Boost |  1.804
Proposed with $\kappa_{1}=8,\kappa_{2}=8$ |  1.045
Energy Y2 |  Random Forest |  1.501
Decision Tree |  1.678
k-NN |  2.055
XG Boost |  1.275
Proposed with $\kappa_{1}=12,\kappa_{2}=5$ |  1.118
Energy Y1 |  Random Forest |  0.587
Decision Tree |  0.545
k-NN |  1.523
XG Boost |  0.461
Proposed with $\kappa_{1}=2,\kappa_{2}=24$ |  0.347
Iris |  Random Forest |  0.356
Decision Tree |  0.413
k-NN |  0.348
XG Boost |  0.416
Proposed with $\kappa_{1}=7,\kappa_{2}=50$ |  0.269
Concrete |  Random Forest |  5.009
Decision Tree |  6.944
k-NN |  7.970
XG Boost |  4.818
Proposed with $\kappa_{1}=7,\kappa_{2}=10$ |  4.382
Wine Quality |  Random Forest |  0.586
Decision Tree |  0.755
k-NN |  0.726
XG Boost |  0.655
Proposed with $\kappa_{1}=0.5,\kappa_{2}=13$ |  0.584
Average |  Random Forest |  1.659
Decision Tree |  2.105
k-NN |  2.537
XG Boost |  1.653
Proposed |  1.380
Figure 2: Mean Absolute Error for each Model across all Datasets.

###  4.4 Run-Time Analysis

The proposed model performs well against all other models used in this
comparison but it requires more computational time. We did some run time
testing for the datasets for the various algorithms using time.perf_counter()
which is a function in the ‘time’ module of Python’s standard library. It is
used to measure the time a block of code typically takes to run. The function
returns the value, in fractional seconds, of a high-resolution timer. On
average, the proposed approach took approximately 1.56 times longer than XG
Boost and about 274 times longer than the other algorithms. However, we can
perform computation optimizations that will reduce the run-time significantly.
We plan to explore ways to more efficiently determine the optimal $\kappa$
values.

###  4.5 Parameter Optimization

According to [Hosein, 2022], the error appears to be a convex function of
$\kappa$. As $\kappa$ increases, the error decreases up to a certain point and
then the error increases after this point. Therefore, an optimal value can be
found that minimizes the error. Our approach uses two different $\kappa$
values, one for imputing the values and one for predicting. This allows for
more flexibility in the model since each value serves a different purpose and
allows them to optimize both aspects independently to minimize the error. In
some cases, the optimal $\kappa$ value for the first part may not be best for
the predicting part and vice versa. In Figure 3, we can clearly see the
pattern of the error. The MAE decreases to a minimum and then increases after
$\kappa_{2}$ = 8. In our approach, we used a single value for all the features
when imputing. However, we note that the parameter can be further optimised by
using different values for different features. This further optimisation will
lead to an even more complex model but may yield better results so we intend
to explore this in the future.

Figure 3: MAE vs $\kappa$ values for the Student Dataset

##  5 DISCUSSION

We compared the Random Forest, Decision Tree, k-Nearest Neighbours (k-NN), XG
Boost and the proposed algorithm on seven diverse datasets from UCI. The
datasets were from various fields of study and consist of a combination of
categorical and ordinal features. Our proposed approach uses two hyper-
parameters to optimize predictions. The average mean absolute error of the
proposed approach is 45.6 % lower than k-NN, 34.4% lower than Decision Tree,
16.8% lower than the Random Forest and 16.5% lower than XG Boost. The proposed
approach achieves the lowest MAE for all datasets. These results illustrate
the value and potential of the proposed approach.

##  6 CONCLUSIONS AND FUTURE WORK

We present a robust approach that can be used for any regression problem. The
approach is based on a weighted average of the target values of the training
points where the weights are determined by the inverse of the Euclidean
distance between the test point and the training points raised to the power of
a parameter $\kappa$. As shown in Figure 2 the proposed algorithm surpasses
the traditional algorithms in each dataset. Its performance indicates that the
proposed method is a promising approach for solving regression tasks and
should be considered as a strong candidate for future applications. However,
there is significant room for improvement of this algorithm. Future work can
include using different $\kappa$ values for each feature and exploring
heuristic methods to determine these values which may result in even better
performance. Also, since the algorithm’s computations can be done in parallel
(i.e., the grid search over the $\kappa$ space), the run time can also be
considerably decreased.

## REFERENCES

  * Banerjee, 2020 Banerjee, P. (2020).  Comprehensive guide on feature selection.
  * Breiman, 2001 Breiman, L. (2001).  Random forests.  Machine learning, 45(1):5–32.
  * Brownlee, 2021 Brownlee, J. (2021).  Recursive feature elimination (rfe) for feature selection in python.
  * Chen and Guestrin, 2016 Chen, T. and Guestrin, C. (2016).  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 785–794. ACM.
  * Cortez, 2014 Cortez, P. (2014).  Student Performance.  UCI Machine Learning Repository.  DOI: 10.24432/C5TG7T.
  * Cortez et al., 2009 Cortez, P., Cerdeira, A., Almeida, F., Matos, T., and Reis, J. (2009).  Wine Quality.  UCI Machine Learning Repository.  DOI: 10.24432/C56S3T.
  * Fisher, 1988 Fisher, R. A. & Fisher, R. (1988).  Iris.  UCI Machine Learning Repository.  DOI: 10.24432/C56C76.
  * Gooljar, 2023 Gooljar, S. (2023).  Comparison of the performance of a proposed algorithm to other algorithms.  GitHub repository.
  * Hosein, 2022 Hosein, P. (2022).  A data science approach to risk assessment for automobile insurance policies.  International Journal of Data Science and Analytics, pages 1–12.
  * Quinlan, 1986 Quinlan, J. R. (1986).  Induction of decision trees.  Machine learning, 1(1):81–106.
  * Quinlan, 1993 Quinlan, R. (1993).  Auto MPG.  UCI Machine Learning Repository.  DOI: 10.24432/C5859H.
  * Tsanas, 2012 Tsanas, Athanasios & Xifara, A. (2012).  Energy efficiency.  UCI Machine Learning Repository.  DOI: 10.24432/C51307.
  * Yeh, 2007 Yeh, I.-C. (2007).  Concrete Compressive Strength.  UCI Machine Learning Repository.  DOI: 10.24432/C5PK67.