# Real-Time Tube-Based Non-Gaussian Risk Bounded Motion Planning for
Stochastic Nonlinear Systems in Uncertain Environments via Motion Primitives

Weiqiao Han∗1, Ashkan Jasour∗2, Brian Williams1 1Computer Science and
Artificial Intelligence Laboratory, Massachusetts Institute of Technology,
<EMAIL_ADDRESS>2 Massachusetts Institute of Technology,
<EMAIL_ADDRESS>Now with NASA Jet Propulsion Lab, California Institute of
Technology. ∗These authors contributed equally to the paper.

###### Abstract

We consider the motion planning problem for stochastic nonlinear systems in
uncertain environments. More precisely, in this problem the robot has
stochastic nonlinear dynamics and uncertain initial locations, and the
environment contains multiple dynamic uncertain obstacles. Obstacles can be of
arbitrary shape, can deform, and can move. All uncertainties do not
necessarily have Gaussian distribution. This general setting has been
considered and solved in [1]. In addition to the assumptions above, in this
paper, we consider long-term tasks, where the planning method in [1] would
fail, as the uncertainty of the system states grows too large over a long time
horizon. Unlike [1], we present a real-time online motion planning algorithm.
We build discrete-time motion primitives and their corresponding continuous-
time tubes offline, so that almost all system states of each motion primitive
are guaranteed to stay inside the corresponding tube. We convert probabilistic
safety constraints into a set of deterministic constraints called risk
contours. During online execution, we verify the safety of the tubes against
deterministic risk contours using sum-of-squares (SOS) programming. The
provided SOS-based method verifies the safety of the tube in the presence of
uncertain obstacles without the need for uncertainty samples and time
discretization in real-time. By bounding the probability the system states
staying inside the tube and bounding the probability of the tube colliding
with obstacles, our approach guarantees bounded probability of system states
colliding with obstacles. We demonstrate our approach on several long-term
robotics tasks.

##  I Introduction

When a robot moves around in the real world, it encounters all kinds of
stochasticity from nature. The uneven terrain, the volatile weather, and
surrounding uncertain moving obstacles all could pose challenges to the robot
navigation. How to plan motions for the robot to navigate safely and reach the
goal is a long-standing problem.

Traditional planning algorithms, such as rapidly exploring random tree (RRT),
probabilistic roadmap (PRM), and virtual potential field methods, plans paths
in deterministic environments. Planning algorithms assuming stochasticity and
uncertainty are mainly limited to Gaussian uncertainty and convex obstacles
[2, 3, 4, 5, 6, 7, 8]. There are also sampling based methods that do not
necessarily assume Gaussian uncertainty [9, 10, 11, 12].

[1] considered the most general setting so far, where (a) the system dynamics
is stochastic and nonlinear, and the uncertainty is not necessarily Gaussian;
(b) the initial position of the system is uncertain and is necessarily
Gaussian; (c) the obstacles can be of arbitrary shape, can deform, can move,
and have arbitrary uncertainty, not necessarily Gaussian. [1] proposed a
trajectory optimization method to plan a path with bounded risk, where risk is
defined to be the probability of colliding with obstacles or not reaching the
goal. However, since the planned path is open loop, the uncertainty grows as
the system evolves. Therefore, the planning horizon is limited. We want to
plan longer horizons for the robot to accomplish more complex tasks.

Figure 1: Motion primitives and tubes

One would think of using closed loop control to track the nominal trajectory
planned by the trajectory optimization. However, currently closed loop
controllers are usually designed for linear systems. For example, in [13],
controllers are designed to steer the stochastic time-varying linear system to
the goal with desired covariance, while satisfying chance constraints. It is
not obvious how such controller would behave on a linearized stochastic
nonlinear system. It is also unlikely to incorporate the design of such closed
loop controller into the trajectory optimization for long-term motion
planning.

One approach to long-term tasks is online planning with motion primitives [14,
15, 16, 17, 18, 19]. For example, [14] considers state space sampling of
motion primitives, and then use numerical methods to solve for the controller.
They assume the system dynamics is deterministic. In contrast, we are going to
use motion primitives for stochastic systems.

In this paper, we consider the motion planning problem in the same setting as
in [1], satisfying conditions (a)–(c) listed above, and in addition, we
consider (d) long-term tasks, where the short-term planning method in [1]
would fail, as the uncertainty of the system states grows too large over a
long time horizon. We present a real-time online motion planning approach to
this general problem. We first generate motion primitives in control space
over a discrete time horizon, and then build a continuous-time tube in state
space for each motion primitive, such that almost all system states remain in
the tube over the discrete time horizon. We provide theoretical guarantees
that the probability of system states outside of the tube is bounded above by
a very small number. For the uncertain environment, we use concentration
inequalities to convert the uncertain obstacles into deterministic risk
contours. During online execution, we verify the safety of the tubes against
deterministic risk contours using sum-of-squares (SOS) programming. Among all
feasible tubes that stay inside the low-risk contour, an objective is used to
score their corresponding motion primitives, and the one with the highest
score is selected and executed. Our approach bridges the gap between the low-
level discrete-time control sequence planning of stochastic nonlinear systems,
and the high-level continuous-time state space planning, where safety can be
verified via SOS programming. Our approach provides theoretical guarantees
that the probability of system states colliding with any obstacle is bounded
above by a small number. We demonstrate our approach on several long-term
robotics tasks.

##  II Problem Definition

Consider the stochastic nonlinear dynamical system represented by

$\mathbf{x}_{t+1}=f(\mathbf{x}_{t},\mathbf{u}_{t},\mathbf{\omega}_{t})$ (1)

where $\mathbf{x}_{t}\in\mathcal{X},\mathbf{u}_{t}\in\mathcal{U}$ are the
state and control input at time $t$, respectively,
$\mathcal{X}\subseteq\mathbb{R}^{n}$ is the state space,
$\mathcal{U}\subseteq\mathbb{R}^{m}$ is the control space, $n,m\in\mathbb{N}$,
$\mathbf{\omega}_{t}$ is the noise at time $t$ and it may not necessarily have
Gaussian distribution. The initial state $\mathbf{x}_{0}$ can either be
deterministic, or can be a random variable with some known probability
distribution, not necessarily Gaussian. In the environment, there are multiple
obstacles, denoted by $\mathcal{O}_{i},i=1,\ldots,M$, where
$M\in\mathbb{N}^{+}$, and each obstacle $\mathcal{O}_{i}$ has a polynomial
representation

$\displaystyle\mathcal{O}_{i}(\tilde{w}_{i},t)=\\{\mathbf{x}\in\mathcal{X}:\
p_{i}(\mathbf{x},\tilde{\omega}_{i},t)\leq 0\\},\ i=1,...,M$ (2)

where $p_{i}$ is a polynomial and $\tilde{\omega}_{i}$ is a random variable
with some known probability distribution, not necessarily Gaussian. The
obstacle can change its position and shape over time, and hence the time
variable $t$ in the polynomial representation. We define the risk to be the
probability of collision with any uncertain obstacles at any time step. The
goal is to plan control sequences $\mathbf{u}_{0},\mathbf{u}_{1},\ldots$
online to steer the system into the goal region
$\mathbf{x}_{goal}\subseteq\mathcal{X}$. More precisely, we want to solve the
following planning problem:

Risk Bounded Motion Planning Problem

$\displaystyle\min_{\mathbf{u}}$ $\displaystyle
E[l_{f}(\mathbf{x}_{T})+\sum_{t=0}^{T-1}l(\mathbf{x}_{t},\mathbf{u}_{t},\omega_{t})]$
(3) s.t.
$\displaystyle\mathbf{x}_{t+1}=f(\mathbf{x}_{t},\mathbf{u}_{t},\omega_{t}),\ \
|_{t=0}^{T-1},$
$\displaystyle\text{Prob}(\mathbf{x}_{t}\in\mathcal{O}_{i}(\tilde{\omega}_{i},t))\leq\Delta_{o},\
\ |_{t=0}^{T-1},\ |_{i=1}^{M}$
$\displaystyle\text{Prob}(\mathbf{x}_{T}\not\in\mathbf{x}_{goal})\leq\Delta_{goal},$
$\displaystyle\mathbf{x}_{0}\sim pr(\mathbf{x}_{0})$

where $\Delta_{o},\Delta_{goal}\in[0,1]$ are the given acceptable risk levels,
and $pr(\mathbf{x}_{0})$ is the given probability distribution of the initial
system states. We assume that either there is a global planner, obtained from
the method in [20], for example, which roughly traces out a general path from
the initial position to the goal region, or there is a high-level objective
function that guides the system to the goal region.

###  II-A Notations and Definitions

Let $\mathbb{R}[x]$ be the polynomial ring in the variables
$\mathbf{x}=(x_{1},\ldots,x_{n})$ with real coefficients. A polynomial
$p(\mathbf{x})\in\mathbb{R}[x]$ can be written as
$p(\mathbf{x})=\sum_{\bm{\alpha}\in\mathbb{N}^{n}}p_{\bm{\alpha}}\mathbf{x}^{\bm{\alpha}}$,
where $\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{n})\in\mathbb{N}^{n}$ and
$\mathbf{x}^{\bm{\alpha}}=\prod_{i=1}^{n}x_{i}^{\alpha_{i}}$ is a monomial in
standard basis.

Let $(\Omega,\Sigma,\mu)$ be a probability space, where $\Omega$ is the sample
space, $\Sigma$ is the $\sigma$-algebra of $\Omega$, and $\mu:\Sigma\to[0,1]$
is the probability measure on $\Sigma$. Suppose
$\mathbf{x}\in\Omega\subseteq\mathbb{R}^{n}$ is an $n$-dimensional random
vector. Let $\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{n})\in\mathbb{N}^{n}$. The
expectation of $\mathbf{x}^{\bm{\alpha}}$ defined as
$m_{\bm{\alpha}}=E[\mathbf{x}^{\bm{\alpha}}]$ is a moment of order $\alpha$,
where $\alpha=\sum_{i}\alpha_{i}$. The sequence of all moments of order
$\alpha$, denoted by $\mathbf{m}_{\alpha}$, is the expectation of all
monomials of order $\alpha$ sorted in graded reverse lexicographic order
(grevlex).

Sum-of-Squares Polynomial: Polynomial $\mathcal{P}(x)$ is a sum-of-squares
polynomial if it can be written as a sum of finitely many squared polynomials,
i.e., $\mathcal{P}(x)=\sum_{j=1}^{m}h_{j}(x)^{2}$ for $h_{j}(x)$ for $1\leq
j\leq m$. Checking if a polynomial is SOS can be cast into a semidefinite
programming program. One can use the packages like Yalmip [21] and Spotless
[22] to check if a polynomial is SOS.

##  III Method

In this section, we first review how to convert uncertain obstacles in the
environment into risk contours [23]. Second, we review how to verify the
safety of the tube using risk contours and SOS programming [23]. Next we
present how to generate motion primitives and their corresponding tubes, and
discuss several points to optimize implementation in online execution.
Finally, we prove the theoretical guarantees on the bounded risk.

###  III-A Risk Contour Representation of Uncertain Environments

In this subsection, we review risk contours [24], which have been used in [1,
20, 23].

Suppose we want to bound the risk of colliding with the obstacle $\mathcal{O}$
defined in (2) by $\Delta$, i.e.,
$\text{Prob}(\textbf{x}\in\mathcal{O})=\text{Prob}(p(\mathbf{x},\tilde{\omega},t)\leq
0)\leq\Delta$. Then, the associated risk contour is defined as

$\mathcal{C}_{r}^{\Delta}(t)=\\{\mathbf{x}\in\mathcal{X}:\hbox{Prob}\left(p(\mathbf{x},\tilde{\omega},t)\leq
0\right)\leq\Delta\\}$ (4)

The following result holds true.

  * Theorem(Theorem 1 in [23]).

The set

$\hat{\mathcal{C}}_{r}^{\Delta}(t)=\left\\{\mathbf{x}\in\mathcal{X}:\frac{E[g^{2}]-E[g]^{2}}{E[g^{2}]}\leq\Delta,E[g]\leq
0,\right\\}$ (5)

where $g=-p(\mathbf{x},\tilde{\omega},t)$, is the inner approximations of the
risk contour in (4).

In the proof of the Theorem above we use concentration inequalities to provide
bounds on random variables [23]. More precisely, by Cantelli’s inequality, for
any random variable $z$,

$\displaystyle\text{Prob}(z\geq 0)\leq\frac{E[z^{2}]-E[z]^{2}}{E[z^{2}]},$ (6)

whenever $E[z]\leq 0$. Any upper bound on the right-hand side of (6) would be
an upper bound of the probability on the left-hand size.

The set $\hat{\mathcal{C}}_{r}^{\Delta}(t)$ is a region in state space where
the probability of collision with the uncertain obstacle is $\leq\Delta$. For
one obstacle, $\hat{\mathcal{C}}_{r}^{\Delta}(t)$ is defined by 2 polynomials.
In general, for $M$ obstacles, the risk-bounded set
$\hat{\mathcal{C}}_{r}^{\Delta}(t)$ is defined by $2M$ polynomials, $E[g_{i}]$
and $E[g_{i}^{2}]$, where $g_{i}=-p_{i}(\mathbf{x},\tilde{\omega},t)$,
$i=1,\ldots,M$, and $\hat{\mathcal{C}}_{r}^{\Delta}(t)$ is the region in state
space where the probability of collision with any of the $2M$ uncertain
obstacles is $\leq\Delta$.

###  III-B Tube Verification

Here we briefly review tube safety verification. The provided method verifies
the safety of the tube in the presence of uncertain obstacles without the need
for uncertainty samples and time discretization. Curious readers are suggested
to refer to [23] for more details. Given a polynomial trajectory
$\mathcal{P}(t)$ parameterized by $t$, a tube around it takes the form

$\mathcal{Q}(\mathcal{P}(t))=\\{\mathbf{x}\in\mathbb{R}^{n}:(\mathbf{x}-\mathcal{P}(t))^{T}{Q}(\mathbf{x}-\mathcal{P}(t))\leq
1\\}$ (7)

where $Q\in\mathbb{R}^{n\times n}$ is the given positive definite matrix and
$t\in[t_{0},t_{f}]$.

We can use SOS programming to verify if a tube is inside the risk-bounded set
$\hat{\mathcal{C}}_{r}^{\Delta}(t)$ in (5) over the entire planning time
horizon $t\in[t_{0},t_{f}]$. This is the following theorem.

  * Theorem(Theorem 3 in [23]).

Denote $P_{1i}(\textbf{x},t)={E}[g_{i}^{2}]$ and
$P_{2i}(\textbf{x},t)={E}[g_{i}]$. The given tube
$\mathcal{Q}(\mathcal{P}(t))$ over $t\in[t_{0},t_{f}]$ is inside the risk-
bounded set $\hat{\mathcal{C}}_{r}^{\Delta}(t)$ in (5) if the polynomials of
$\hat{\mathcal{C}}_{r}^{\Delta}(t)$ take the following SOS representation:

$\displaystyle
P^{2}_{2i}\left(\mathcal{P}(t)+\hat{\mathbf{x}}_{0},t\right)-(1-\Delta)P_{1i}(\mathcal{P}(t)+\hat{\mathbf{x}}_{0},t)=$
(8)
$\displaystyle\quad\quad\quad\quad{\sigma_{0}}_{i}(t,\hat{\mathbf{x}}_{0})+{\sigma_{1}}_{i}(t,\hat{\mathbf{x}}_{0})(t-t_{0})(t_{f}-t)$
$\displaystyle\quad\quad\quad\quad+{\sigma_{2}}_{i}(t,\hat{\mathbf{x}}_{0})(1-\hat{\mathbf{x}}_{0}^{T}{Q}\hat{\mathbf{x}}_{0})|_{i=1}^{M}$
$\displaystyle P_{2i}(\mathcal{P}(t)+\hat{\mathbf{x}}_{0},t)=$ (9)
$\displaystyle\quad\quad\quad\quad{\sigma_{3}}_{i}(t,\hat{\mathbf{x}}_{0})+{\sigma_{4}}_{i}(t,\hat{\mathbf{x}}_{0})(t-t_{0})(t_{f}-t)$
$\displaystyle\quad\quad\quad\quad+{\sigma_{5}}_{i}(t,\hat{\mathbf{x}}_{0})(1-\hat{\mathbf{x}}_{0}^{T}{Q}\hat{\mathbf{x}}_{0})|_{i=1}^{M}$

where $\hat{\mathbf{x}}_{0}\in\mathbb{R}^{n}$ is the variable vector,
${\sigma_{j}}_{i}(t,\hat{\mathbf{x}}_{0}),j=0,...,5,i=1,...,M$ are SOS
polynomials with appropriate degrees.

As shown in the experiment section, one can verify the obtained safety SOS
conditions in real-time. Note that the complexity of the provided safety SOS
conditions is independent of the size of the planning time horizon
$[t_{0},t_{f}]$ and the length of the polynomial trajectory $\mathcal{P}(t)$.
Hence, they can be easily used to verify the safety of trajectories and their
neighborhoods represented by tubes in uncertain environments over the long
planning time horizon.

###  III-C Motion Primitives and Their Corresponding Tubes

A motion primitive is defined to be a pre-computed sequence of control inputs
$\\{\mathbf{u}_{0},\ldots,\mathbf{u}_{T-1}\\}$ of length $T\in\mathbb{N}$.
When working with stochastic systems, each control $\mathbf{u}_{i}$ can depend
on the state random variable $\mathbf{x}_{i}$, i.e.,
$\mathbf{u}_{i}=\mathbf{u}_{i}(\mathbf{x}_{i})$. In general, there are two
ways to generate motion primitives for deterministic systems [14]. One way is
to sample the state space. In this way, the state constraints imposed by the
environment, such as obstacle avoiding and goal reaching, can be easily
satisfied. However, it is hard to satisfy dynamics constraints imposed by the
underlying dynamical system. The other way is to sample the control space. The
states are formed by simulating forward in dynamics. The state constraints in
this case are not easily satisfied. In light of this, we consider two ways to
generate motion primitives and their associated tubes. One is based on control
space sampling and the other is based on state space sampling.

####  III-C1 Method 1

Control space sampling. In this method, we first sample control sequences from
the control space to get a desired motion primitive. Next, we fit a polynomial
nominal trajectory $\bar{\textbf{x}}_{\theta}(t)$, $t\in[0,1]$, parameterized
by $\theta$ in some parameter space $\Theta$, and fit a polynomial tube to the
motion primitive. Here, fitting a polynomial nominal trajectory rigorously
means solving the following optimization problem

$\displaystyle\min_{\theta\in\Theta}$
$\displaystyle\sum_{k=1}^{T}|E[\mathbf{x}_{k}]-\bar{\textbf{x}}_{\theta}(k/T)|^{2}$
subject to $\displaystyle\quad\text{system dynamics }(\ref{eq:sysdyn})$
$\displaystyle\quad\text{initial state distribution }\mathbf{x}_{0}\sim
p(\mathbf{x}_{0})$

If $E[\mathbf{x}_{k}]$, $k=1,\ldots,T$, can be calculated analytically or
estimated using sampling, then the problem amounts to minimizing a polynomial
function of $\theta$. Besides solving the optimization problem, we can also
use sampling to approximately find a good enough $\theta$. The construction of
the tube is to be explained in the next subsection.

####  III-C2 Method 2

State space sampling. In this method, we first generate a desired polynomial
nominal trajectory $\bar{\textbf{x}}(t)$, $t\in[0,1]$, for some
$\theta\in\Theta$. Next we apply trajectory optimization or sampling to look
for the control sequence to follow the nominal trajectory so that the expected
states are on or near the nominal trajectory. The trajectory optimization has
the form of

$\displaystyle\min_{\mathbf{u}_{0},\ldots,\mathbf{u}_{T-1}}$
$\displaystyle\sum_{k=1}^{T}|E[\mathbf{x}_{k}]-\bar{\textbf{x}}_{\theta}(k/T)|^{2}$
subject to $\displaystyle\quad\text{system dynamics }(\ref{eq:sysdyn})$
$\displaystyle\quad\text{initial state distribution }\mathbf{x}_{0}\sim
p(\mathbf{x}_{0})$

If this optimization problem can be solved, we get a sequence of control
inputs $\\{\mathbf{u}_{0},\ldots,\mathbf{u}_{T-1}\\}$, which is the motion
primitive. Besides solving the optimization problem, we can also use sampling
to approximately find a good enough control sequence
$\\{\mathbf{u}_{0},\ldots,\mathbf{u}_{T-1}\\}$. Finally, we fit the trajectory
with a tube, which is to be explained in the following subsection.

Figure 2: Motion primitives and their corresponding tubes. The top two plots
are motion primitives and tubes for the underwater vehicle model. The bottom
two plots are motion primitives and tubes for the ground vehicle model. In
each plot, 10,000 samples are simulated forward using each motion primitive.
Each black circle surrounding particles is the corresponding tube at that
particular time step. In the left two plots, the tube sizes are computed using
sampling, and in the right two plots, the tube sizes are computed using
analytical methods. Although the tubes on the right are larger in size, they
provide theoretical guarantees that at least 99.9% of the system states stay
inside the tube.

###  III-D Tube Construction

A tube along the polynomial nominal trajectory $\bar{\textbf{x}}(t)$ in the
form of

$\displaystyle\mathcal{Q}_{\phi}(\bar{\textbf{x}}(t))=\\{\mathbf{x}\in\mathbb{R}^{n_{x}}:(\mathbf{x}-\bar{\textbf{x}}(t))^{T}{Q}_{\phi}(\mathbf{x}-\bar{\textbf{x}}(t))\leq
1\\}$ (10)

for $t\in[0,1]$ is constructed so that the probability of system states
staying inside the tube is $\geq(1-\Delta_{tube})$, where $\Delta_{tube}>0$ is
a small positive number, i.e.,

$\displaystyle\text{Prob}(\textbf{x}_{t}\in\mathcal{Q}_{\phi}(\bar{\textbf{x}}(t)))\geq
1-\Delta_{tube}$ (11)

Here the parameter to be determined is $\phi\in\Phi$, which controls the size
of the tube. For example, in 2D the tube at time $t$ can be a disc

$\displaystyle\mathcal{Q}_{\phi}(\bar{\textbf{x}}(t))=\\{\mathbf{x}\in\mathbb{R}^{2}:|(\mathbf{x}-\bar{\textbf{x}}(t))|^{2}\leq
r^{2}\\}$ (12)

and the parameter is $\phi=r^{-1}$ for ${Q}_{\phi}=\phi*I_{2}$, where $I_{2}$
is the $2\times 2$ identity matrix.

Note that tube is continuous-time, while motion primitives are discrete-time.
The discrete time steps in motion primitives correspond to $t=0,{1\over
T},\ldots,{T-1\over T},1$, in the interval $t\in[0,1]$. We want to ensure that
Inequality (11) holds for any $t=0,{1\over T},\ldots,{T-1\over T},1$. Clearly
a larger tube size would allow more system states to stay inside the tube, but
it would be more conservative for high level state space planning when we
check the safety of the tube against obstacles. So we want the tube size to be
tight. For certain parameterizations of tubes, such as discs as in Equation
12, we can use binary search to determine the min-sized tube, so that a little
smaller in size would break Inequality (11). In order for binary search to
work, we need a procedure to verify if the tube of size $\phi$ satisfies
Inequality (11) or not. There are two ways – one is sampling, and the other is
an analytical method, which provides theoretical guarantees.

####  III-D1 Method 1

Sampling. We sample $N$ states, and simulate forward using stochastic
dynamics. Given a $\phi$, at each time step $k$ over the horizon
$0,1,\ldots,T$, the number of samples inside the tube
$\mathcal{Q}_{\phi}(\bar{\textbf{x}}(k/T))$ is $N_{in}(k)$. If
$N_{in}(k)/N<1-\Delta_{tube}$ for some $k$, then the parameter $\phi$ is
rejected. Otherwise $\phi$ is feasible. Different from the analytical method,
the sampling method does not provide any theoretical guarantees, but it
usually generates smaller sized tubes than the analytical method Figure 2.

####  III-D2 Method 2

Analytical method. The analytical method applies to systems whose future state
moments can be calculated analytically. This class of systems includes most
robotic systems, such as vehicles and drones, as well as the examples in [1,
25] and examples in this paper, with state-independent control inputs or
certain state-dependent controls, which are to be shown in the following
examples. Given a $\phi$, at each time step $k$ over the horizon
$0,1,\ldots,T$, we want

$\displaystyle\text{Prob}(\textbf{x}_{k}\not\in\mathcal{Q}_{\phi}(\bar{\textbf{x}}(k/T)))\leq\Delta_{tube}$
(13)

i.e.,

$\displaystyle\text{Prob}((\textbf{x}_{k}-\bar{\textbf{x}}(t))^{T}{Q}_{\phi}(\textbf{x}_{k}-\bar{\textbf{x}}(t))\geq
1,t=k/T)\leq\Delta_{tube}$ (14)

The probability on the left-hand side of (14) can be relaxed using
concentration inequalities, such as Cantelli’s inequality (6), and the
probability inequality (14) can be converted to an inequality involving
moments of $\textbf{x}_{k}$, as in Section III-A. If the moment inequalities
are satisfied for all $k$, then $\phi$ is feasible, otherwise $\phi$ is
rejected. In contrast to sampling, the analytical method find the min-sized
tube that guarantees Inequality (11).

Example 1. Consider the underwater robot model given by

$\displaystyle x_{t+1}=x_{t}+\Delta
T(v_{t}+{\omega_{v}}_{t})\cos(\theta_{t}+{\omega_{\theta}}_{t})$ (15)
$\displaystyle y_{t+1}=y_{t}+\Delta
T(v_{t}+{\omega_{v}}_{t})\sin(\theta_{t}+{\omega_{\theta}}_{t})$

where $(x_{t},y_{t})$ is the 2D position at time $t$, and the discrete time
interval $\Delta T=0.1$. The noise terms ${\omega_{v}}_{t}$ and
${\omega_{\theta}}_{t}$ have uniform distribution on $[-0.1,0.1]$ at any time
$t$. We generate 5 motion primitives with horizon $T=5$, and fit each of them
with a quadratic function of the form

$\displaystyle\bar{\textbf{x}}(t)=(x(t),y(t))$ (16) $\displaystyle
x(t)=\theta_{1}t$ $\displaystyle y(t)=\theta_{2}x(t)^{2}$

for $t\in[0,1]$, where $\theta=(\theta_{1},\theta_{2})$ is the parameter. We
choose tube to be discs as in Equation 12, and we use binary search to find
the min-sized tube. The control inputs of the motion primitives are state
independent and hence the analytical method applies. We use both sampling
(Figure 2 top left) and the analytical method (Figure 2 top right) to
construct the tube for each motion primitive. The analytical method guarantees
that at least $99.9\%$ of the system states stay inside the tube.

Example 2. Consider the ground vehicle model whose dynamics is in the form of

$\displaystyle x_{t+1}=x_{t}+\Delta Tv_{t}\cos(\theta_{t})$ (17)
$\displaystyle y_{t+1}=y_{t}+\Delta Tv_{t}\sin(\theta_{t})$ $\displaystyle
v_{t+1}=v_{t}+\Delta T(a_{t}+{\omega_{v_{t}}})$
$\displaystyle\theta_{t+1}=\theta_{t}+\Delta T(u_{t}+\omega_{\theta_{t}})$

where the state is $(x_{t},y_{t},v_{t},\theta_{t})$, $(x_{t},y_{t})$ is the 2D
position at time $t$, $v_{t}$ is the velocity at time $t$, and $\theta_{t}$ is
the angle at time $t$. The discrete time interval is $\Delta T=0.1$. The noise
$\omega_{v_{t}}$ has normal distribution with mean 0 and variance $0.09$,
while $\omega_{\theta_{t}}\sim 3B$, where $B$ has Beta distribution with
parameters $(1,3)$ over $[0,1]$. The initial distribution of $x_{0}$ and
$y_{0}$ are both normal distribution with mean 0 and variance $0.01^{2}$,
while the initial deterministic states of $v_{0}$ and $\theta_{0}$ are $1$ and
0, respectively. Similar to Example 1, we generate 5 motion primitives with
horizon $T=5$, and fit each primitive with a quadratic nominal trajectory in
the form of Equation 16. The tubes are discs as in Equation 12. The control
inputs of the motion primitives are state dependent, and they are to track
certain nominal velocities $\bar{v}_{t}$ and nominal angles
$\bar{\theta}_{t}$, i.e., the control inputs are
$a_{t}=(\bar{v}_{t+1}-v_{t})/\Delta T$,
$u_{t}=(\bar{\theta}_{t+1}-\theta_{t})/\Delta T$. The dynamics for $v$ and
$\theta$ essentially becomes $v_{t+1}=\bar{v}_{t+1}+\Delta T{\omega_{v_{t}}}$,
and $\theta_{t+1}=\bar{\theta}_{t+1}+\Delta T\omega_{\theta_{t}}$. For this
system, the moments of system states can be calculated analytically. We use
both sampling (Figure 2 bottom left) and the analytical method (Figure 2
bottom right) to construct the tube for each motion primitive. The analytical
method guarantees that at least $99.9\%$ of the system states stay inside the
tube.

###  III-E Online Execution

So far we have constructed various motion primitives and their corresponding
tubes offline. During online execution, the system, based on the current
state, selects a motion primitive that is risk bounded, verified via SOS
programming in Section III-B, with respect to its surrounding obstacles. The
system executes the first one or few steps of the motion primitive and re-
plans based on the new system state, much like model predictive control (MPC)
style planning. One planning cycle consists of selecting a risk-bounded motion
primitive based on the current state and executing the first one or few steps
of the motion primitive.

The most naive implementation for selecting a risk-bounded motion primitive is
to verify all tubes against all obstacles, and select one of the tubes that
are risk bounded with respect to all obstacles, and execute its corresponding
motion primitive. Since online computation resources are limited, we can
optimize the implementation in a number of ways:

  * •

If the motion primitives are ordered from left to right, viewed in the
direction of the current velocity, then binary search algorithm can be used to
look for feasible tubes against a particular obstacle.

  * •

We only check the safety of the tubes against obstacles within certain
distance of the system. Obstacles far away can clearly be ignored.

  * •

Use lower order polynomials, such as quadratics, circles, and ellipses, to
represent the obstacles would make verification faster, though they might be
more conservative.

Finally, instead of picking feasible motion primitives randomly, we can use an
objective to rank all feasible motion primitives and the one with the highest
score will be selected.

###  III-F Theoretical Guarantees and Summary

Theoretical Guarantees of Bounded Risk. During offline planning, we generate
motion primitives, fit polynomial nominal trajectories, and build tight tubes
so that if the tube starts from the current state, then the probability of
system states staying inside the tube $\mathcal{Q}(\bar{\textbf{x}}(t))$ with
$t\in[0,1]$, for any $t$ in the planning horizon of the motion primitive, is
at least $1-\Delta_{tube}$, i.e., for $t=0,{1\over T},\ldots,1$,

$\displaystyle\text{Prob}(\textbf{x}_{t}\in\mathcal{Q}(\bar{\textbf{x}}(t)))\geq
1-\Delta_{tube}.$ (18)

We convert the uncertain environment into risk contours so that we can check
if a tube $\mathcal{Q}(\bar{\textbf{x}}(t))$ is in the risk-bounded set
$\hat{\mathcal{C}}_{r}^{\Delta}(t))$ using SOS programming. If a tube
$\mathcal{Q}(\bar{\textbf{x}}(t))$ is inside the risk-bounded set
$\hat{\mathcal{C}}_{r}^{\Delta}(t))$, then the probability of the system
states inside the tube colliding with any obstacle $\mathcal{O}$ is bounded by
$\leq\Delta_{o}$, i.e.,

$\displaystyle\text{Prob}(\textbf{x}_{t}\in\mathcal{O}|\textbf{x}_{t}\in\mathcal{Q}(\bar{\textbf{x}}(t))\subseteq\hat{\mathcal{C}}_{r}^{\Delta}(t))\leq\Delta_{o}.$
(19)

During online execution, we select a motion primitive whose tube is inside the
risk-bounded set $\hat{\mathcal{C}}_{r}^{\Delta}(t)$. Therefore, the
probability of system states colliding with any obstacle is bounded by
$\Delta_{o}+\Delta_{tube}$, since

$\displaystyle\quad\text{Prob}(\textbf{x}_{t}\in\mathcal{O}|\mathcal{Q}(\bar{\textbf{x}}(t))\subseteq\hat{\mathcal{C}}_{r}^{\Delta}(t))$
$\displaystyle\leq\text{Prob}(\textbf{x}_{t}\in\mathcal{O}|\textbf{x}_{t}\in\mathcal{Q}(\bar{\textbf{x}}(t))\subseteq\hat{\mathcal{C}}_{r}^{\Delta}(t))\text{Prob}(\textbf{x}_{t}\in\mathcal{Q}(\bar{\textbf{x}}(t)))$
$\displaystyle\quad+\text{Prob}(\textbf{x}_{t}\not\in\mathcal{Q}(\bar{\textbf{x}}(t)))\quad(\text{from
Law of Total Probability})$
$\displaystyle\leq{\Delta_{o}+\Delta_{tube}}\quad\quad\quad\quad\ \
(\text{from
}(\ref{eq:theoretical_guarantee_eq1}),(\ref{eq:theoretical_guarantee_eq2}))$
(20)

The bound in (III-F) is only for one tube, or one planning cycle. Suppose the
system reaches the goal in $N$ planning cycles. Then there are $N$ tubes
connected consecutively, paving the path from the initial position to the goal
region. Let $A_{N}$ denote the event that the system stays in the tubes over
all the $N$ planning cycles. Let $B_{N}$ denote the event that the system
collides with any obstacle over all the $N$ planning cycles. Then

$\displaystyle\text{Prob}(A_{N})=\text{Prob}(\textbf{x}_{t}\in\mathcal{Q}(\bar{\textbf{x}}(t)))^{N}\geq(1-\Delta_{tube})^{N}.$
(21)

By the same reasoning as in (III-F),

$\displaystyle\text{Prob}(B_{N})$
$\displaystyle\leq\Delta_{o}+(1-(1-\Delta_{tube})^{N})$ (22)
$\displaystyle\leq\Delta_{o}+N\Delta_{tube}$ (23)

Both (22) and (23) are bounds on the probability of colliding with any
obstacle over the entire $N$ planning cycle. When $\Delta_{tube}$ is very
small, e.g., $0.001$ in our examples, both bounds are close to each other.

User-side Risk Allocation. Suppose $\Delta$ is the total risk the user
considers. The user can first estimate an upper bound $M$ of the total
planning cycle, and then choose $\Delta_{o}$ and $\Delta_{tube}$ so that

$\displaystyle\Delta_{o}+M\Delta_{tube}\leq\Delta.$ (24)

Suppose there are $N$ actual planning cycles, where $N\leq M$. Then the risk
is bounded by $\Delta$, i.e.,

$\displaystyle\text{Prob}(B_{N})\leq\Delta_{o}+N\Delta_{tube}\leq\Delta_{o}+M\Delta_{tube}\leq\Delta.$
(25)

Algorithm Summary. We summarize the algorithm in Algorithm 1. Note that Line 5
says that we generate risk contours of the uncertain environment online at
each time step. This can vary depending on the robot platform and the
environment. If the environment is static, then risk contours can be generated
offline. If the environment is highly dynamic, such as a traffic scene, and if
the robot is an autonomous driving car with strong computing power, then the
prediction module of the car predicts the future trajectories of surrounding
vehicles at high frequency [26, 27] and risk contours can be form at high
frequency, too.

Offline:

1 Given a total risk $\Delta$, estimate an upper bound $M$ of the total
planning cycle, and choose $\Delta_{o}$ and $\Delta_{tube}$ satisfying
Equation 24.

2 Design motion primitives (Section III-C).

3 Build tubes for the motion primitives (Section III-D).

Online:

4 For each time step until goal is reached:

5 (1) Generate risk contours of the uncertain environment (Section III-A).

6 (2) Check tubes against risk contours (Section III-B).

7 (3) Rank the tubes that are safe using a user-defined objective, and choose
the one with the highest score (Section III-E).

8 (4) Execute the first one or few time steps of the chosen motion primitive
(Section III-E).

Algorithm 1 Real-Time Tube Based Non-Gaussian Risk Bounded Motion Planning

##  IV Experiments

Figure 3: Underwater Robot in Clustered Environment. Left: Our method. Right:
CC-RRT. The red curves are boundaries of obstacles’ risk contours, which are
4th order polynomials. The blue piecewise-linear trajectory is the general
path given by the high-level planner. The green circle represents the goal
region. The plot on the left shows two trajectories given by our method,
colored in magenta and cyan, starting from two different initial positions,
marked by two black stars, and reaching the goal in the end. The plot on the
right shows a cyan trajectory given by CC-RRT. Figure 4: Histogram of the
number of planning cycles to reach the goal from 100 initial states.

###  IV-A Underwater Robot in Clustered Environment

In this example, we consider the underwater vehicle whose dynamics is given by
Equation 15. In the environment there are multiple static obstacles, and we
approximate them, in the spirit of Section III-E, in the form of

$\displaystyle\\{(x,y)\in\mathbb{R}^{2}:(x-x_{i})^{2}+(y-y_{i})^{2}\leq
w^{2}\\}$

where $(x_{i},y_{i})$ is the center for the $i$-th obstacle, and $w$ is a
random variable with uniform distribution on $[0.3,0.4]$. The risk contours of
the obstacles are 4th order polynomials in 2 variables and we plotted only
their boundaries in (Figure 3). We want the total risk level to be
$\Delta=0.2$. We estimate the total number of planning cycles is bounded by
100. So we choose risk contours with the risk level of $\Delta_{o}=0.1$, and
choose $\Delta_{tube}=0.001$. We design motion primitives using the analytical
method. The goal region is given by
$\\{(x,y)\in\mathbb{R}^{2}:(x-x_{g})^{2}+(y-y_{g})^{2}\leq r_{g}^{2}\\}$,
where $(x_{g},y_{g})=(5.5,2)$ and $r_{g}=0.09$. The high-level planner
provides a general path represented by the blue piecewise linear trajectory in
Figure 3. The objective to rank the feasible motion primitives is the sum of
squared distance between expected future states of the motion primitive and
the path given by the high-level planner. The feasible motion primitive that
minimizes the objective is selected.

We assume the initial state has uniform distribution on the square with side
length 1 centered at the point $(0,3)$. We sample 100 initial states and apply
our method to those samples. The system re-plans every 2 time steps. All
samples reach the goal while staying inside the tube. The success rate is
100%. We plotted the histogram of the number of planning cycles to reach the
goal in Figure 4. The max number of planning cycles is 37, which is less than
our estimate of 100. This verifies that our estimate is indeed a good upper
bound.

In Figure 3 Left, we plotted the trajectory starting from two initial
positions, $(0,3)$ and $(-0.5,4)$, both with initial velocity $(1,0)$. The
initial state $(-0.5,4)$ is actually outside of the initial distribution. The
trajectory starting from $(0,3)$ finishes in $N_{1}=34$ planning cycles, and
hence the risk is bounded by $\Delta_{o}+N_{1}\Delta_{tube}=0.134$, or
$\Delta_{o}+(1-(1-\Delta_{tube})^{N_{1}})=0.1335$. The trajectory starting
from $(-0.5,4)$ finishes in $N_{2}=40$ planning cycles, and hence the risk is
bounded by $\Delta_{o}+N_{2}\Delta_{tube}=0.14$, or
$\Delta_{o}+(1-(1-\Delta_{tube})^{N_{2}})=0.1393$.

Figure 5: Outer approximation of the risk contour of a vehicle. The red curve
is the boundary of the risk contour, which is a 4th order polynomial in 3
variables. The blue ellipse is an outer approximation of the risk contour,
which is a quadratic function in 3 variables.

Figure 6: Ground vehicle changing lanes. The plots are the snapshots of the
state space at time $t=0.2,2.2,5.2,5.5,5.8,6.2,6.8,10$, and $13.3$. All
vehicles are moving to the right. Blue ellipses are risk contours of the
surrounding vehicles. The camera frame is relatively static to the top left
vehicle, and hence the top left vehicle looks static across the images. In
this scene, two vehicles on the top have almost the same velocity. The cyan
curve represents the trajectory the system has taken so far. The blue dot at
the right end of the cyan curve is the current system position. The system
reaches the goal at time 13.3.

Comparison with CC-RRT. We compare our method with Chance-Constrained RRT (CC-
RRT) method proposed in [7]. In their setting, the system dynamics is linear
and has Gaussian noise. The initial position of the system has Gaussian
distribution. The obstacles are convex polyhedra.

We apply CC-RRT to the clustered environment and a trajectory is plotted as a
cyan curve in Figure 3 Right. The system dynamics is simplified as

$\displaystyle x_{t+1}$ $\displaystyle=x_{t}+\Delta T*v_{x,t}+w_{x,t}$
$\displaystyle y_{t+1}$ $\displaystyle=y_{t}+\Delta T*v_{y,t}+w_{y,t}$

where $w_{x,t}$ and $w_{y,t}$ are Gaussian distributions with mean 0 and
standard deviation $0.02$. Similar to our online planning method, we do MPC
style planning and at each planning cycle, we do RRT search and the safe
candidates are ranked by an objective function. The objective function is the
distance of the future state to the path given by the high level planner, plus
a weighted distance of the future state to the goal. The safe candidate with
the lowest score is selected and executed.

Our method has several advantages over CC-RRT: (i) Our method can bound the
risk using tubes over the continuous state space, while CC-RRT can only check
safety at certain discrete points via sampling. (ii) Our method can work with
nonlinear systems, while CC-RRT works with linear systems. (iii) Our method
can deal with more general probabilistic distributions, while CC-RRT only
works with Gaussian distribution. (iv) Our method can deal with more general
obstacles, and does not assume the obstacles are convex, though using simple
convex obstacles would speed up our method.

Figure 7: Histogram of the number of planning cycles to reach the goal across
20 scenes.

###  IV-B Vehicle Changing Lanes

In this example, we consider the ground vehicle whose dynamics is given by
Equation 17. We randomly generate 20 scenes. In each scene, there are three
uncertain vehicles moving towards the right. Two are moving on the top lane,
and one is moving on the bottom lane. The initial positions and velocities of
the vehicles vary across the 20 scenes. Each vehicle has the form

$\displaystyle\\{(x,y)\in\mathbb{R}^{2}:(x-x_{i})^{2}+4(y-y_{i})^{2}\leq
w^{2}\\}$

where $w$ has uniform distribution over $[0.3,0.4]$. The system starts from
the position $(0,0)$ moving with initial velocity $(1,0)$, and the goal is to
switch lanes from $y=0$ to $y=1$.

We are given a total risk level of 0.3. We estimate an upper bound on the
number of planning cycles is 200. We use risk contours with the risk level of
$\Delta_{o}=0.1$ and $\Delta_{tube}=0.001$. The risk contour is a polynomial
of order 4. We outer approximate the risk contour by an ellipse, which is a
quadratic function, reducing the order of the risk contour to 2 (Figure 5). We
build motion primitives using the analytical method. There is no high level
planner in this example. Instead there is an cost function
$(y-1)^{2}+10\theta^{2}+10^{7}*I(|\theta-\pi/6|>0)$ to minimize, where $y$ and
$\theta$ are the states at the last time step of the motion primitive, and $I$
is the indicator function. So a motion primitive that approaches the lane
$y=1$ while keeping $\theta$ close to 0 is preferable. The task is
accomplished once $y$ is close to 1 and $\theta$ is close to 0. The system re-
plans at every time step.

Among the 20 scenes, the max number of planning cycles is 150 (Figure 7),
which is less than our estimate of 200, and hence the total risk of 0.3 is
guaranteed. We plotted an example trajectory at several different time steps
in Figure 6. In this example, $N=133$, and hence the entire trajectory has
bounded risk of $\Delta_{o}+N\Delta_{tube}=0.233$ or
$\Delta_{o}+(1-(1-\Delta_{tube})^{N})=0.2246$.

##  V Conclusion and Future Work

We have presented a real-time tube-based motion planning approach for
stochastic nonlinear systems in uncertain environments via motion primitives.
Our approach works for long-term tasks, which trajectory optimization methods,
such as the one in [1], cannot be directly applied. Our approach guarantees
the probability of system states colliding with any obstacle is bounded above.
Our approach is very practical and can be deployed on various robotics
systems. Future work includes implementation on real robots and autonomous
driving cars.

## References

  * [1] W. Han, A. Jasour, and B. Williams, “Non-gaussian risk bounded trajectory optimization for stochastic nonlinear systems in uncertain environments,” in _2022 International Conference on Robotics and Automation (ICRA)_ , 2022, pp. 11 044–11 050.
  * [2] B. Axelrod, L. P. Kaelbling, and T. Lozano-Pérez, “Provably safe robot navigation with obstacle uncertainty,” _The International Journal of Robotics Research_ , vol. 37, no. 13-14, pp. 1760–1774, 2018.
  * [3] W. Schwarting, J. Alonso-Mora, L. Pauli, S. Karaman, and D. Rus, “Parallel autonomy in automated vehicles: Safe motion generation with minimal intervention,” in _2017 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 1928–1935.
  * [4] S. Dai, S. Schaffert, A. Jasour, A. Hofmann, and B. Williams, “Chance constrained motion planning for high-dimensional robots,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2019.
  * [5] T. Lew, R. Bonalli, and M. Pavone, “Chance-constrained sequential convex programming for robust trajectory optimization,” in _2020 European Control Conference (ECC)_. IEEE, 2020, pp. 1871–1878.
  * [6] C. Dawson, A. Jasour, A. Hofmann, and B. Williams, “Provably safe trajectory optimization in the presence of uncertain convex obstacles,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2020.
  * [7] B. Luders, M. Kothari, and J. How, “Chance constrained rrt for probabilistic robustness to environmental uncertainty,” in _AIAA guidance, navigation, and control conference_ , 2010, p. 8160.
  * [8] L. Blackmore and M. Ono, “Convex chance constrained predictive control without sampling,” in _AIAA Guidance, Navigation, and Control Conference_ , 2009, p. 5876.
  * [9] L. Blackmore, M. Ono, A. Bektassov, and B. C. Williams, “A probabilistic particle-control approximation of chance-constrained stochastic predictive control,” _IEEE transactions on Robotics_ , vol. 26, no. 3, pp. 502–517, 2010.
  * [10] L. Janson, E. Schmerling, and M. Pavone, “Monte carlo motion planning for robot trajectory optimization under uncertainty,” in _Robotics Research_. Springer, 2018, pp. 343–361.
  * [11] G. C. Calafiore and M. C. Campi, “The scenario approach to robust control design,” _IEEE Transactions on automatic control_ , vol. 51, no. 5, pp. 742–753, 2006.
  * [12] M. Cannon, “Chance-constrained optimization with tight confidence bounds,” _arXiv preprint arXiv:1711.03747_ , 2017.
  * [13] K. Okamoto, “Optimal covariance steering: Theory and its application to autonomous driving,” Ph.D. dissertation, Georgia Institute of Technology, 2019\.
  * [14] T. M. Howard, C. J. Green, A. Kelly, and D. Ferguson, “State space sampling of feasible motions for high-performance mobile robot navigation in complex environments,” _Journal of Field Robotics_ , vol. 25, no. 6-7, pp. 325–345, 2008.
  * [15] T. M. Howard and A. Kelly, “Optimal rough terrain trajectory generation for wheeled mobile robots,” _The International Journal of Robotics Research_ , vol. 26, no. 2, pp. 141–166, 2007.
  * [16] X. Yang, A. Agrawal, K. Sreenath, and N. Michael, “Online adaptive teleoperation via motion primitives for mobile robots,” _Autonomous Robots_ , vol. 43, no. 6, pp. 1357–1373, 2019.
  * [17] C. L. Bottasso, D. Leonello, and B. Savini, “Path planning for autonomous vehicles by trajectory smoothing using motion primitives,” _IEEE Transactions on Control Systems Technology_ , vol. 16, no. 6, pp. 1152–1168, 2008\.
  * [18] T. Löw, T. Bandyopadhyay, J. Williams, and P. V. Borges, “Prompt: Probabilistic motion primitives based trajectory planning.” in _Robotics: Science and Systems_ , 2021.
  * [19] A. Paraschos, C. Daniel, J. R. Peters, and G. Neumann, “Probabilistic movement primitives,” _Advances in neural information processing systems_ , vol. 26, 2013.
  * [20] A. Jasour, W. Han, and B. Williams, “Convex risk bounded continuous-time trajectory planning in uncertain nonconvex environments,” in _Robotics: Science and Systems_ , 2021.
  * [21] J. Lofberg, “Yalmip: A toolbox for modeling and optimization in matlab,” in _International conference on robotics and automation_ , 2004.
  * [22] M. M. Tobenkin, F. Permenter, and A. Megretski, “spotless: Polynomial and conic optimization,” 2013. [Online]. Available: github.com/spot-toolbox/spotless
  * [23] A. Jasour, W. Han, and B. Williams, “Real-time risk-bounded tube-based trajectory safety verification,” in _IEEE Conference on Decision and Control (CDC)_ , 2021.
  * [24] A. Jasour and B. Williams, “Risk contours map for risk bounded motion planning under perception uncertainties.” in _Robotics: Science and Systems_ , 2019\.
  * [25] A. Jasour, A. Wang, and B. C. Williams, “Moment-based exact uncertainty propagation through nonlinear stochastic autonomous systems,” _arXiv preprint arXiv:2101.12490_ , 2021.
  * [26] J. Gu, C. Sun, and H. Zhao, “Densetnt: End-to-end trajectory prediction from dense goal sets,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 15 303–15 312.
  * [27] Q. Lu, W. Han, J. Ling, M. Wang, H. Chen, B. Varadarajan, and P. Covington, “Kemp: Keyframe-based hierarchical end-to-end deep model for long-term trajectory prediction,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 646–652.


Department od Mathematics, University of Zagreb, Bijenička cesta 30, Zagreb,
Croatia<EMAIL_ADDRESS>https://web.math.pmf.unizg.hr/ tguzvic/ Department of
Mathematics, University of Zagreb, Bijenička cesta 30, Zagreb, Croatia
<EMAIL_ADDRESS>https://www.math.pmf.unizg.hr/hr/borna-vukorepa Elliptic
curves [2010]11G05

# Torsion groups of elliptic curves over $\mathbb{Q}(\mu_{p^{\infty}})$

Tomislav Gužvić  Borna Vukorepa

###### Abstract

Let $E/\Q$ be an elliptic curve and $p\in\\{5,7,11\\}$ be a prime. We
determine the possibilities for $E(\mathbb{Q}(\zeta_{p}))_{tors}$.
Additionally, we determine all the possibilities for
$E(\mathbb{Q}(\zeta_{16}))_{tors}$ and $E(\mathbb{Q}(\zeta_{27}))_{tors}$.
Using these results we are able to determine the possibilities for
$E(\mathbb{Q}(\mu_{p^{\infty}}))_{tors}$.

##  1 Introduction

Let $K$ be a number field such that $[K:\mathbb{Q}]=d$ and let $E/K$ be an
elliptic curve. A celebrated theorem of Mordell and Weil shows that $E(K)$ is
a finitely generated abelian group. Therefore this group can be decomposed as
$E(K)=E(K)_{tors}\oplus\mathbb{Z}^{r}$, $r\geq 0$. It is known that
$E(K)_{tors}$ is of the form $C_{m}\oplus C_{n}$ for two positive integers
$m,n$ such that $m$ divides $n$, where $C_{m}$ and $C_{n}$ denote cyclic
groups of order $m$ and $n$, respectively.

One of the goals in the theory of elliptic curves is the classification of
torsion groups of elliptic curves defined over various fields.

Let $d$ be a positive integer. Define $\Phi(d)$ to be the set of possible
isomorphism classes of groups $E(K)_{tors}$, where $K$ runs through all number
fields $K$ of degree $d$ and $E$ runs through all elliptic curves over $K$. In
[21], Merel proved that $\Phi(d)$ is finite for all positive integers $d$. The
set $\Phi(1)$ can be seen in Theorem 1 and was determined by Mazur
[mazurtorzija]. [Mazur, [mazurtorzija]] Let $E/\mathbb{Q}$ be an elliptic
curve. Then

$E(\mathbb{Q})_{tors}\cong\cases{C}_{m},&m=1,...,10,12,\\\ C_{2}\oplus
C_{2m},m=1,...,4.$

The set $\Phi(2)$ has been determined by Kenku, Momose and Kamienny
[kenkumomose], [kamienny]. Derickx, Etropolski, Hoeij, Morrow and Zureick-
Brown have determined $\Phi(3)$ in [hoeij].

Define $\Phi^{\mathrm{CM}}(d)$ to be the set of possible isomorphism classes
of groups $E(K)_{tors}$, where $K$ runs through all number fields $K$ of
degree $d$ and $E$ runs through all elliptic curves with complex
multiplication ($\mathrm{CM}$). The set $\Phi^{\mathrm{CM}}(1)$ has been
determined by Olson in [32] and $\Phi^{\mathrm{CM}}(d)$ for $d=2,3$ by Zimmer
and his collaborators in [80], [81] and [82]. The sets
$\Phi^{\mathrm{CM}}(d)$, for $4\leq d\leq 13$ have been determined by Clark,
Corn, Rice and Stankewicz in [30]. Bourdon, Pollack and Stankewicz have
determined torsion groups of $\mathrm{CM}$ elliptic curves over odd degree
number fields in [bourdonpollack].

Define $\Phi_{\mathbb{Q}}(d)\subseteq\Phi(d)$ to be the set of possible
isomorphism classes of groups $E(K)_{tors}$, where $K$ runs through all number
fields $K$ of degree $d$ and $E$ runs through all elliptic curves defined over
$\mathbb{Q}$. For $d=2,3$, the sets $\Phi_{\mathbb{Q}}(d)$ have been
determined by Najman [najman_ratcubic].

Let $E/\mathbb{Q}$ be an elliptic curve and $K/\mathbb{Q}$ a quadratic
extension. Then $E(K)_{tors}$ is isomorphic to the one of the following
groups: C_m, m = 1, 2, …, 9, 10, 12, 15, 16
C_2 ⊕C_2m, m = 1, 2, 3, 4, 5, 6
C_3 ⊕C_3m, m = 1, 2
C_4 ⊕C_4. $C_{15}$ is the only group which appears in only finitely many
cases, and only over the extensions $\mathbb{Q}(\sqrt{5})$ and
$\mathbb{Q}(\sqrt{-15})$.

Let $E/\mathbb{Q}$ be an elliptic curve and $K/\mathbb{Q}$ a cubic extension.
Then $E(K)_{tors}$ is isomorphic to the one of the following groups: C_m, m =
1, 2, …, 10, 12, 13, 14, 18, 21
C_2 ⊕C_2m, m = 1, 2, 3, 4, 7. $C_{21}$ is the only group which appears in only
finitely many cases, and only over the extension $\mathbb{Q}(\zeta_{9})^{+}$.

The set $\Phi_{\mathbb{Q}}(4)$ has been determined by Chou [chou2] and
González-Jiménez and Najman [growth]. The set $\Phi_{\mathbb{Q}}(5)$ has been
determined by González-Jiménez in [38]. González-Jiménez and Najman have also
proved that $\Phi_{\mathbb{Q}}(p)=\Phi(1)$ for primes $p\geq 7$ in [growth].
For an odd prime $\ell$ and a positive integer $d$, Propp [90] has determined
when there exists a degree $d$ number field $K$ and an elliptic curve $E/K$
with $j(E)\in\mathbb{Q}\setminus\\{0,1728\\}$ such that $E(K)_{tors}$ contains
a point of order $\ell$.

Let $\mu_{n}$, for positive integer $n$, be the set of all complex numbers
$\omega$ such that $\omega^{n}=1$. Note that for a prime number $p$ we have
that $\Q(\mu_{p})=\Q(\zeta_{p})$, where $\zeta_{p}$ is, as usual, $p$th
primitive root of unity.

For a prime number $p$, we define a set $\mu_{p^{\infty}}$ as the set of all
complex numbers $\omega$ for which there exists non-negative integer $k$ such
that $\omega^{p^{k}}=1$. Note that $\Qz$ is the set $\Q$ extended with all
$p^{n^{\text}{th}}$primitive roots of unity.

In [guzkri], the authors considered the following problem: assume that
$E/\mathbb{Q}$ is an elliptic curve, $p$ a prime number and
$K=\Q(\mu_{p}^{\infty})$. They show that the torsion subgroup of $E$ grows
only over small subfields of $K$. More precisely, they showed the following:

Let $E/\Q$ be an elliptic curve, then for a prime number $p\geq 5$ it holds
that

$E(\Qz)_{\tors}=E(\Qz[p])_{\tors}.$

Furthermore,

$E(\Qz[3][\infty])_{\tors}=E(\Qz[3^{3}])_{\tors}\qquad\text{and}\qquad
E(\Qz[2][\infty])_{\tors}=E(\Qz[2^{4}])_{\tors}.$

This result is “the best possible”. For
$E=\href{http://www.lmfdb.org/EllipticCurve/Q/27a4}{27a4}$ we have that

$E(\Qz[3^{2}])_{\tors}=C_{9}\subsetneq C_{27}=E(\Qz[3^{3}])_{\tors}$

and for $E=\href{http://www.lmfdb.org/EllipticCurve/Q/32a4}{32a4}$ it holds
that

$E(\Qz[2^{3}])_{\tors}=C_{2}\oplus C_{4}\subsetneq C_{2}\oplus
C_{8}=E(\Qz[2^{4}])_{\tors}.$

It becomes natural to ask how can the torsion group of $E/\mathbb{Q}$ grow
when we consider the base change $E/\mathbb{Q}(\zeta_{p})$. This becomes much
harder then it seems as $p$ grows because our methods sometimes rely on pure
computation.
Our results are the following theorems.

Let $E/\Q$ be an elliptic curve. Then $E(\Q(\zeta_{16}))_{\tors}$ is one of
the following groups (apart from those in Mazur’s theorem):

$C_{4}\oplus
C_{4}\quad\href{https://www.lmfdb.org/EllipticCurve/Q/15a1}{(15a1)},\qquad
C_{2}\oplus
C_{10}\quad\href{https://www.lmfdb.org/EllipticCurve/Q/2112bd2}{(2112bd2)}.$

Let $E/\Q$ be an elliptic curve. Then $E(\Q(\zeta_{27}))_{\tors}$ can be only
one of the next two groups (apart from those in Mazur’s theorem):

$C_{3}\oplus C_{3},\qquad C_{3}\oplus C_{6},\qquad C_{3}\oplus C_{9},\qquad
C_{21},\qquad C_{27}.$

Additionally, in Lemma 5 we give a description of the set of possible group
structures $E(\mathbb{Q}(\zeta_{p}))_{tors}$ can be isomorphic to for $p=5,7$
and $11$, where $E/\mathbb{Q}$ is an elliptic curve.
We discuss further attempts to classify $E(\mathbb{Q}(\zeta_{p}))_{tors}$, for
arbitrary prime number $p$.
Magma [3] code used in this paper can be found https://github.com/brutalni-
vux/TorsionCyclotomic?fbclid=IwAR32kNBmGbzdULcwWpmXhP5QZ0Ap4QjZf4vqSdbFJ_9dWko6d16W3YtRQpQhere.

##  2 Notation and auxiliary results

Let $E/F$ be an elliptic curve defined over a number field $F$. There exists
an $F$-rational cyclic isogeny $\phi\colon E\to E^{\prime}$ of degree $n$ if
and only if $\langle P\rangle$, where $P\in E(\overline{F})$ is a point of
order $n$, is a $\Gal(\overline{F}/F)$-invariant group; in this case we say
that $E$ has an $F$-rational $n$-isogeny. When $F=\mathbb{Q}$, the possible
degrees of $n$-isogenies of elliptic curves over $\mathbb{Q}$ are known by the
following theorem.

[Mazur [mazurizog], Kenku [kenku1], [kenku2], [kenku3], [kenku4]] Let
$E/\mathbb{Q}$ be an elliptic curve with a rational $n$-isogeny. Then

$n\in\\{1,...,19,21,25,27,37,43,67,163\\}.$

There are infinitely many elliptic curves (up to
$\overline{\mathbb{Q}}$-isomorphism) with a rational $n$-isogeny over
$\mathbb{Q}$ for

$n\in\\{{1,...,10,12,13,16,18,25}\\}$

and only finitely many for all the other $n$. If $E$ does not have complex
multiplication, then $n\leq 18$ with $n\neq 14$ or $n\in\\{21,25,37\\}.$

### Galois representations

Let $E/\mathbb{Q}$ be an elliptic curve and let $n$ a positive integer. The
field $\mathbb{Q}(E[n])$ is the number field obtained by adjoining to
$\mathbb{Q}$ all the $x$ and $y$-coordinates of the points of $E[n]$. The
absolute Galois group $\Gal(\overline{\mathbb{Q}}/\mathbb{Q})$ acts on $E[n]$
by its action on the coordinates of the points, inducing a mod $n$ Galois
representation attached to $E$:

$\rho_{E,n}:\Gal(\overline{\mathbb{Q}}/\mathbb{Q})\to\Aut(E[n]).$

After we fix a basis for the $n$-torsion, we can identify $\Aut(E[n])$ with
$\GL_{2}(\mathbb{Z}/n\mathbb{Z})$. This means that we can consider
$\rho_{E,n}(\Gal(\overline{\mathbb{Q}}/\mathbb{Q}))$ as a subgroup of
$\GL_{2}(\mathbb{Z}/n\mathbb{Z})$, uniquely determined up to conjugacy. We
shall denote $\rho_{E,n}(\Gal(\overline{\mathbb{Q}}/\mathbb{Q}))$ by
$G_{E}(n)$. Moreover, since $\mathbb{Q}(E[n])$ is a Galois extension of
$\mathbb{Q}$ and
$\ker\rho_{E,n}=\Gal(\overline{\mathbb{Q}}/\mathbb{Q}(E[n]))$, by the first
isomorphism theorem we have $G_{E}(n)\cong\Gal(\mathbb{Q}(E[n])/\mathbb{Q})$.

We would like to know what are the possibilities for $G_{E}(n)$ as a subgroup
of $\GL(C_{n})$. For some values of $n$, this can be seen in Tables 1 and 2.
For most values of $n$ we do not have a list of possibilities of $G_{E}(n)$,
but we have a result that helps us see if for a given matrix subgroup $M$ of
$\GL(C_{n})$ there exists an elliptic curve $E/\mathbb{Q}$ such that
$\rho_{E,n}(\Gal(\overline{\mathbb{Q}}/\mathbb{Q}))=M$ (up to conjugation).
The following lemma is well-known and it will be useful in cases where we
analyze torsion group of elliptic curve $E/\mathbb{Q}$ over a maximal real
subfield of $\mathbb{Q}(\zeta_{p})$. Let $E/\Q$ be an elliptic curve and $L/K$
a quadratic extension of number fields with $L=K(\sqrt{d})$. Let
$E(K)_{(2^{\prime})}$ be the group of $K$-rational points of $E$ of odd order.
Then we have:

$E(K(\sqrt{d}))_{(2^{\prime})}\cong E(K)_{(2^{\prime})}\oplus
E^{d}(K)_{(2^{\prime})}.$

Since all cyclotomic extensions are Galois over $\mathbb{Q}$, the following
result imposes restrictions on the possibilities for torsion subgroup of
$E/\mathbb{Q}$ over a cyclotomic fields. Let $E/\mathbb{Q}$ be an elliptic
curve, $m,n\in\mathbb{N}$ and $K$ a finite Galois extension of $\mathbb{Q}$.
Let $E(K)[mn]\cong C_{m}\oplus C_{mn}$ and $P\in E(K)$ point of order $mn$.
Then we have:

$[\mathbb{Q}(mP):\mathbb{Q}]\mid M(\phi(n),[K:\mathbb{Q}]),$

where $M(\cdot,\cdot)$ is the greatest common divisor and $\phi$ is the Euler
function. Let $P$ be a point of order $mn$ with coordinates in $K$. Then we
can take $Q\in E[mn]$ such that $\\{P,Q\\}$ is a basis for $E[mn]$. Consider
the Galois representation modulo $n$ with respect to $E$:

$\rho:\Gal(\overline{\mathbb{Q}}/\mathbb{Q})\to\GL_{2}(\mathbb{Z}/n\mathbb{Z}).$

Let $\sigma\in\Gal(K/\mathbb{Q})$. Then we have $P^{\sigma}=\alpha P+\beta Q$
for some $\alpha,\beta\in C_{mn}$ because the action of $\sigma$ on $P$
preserves the order of a point.

Now we have $P^{\sigma}-\alpha P=\beta Q$, so $\beta Q\in E(K)$. From that
follows $m\beta\equiv 0\pmod{mn}$. Multiplying by $m$ gives us
$(mP)^{\sigma}=\alpha(mP)$ so $(mP)^{\sigma}\in\langle mP\rangle$ for all
$\sigma\in\Gal(K/\mathbb{Q})$. Because of preserving the order, $\alpha$ has
to be in $(\mathbb{Z}/n\mathbb{Z})^{\times}$.

Since by considering the restriction map we get
$\Gal(K/\mathbb{Q})\cong\Gal(\overline{\mathbb{Q}}/\mathbb{Q})/\Gal(\overline{\mathbb{Q}}/K)$,
we have $(mP)^{\sigma}\in\langle mP\rangle$ for all
$\sigma\in\Gal(\overline{\mathbb{Q}}/\mathbb{Q})$. Therefore, for all
$\sigma\in\Gal(\overline{\mathbb{Q}}/\mathbb{Q})$:

$\rho(\sigma)=\pmatrix{\varphi}(\sigma)&\tau(\sigma)\\\ 0\psi(\sigma),$

where $\varphi,\psi,\tau:\Gal(\overline{\mathbb{Q}}/\mathbb{Q})\rightarrow
C_{n}$, and $\varphi,\psi$ are homomorphisms with image in $(C_{n})^{\times}$.
We know that $(mP)^{\sigma}=g(mP)\Leftrightarrow\varphi(\sigma)=g$, for all
$\sigma\in\Gal(\overline{\mathbb{Q}}/\mathbb{Q})$. Therefore, we have:

$|Im(\varphi)|=|\\{(mP)^{\sigma}:\sigma\in\Gal(K/\mathbb{Q})\\}|=|Orb(mP)|.$

It is clear that Stab$(mP)=\Gal(K/\mathbb{Q}(mP))$, so by orbit and stabilizer
theorem we have:

$|Im(\varphi)|=\frac{|\Gal(K/\mathbb{Q})|}{|\Gal(K/\mathbb{Q}(mP))|}=[\mathbb{Q}(mP):\mathbb{Q}].$

On the other hand, we have $Im(\varphi)\leq(\mathbb{Z}/n\mathbb{Z})^{\times}$,
so we have:

$[\mathbb{Q}(mP):\mathbb{Q}]\mid\phi(n).$

$[\mathbb{Q}(mP):\mathbb{Q}]\mid[K:\mathbb{Q}]$ is obvious and the proof is
complete. One of the crucial results that we will need is the main result from
[chou]: Let $E/\mathbb{Q}$ be a rational elliptic curve. Then
$E(\mathbb{Q}^{ab})_{tors}$ is isomorphic to the one of the following groups:
C_m, m = 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 25, 27, 37, 43, 67, 163
C_2 ⊕C_2m, m = 1, 2, …, 8, 9
C_3 ⊕C_3m, m = 1, 3
C_4 ⊕C_4m, m = 1, 2, 3, 4
C_5 ⊕C_5,
C_6 ⊕C_6,
C_8 ⊕C_8.

This means that all of our candidate torsion subgroups are the subgroups of
the groups in the above list. Our approach will mainly consist of eliminating
a certain set of possibilities from the list above in order to classify
torsion groups of elliptic curves over a specific cyclotomic field.

##  3 Torsion growth over $\mathbb{Q}(\zeta_{16})$Q16

Assume that $E/\mathbb{Q}$ is an elliptic curve and that $C_{m}\oplus
C_{mn}\subseteq E(\mathbb{Q}(\zeta_{16}))_{tors}$. By the properties of the
Weil pairing, we have $\mathbb{Q}(\zeta_{m})\subseteq\mathbb{Q}(\zeta_{16})$.
It follows that $m\in\\{1,2,4,8\\}.$ We first eliminate a certain amount of
cyclic groups listed in Theorem 2.

Let $E/\Q$ be an elliptic curve. Then $E(\Q(\zeta_{16}))_{\tors}$ is not
isomorphic to $C_{n}$ if

$n\in\\{11,14,17,18,19,21,25,27,37,43,67,163\\}.$

Lemma 2 gives us that if $P_{n}$ is a point of order
$n\not\in\\{17,21,25,37\\}$, we have $[\Q(P_{n}):\Q]\mid 2$, which is
impossible by Theorem 1. By the same lemma we get that if $P_{n}$ is a point
of order $n\in\\{21,25,37\\}$, then we have $[\Q(P_{n}):\Q]\mid 4$, which is
impossible by [chou2, Theorem 1.4].
It remains to consider the case $n=17$. By [growth, Theorem 5.8] we conclude
that the point $P_{17}$ of order $17$ cannot be defined over some strictly
smaller subfield of $\Q(\zeta_{16})$. That means that all
$\sigma\in\Gal(\Q(\zeta_{16})/\Q)$ act differently on $P_{17}$. Since
$\Gal(\Q(\zeta_{16}))$ has four elements $\sigma$ such that $\sigma^{2}=id$,
we have that $P_{17}^{\sigma^{2}}=k^{2}P_{17}=P_{17}$ for four different
$\sigma$. That means that we have $k^{2}\equiv 1\pmod{17}$ for four different
$k$, a contradiction.

After eliminating plenty of cyclic groups, we discuss the cases when $E$
obtains full $2$-torsion over $\mathbb{Q}(\zeta_{16}).$ This is done by the
following lemmas:

Let $E/\Q$ be an elliptic curve. Then $E(\Q(\zeta_{16}))_{\tors}$ is not
isomorphic to $C_{2}\oplus C_{14}$ or $C_{2}\oplus C_{18}$. We will prove the
result for $C_{2}\oplus C_{14}$ and the proof for the case $C_{2}\oplus
C_{18}$ is identical.
Let $P_{14}\in E(\Q(\zeta_{16}))$ be the point of order $14$. From Lemma 2 we
get that $[\mathbb{Q}(2P_{14}):\mathbb{Q}]\mid 2$. It is also well-known that
$[\mathbb{Q}(E[2]):\mathbb{Q}]\in\\{1,2,3,6\\}$. Since $E[2]$ is defined over
$\Q(\zeta_{16})$, we have $[\mathbb{Q}(E[2]):\mathbb{Q}]\in\\{1,2\\}$.
Let $Q_{2}$ be a point of order $2$ different from $7P_{14}$. We now have
$[\mathbb{Q}(2P_{14},7P_{14},Q_{2}):\mathbb{Q}]\mid 4$. Since
$2P_{14},7P_{14}$ and $Q_{2}$ generate our torsion subgroup $C_{2}\oplus
C_{14}$, we now know that this torsion subgroup appears over some strictly
smaller subfield of $\Q(\zeta_{16})$.
Now we get a contradiction by using Theorem 1 and [chou2, Theorem 1.4].

Let $E/\Q$ be an elliptic curve. Then $E(\Q(\zeta_{16}))_{\tors}$ is not
isomorphic to $C_{15}$, $C_{2}\oplus C_{12}$, $C_{4}\oplus C_{12}$,
$C_{4}\oplus C_{8}$, $C_{4}\oplus C_{16}$ or $C_{8}\oplus C_{8}$. Both
$X_{1}(15)$ and $X_{1}(2,12)$ are elliptic curves. A computation in Magma
shows that $X_{1}(15)(\Q(\zeta_{16}))$ and $X_{1}(15)(\Q)$ have the same
Mordell-Weil group structure. Therefore, since $C_{15}\not\in\Phi(1)$, it also
cannot appear over $\Q(\zeta_{16})$.
The curves $X_{1}(2,12)(\Q(\zeta_{16}))$ and $X_{1}(2,12)(\Q(i))$ have the
same Mordell-Weil group structure. It was proven in [najCyclo, Lemma 7] that
$C_{2}\oplus C_{12}$ does not appear as a torsion subgroup over $\Q(i)$.
Therefore, it also cannot appear over $\Q(\zeta_{16})$. This also covers the
case $C_{4}\oplus C_{12}$.
We consider the modular curves $X_{1}(4,8)(\Q(\zeta_{16}))$ and
$X_{1}(4,8)(\Q(\zeta_{8}))$ which are actually elliptic curves. A computation
in Magma shows that $X_{1}(4,8)(\Q(\zeta_{16}))$ has rank $0$ and the same
torsion as $X_{1}(4,8)(\Q(\zeta_{8}))$, which contains only cusps, see
[najman_bruin, Case 6.11]. This also covers the cases $C_{8}\oplus C_{8}$ and
$C_{4}\oplus C_{16}$.

The following lemma is a bit more complicated than the previous ones. The idea
is to consider the corresponding modular curve $X_{1}(16)$ and its Jacobian
$J_{1}(16)$ over some cyclotomic fields in order to determine that the
Jacobian has rank $0$. After that, we determine torsion of
$J_{1}(\mathbb{Q}(\zeta_{16}))$ and consequently the number of points on
$X_{1}(16)(\mathbb{Q}(\zeta_{16}))$, all of which turn out to be cusps.

Let $E/\Q$ be an elliptic curve. Then $E(\Q(\zeta_{16}))_{\tors}$ is not
isomorphic to $C_{16}$, $C_{2}\oplus C_{16}$ or $C_{13}$. We consider the
modular curve $X_{1}(16)(\Q(\zeta_{16}))$ and its Jacobian
$J_{1}(16)(\Q(\zeta_{16}))$. We will demonstrate the use of standard methods
for determining points on $X_{1}(16)(\Q(\zeta_{16}))$.
A computation in Magma shows that $r(J_{1}(16)(\Q(\zeta_{16})))=0$. Since

$r(J_{1}(16)(\Q(\zeta_{8})(\sqrt{\zeta_{8}})))=r(J_{1}(16)(\Q(\zeta_{8})))+r(J_{1}^{\zeta_{8}}(16)(\Q(\zeta_{8}))),$

the computation of the rank becomes shorter and we obtain that the rank of our
Jacobian is $0$.
Now we determine $J_{1}(16)(\Q(\zeta_{16}))_{tors}$. Rational prime $p=17$
splits completely in $\Q(\zeta_{16})$ so by reducing modulo some $\mathfrak p$
that lies above $p$ we get an injection

$red_{\mathfrak p}:J_{1}(16)(\Q(\zeta_{16}))_{tors}\to J_{1}(16)(\F_{17}).$

This map is injective due to the result of Katz [katz].
A computation in Magma shows that $|J_{1}(16)(\F_{17})|=400.$ It follows that
$|J_{1}(16)(\Q(\zeta_{16}))|\leq 400$. By using the generators of the
$2$-torsion subgroup of $J_{1}(16)(\Q(\zeta_{16}))$ and some elements of
$J_{1}(16)(\Q(\zeta_{16}))$ that we get from some known points on
$X_{1}(16)(\Q(\zeta_{16}))$, we are able to generate a group with $400$
elements. Therefore, we know exactly how $J_{1}(16)(\Q(\zeta_{16}))$ looks
like.
Now we are able to determine all points on $X_{1}(16)(\Q(\zeta_{16}))$ by
considering the Mumford representations of the elements of
$J_{1}(16)(\Q(\zeta_{16}))$. We easily get that
$|X_{1}(16)(\Q(\zeta_{16}))|=14$ with all points being cusps. Therefore, we
can conclude that there are no elliptic curves $E/\Q(\zeta_{16})$ (and
consequently $E/\Q$) with a point of order $16$ over $\Q(\zeta_{16})$.
It remains to show that $E(\Q(\zeta_{16}))_{\tors}$ can’t be $C_{13}$. We
consider the modular curve $X_{1}(13)(\Q(\zeta_{16}))$ and its Jacobian
$J_{1}(13)(\Q(\zeta_{16}))$.

A computation in Magma shows that $r(J_{1}(13)(\Q(\zeta_{16})))=0$. As in the
previous lemma, we obtain:

$r(J_{1}(13)(\Q(\zeta_{16})))=r(J_{1}(16)(\Q(\zeta_{8})))+r(J_{1}^{\zeta_{8}}(16)(\Q(\zeta_{8})))=0.$

The next step is to determine $J_{1}(13)(\Q(\zeta_{16}))_{tors}$. We determine
the two-torsion subgroup, which turns out to be trivial. Using the result of
Katz [katz], we get an injection:

$red_{\mathfrak p}:J_{1}(13)(\Q(\zeta_{16}))_{tors}\to J_{1}(13)(\F_{17}).$

We also get that rational prime $q=41$ has inertia degree $2$ in
$\Q(\zeta_{16})$ so we have another injection:

$red_{\mathfrak q}:J_{1}(13)(\Q(\zeta_{16}))_{tors}\to
J_{1}(13)(\F_{41^{2}}).$

We notice that $gcd(\\#J_{1}(13)(\F_{17}),\\#J_{1}(13)(\F_{41^{2}}))=76$, so
$\\#J_{1}(13)(\Q(\zeta_{16}))\mid 76$.
Since the two torsion subgroup is trivial, we get that
$\\#J_{1}(13)(\Q(\zeta_{16}))\mid 19$. We can find a point of order $19$ on
our Jacobian. By checking the Mumford representations of those divisors, we
find that all of the points on the Jacobian come from cusps on
$X_{1}(13)(\Q(\zeta_{16}))$ (and actually $X_{1}(13)(\Q))$.
Therefore, we can conclude that there are no elliptic curves
$E/\Q(\zeta_{16})$ (and consequently $E/\Q$) such that
$E(\Q(\zeta_{16}))_{\tors}\cong C_{13}$.

##  4 Torsion growth over $\mathbb{Q}(\zeta_{27})$Q27

In this section we prove Theorem 1 using a series of lemmas. First we
eliminate some possibilities for a cyclic group to appear as the subgroup of
$E(\mathbb{Q}(\zeta_{27}))$. Assume that $E/\mathbb{Q}$ is an elliptic curve
and that $C_{m}\oplus C_{mn}\subseteq E(\mathbb{Q}(\zeta_{27}))_{tors}$. By
the properties of the Weil pairing and taking the Theorem 2 into account, we
have $\mathbb{Q}(\zeta_{m})\subseteq\mathbb{Q}(\zeta_{27})$. It follows that
$m\in\\{1,2,3,6\\}.$

Let $E/\Q$ be an elliptic curve. Then $E(\Q(\zeta_{27}))_{tors}$ is not
isomorphic to $C_{n}$ if $n\in\\{11,13,14,15,16,17,19,25,37,43,67,163\\}$.
$C_{11},C_{25}$: If $n\in\\{11,25\\}$ is a point of order $n$, then then by
Lemma 2 we have $[\Q(P_{n}):\Q]|2$, which is impossible by Theorem 1.
$C_{13}$: Let $P_{13}\in E(\Q(\zeta_{27}))$ be the point of order $13$. By
Lemma 2 we have $[\Q(P_{13}):\Q]\mid 6$. Therefore, this torsion subgroup is
defined over $\Q(\zeta_{9})$. Theorem 1 tells us that this torsion subgroup
cannot be defined over quadratic field. Therefore, it is defined over sextic
or cubic field. Assume it is defined over sextic field (entire
$\Q(\zeta_{9})$). Then we can use Lemma 2 to get:

$C_{13}\cong E(\Q(\zeta_{9}))_{(2^{\prime})}\cong
E(\Q(\zeta_{9}^{+}))_{(2^{\prime})}\oplus
E^{-3}(\Q(\zeta_{9}^{+}))_{(2^{\prime})}.$

This means that either $E$ or $E^{-3}$ has torsion subgroup $\Z/13\Z$ defined
over $\Q(\zeta_{9})^{+}$.

Now we will be finished if we prove that torsion subgroup $\Z/13\Z$ cannot
appear over $\Q(\zeta_{9})^{+}$. To do this, we consider
$X_{1}(13)(\Q(\zeta_{9})^{+})$.
As before, we use Magma to determine that $J_{1}(13)(\Q(\zeta_{9})^{+})\cong
J_{1}(13)(\Q)$ and that all points on the Jacobian come from cusps, which
completes the proof.
$C_{14}$: Let $P_{14}\in E(\mathbb{Q}(\zeta_{27}))$ be a point of order $14$.
By Lemma 2 it follows that $[\mathbb{Q}(2P_{14}):\mathbb{Q}]$ divides $6$, so
$\mathbb{Q}(2P_{14})$ is contained in $\mathbb{Q}(\zeta_{9})$. The point
$7P_{14}$ of order $2$ satisfies
$[\mathbb{Q}(7P_{14}):\mathbb{Q}]\in\\{1,2,3\\}$, which means that it is also
contained in $\mathbb{Q}(\zeta_{9})$. It follows that $P_{14}\in
E(\mathbb{Q}(\zeta_{9})).$ Consider the modular curve $X_{1}(14)$. It is an
elliptic curve with LMFDB label
https://www.lmfdb.org/EllipticCurve/Q/14/a/514.a5. On the LMFDB page of the
mentioned curve we can see that its torsion subgroup does not grow in any
number field contained in $\mathbb{Q}(\zeta_{9})$. It remains to show that
$r(E(\mathbb{Q}))=r(E(\mathbb{Q}(\zeta_{9}))=0$, which turns out to be true by
a computation in Magma [3]. Therefore
$X_{1}(14)(\mathbb{Q})=X_{1}(14)(\mathbb{Q}(\zeta_{9}))$ and there does not
exist an elliptic curve over $\mathbb{Q}$ with a point of order $14$ over
$\mathbb{Q}(\zeta_{9})$ and consequently over $\mathbb{Q}(\zeta_{27})$.
$C_{15}$: Let $P_{15}\in E(\mathbb{Q}(\zeta_{27}))$ be a point of order $15$.
Then $3P_{15}$ is a point of order $5$ and $[\mathbb{Q}(3P_{15}):\mathbb{Q}]$
is a divisor of $[\mathbb{Q}(\zeta_{27}):\mathbb{Q}]=18$. By Table 1 we see
that $[\mathbb{Q}(3P_{15}):\mathbb{Q}]\in\\{1,2\\}$. The same way as in the
Lemma 4, $C_{2}\oplus C_{12}$ case, we see that the point $5P_{15}$ of order
$3$ is also defined over at most quadratic extension contained in
$\mathbb{Q}(\zeta_{27})$. Since there is only one quadratic extension
contained in $\mathbb{Q}(\zeta_{27})$, namely $\mathbb{Q}(\zeta_{3})$, we have
$P_{15}\in E(\mathbb{Q}(\zeta_{3}))$, which contradicts the Theorem 1.
$C_{16}$: Let $P_{16}\in E(\mathbb{Q}(\zeta_{27}))$ be a point of order $16$.
Then the point $8P_{16}$ has order $2$ and we have
$[\mathbb{Q}(8P_{16}):\mathbb{Q}]\in\\{1,2,3\\}$. By [growth, Proposition
4.8.] we have $[\mathbb{Q}(P_{16}):\mathbb{Q}]=2^{a}3^{b}$, where $a\geq 0$ is
an integer and $b\in\\{0,1\\}.$ Since the field $\mathbb{Q}(P_{16})$ is
contained in $\mathbb{Q}(\zeta_{27})$, it follows that $2^{a}3^{b}$ divides
$[\mathbb{Q}(\zeta_{27}):\mathbb{Q}]=18$. We conclude that
$[\mathbb{Q}(P_{16}):\mathbb{Q}]\in\\{1,2,3,6\\}.$ Assume that
$[\mathbb{Q}(8P_{16}):\mathbb{Q}]=3$. This means that $\mathbb{Q}(8P_{16})$ is
cyclic, so the entire $2$-torsion is contained in this field and $P_{16}$ is
defined over a number field of degree $3$ or $6$, but such a field is
contained in $\mathbb{Q}(\zeta_{9})$. Therefore we have $C_{2}\oplus
C_{16}\subseteq E(\mathbb{Q}(\zeta_{9})),$ but this is impossible by [tomi1,
Theorem 1.1]. It remains to consider the case when
$[\mathbb{Q}(8P_{16}):\mathbb{Q}]\in\\{1,2\\}$. It follows that
$[\mathbb{Q}(P_{16}):\mathbb{Q}]\in\\{1,2\\}$ again by [growth, Proposition
4.8.]. Since $\mathbb{Q}(P_{16})$ is at most quadratic extension of
$\mathbb{Q}$ contained in $\mathbb{Q}(\zeta_{27})$, it follows that
$\mathbb{Q}(P_{16})\subseteq\mathbb{Q}(\zeta_{3})$. By [najman_torshypN,
Theorem 1] this turns out to be impossible.
$C_{17},C_{37}$: Assume that $n\in\\{17,37\\}$ and that $P_{n}\in
E(\mathbb{Q}(\zeta_{27}))$ is a point of order $n$. By [growth, Theorem 5.8]
it follows that $[\mathbb{Q}(P_{n}):\mathbb{Q}]$ is divisible by $4$, but
since $\mathbb{Q}(P_{n})\subseteq\mathbb{Q}(\zeta_{27})$, this is impossible.
$C_{19}$: Let us consider the case when $n=19$. If $P_{19}\in
E(\mathbb{Q}(\zeta_{27}))$ is a point of order $19$, then $E$ has a rational
$19$-isogeny. By [loz], we have $j(E)=-2^{15}\cdot 3^{3}$. The $19$th division
polynomial $f_{E,19}$ must have a root over $\mathbb{Q}(\zeta_{27})$. Using
Magma, we check that this is not the case and therefore we arrive at the
contradiction.
$C_{43},C_{67},C_{163}$: Assume that $n\in\\{43,67,163\\}$ and $P_{n}\in
E(\mathbb{Q}(\zeta_{27}))[n]$. By [loz, Theorem 2.1] it follows that
$[\mathbb{Q}(P_{n}):\mathbb{Q}]\geq\frac{n-1}{2}>[\mathbb{Q}(\zeta_{27}):\mathbb{Q}]=18,$
a contradiction.

Let $E/\Q$ be an elliptic curve. Then $E(\Q(\zeta_{27}))_{tors}$ is not
isomorphic to $C_{18}$ or $C_{2}\oplus C_{18}$.

If $E(\Q(\zeta_{27}))_{tors}\cong C_{18}$, then Lemma 2 directly gives us that
this torsion subgroup is defined over a number field of degree $6$ which can
only be $\Q(\zeta_{9})$.
If $E(\Q(\zeta_{27}))_{tors}\cong C_{2}\oplus C_{18}$, then Lemma 2 gives us
that if $P_{18}\in E(\Q(\zeta_{27}))_{tors}$ is a point of order $18$, then
$2P_{18}$ is defined over $\Q(\zeta_{9})$. We know that
$[\Q(E[2]):\Q]\in\\{1,2,3,6\\}$, but all unique subextensions of
$\Q(\zeta_{27})$ of those degrees are contained in $\Q(\zeta_{9})$, so again
our torsion subgroup is defined over $\Q(\zeta_{9})$. Now from Lemma 2 we get
that

$C_{9}\cong E(\Q(\zeta_{9}))_{(2^{\prime})}\cong
E(\Q(\zeta_{9}^{+}))_{(2^{\prime})}\oplus
E^{-3}(\Q(\zeta_{9}^{+}))_{(2^{\prime})}.$

Therefore, one of $E(\Q(\zeta_{9}^{+}))$ and $E^{-3}(\Q(\zeta_{9}^{+}))$ has a
point of order $9$. Let $P_{2}\in E(\Q(\zeta_{9}))_{tors}$ be a point of order
$2$. If $[\Q(P_{2}):\Q]\in\\{1,3\\}$, then $P_{2}$ is on both
$E(\Q(\zeta_{9}^{+}))$ and $E^{-3}(\Q(\zeta_{9}^{+}))$. If $[\Q(P_{2}):\Q]=2$,
then there is another point $Q_{2}$ of order $2$ on $E$ defined over $\Q$. In
any case, both $E(\Q(\zeta_{9}^{+}))$ and $E^{-3}(\Q(\zeta_{9}^{+}))$ have a
point of order $2$. Finally, one of them has a point of order $18$. However,
it was proved in [krijanphd, Lemma 3.4.7] that all the points on
$X_{1}(18)(\mathbb{Q}(\zeta_{9})^{+})$ are cusps, which completes the proof.

Let $E/\Q$ be an elliptic curve. Assume that $C_{2}\oplus C_{2n}\cong
E(\Q(\zeta_{27}))_{tors}$. Then $n\in\\{1,2,3,4\\}$. Additionally,
$C_{6}\oplus C_{6}\not\subseteq E(\Q(\zeta_{27}))$.

From Theorem 2 it follows that $n\leq 9$. We have shown that
$E(\Q(\zeta_{27}))$ cannot contain a point of order $18$ in Lemma 4.

$C_{2}\oplus C_{10}$: Let $P_{5}\in E(\Q(\zeta_{27}))$ be the point of order
$5$. It follows that $E$ has a rational $5$-isogeny and
$[\Q(P_{5}):\Q]\in\\{1,2\\}$. If $G_{E}(2)\subseteq 2B$, then by Table 1 we
see that $[\Q(E[2]):\Q]\in\\{1,2\\}$. Thus we have found two at most quadratic
fields contained in $\Q(\zeta_{27})$. Since $\Q(\zeta_{27})$ has an unique
quadratic subextension $F=\Q(\sqrt{-3})$, it follows that
$F=\Q(E[2])=\Q(P_{5})$. We conclude that $C_{2}\oplus C_{10}\subseteq
E(F)_{tors}$, which is impossible by [najman_torshypN, Theorem 1].
Assume that $G_{E}(2)=2Cn$. By Theorem [zywina, Theorem 1.1] it follows that
$j(E)=t^{2}+1728$, for some $t\in\Q$. Since $E$ has a rational $5$-isogeny, by
[zywina, Theorem 1.3] we have $j(E)=\frac{25(s^{2}+10s+5)^{3}}{s^{5}}$, for
some $s\in\Q\setminus\\{0\\}$. It remains to find rational points on the
induced curve. By [tomi3, Page 61] we see that such rational points do not
exist.

$C_{2}\oplus C_{12}$: Let $P_{3}\in E(\Q(\zeta_{27}))$ be the point of order
$3$. The extension $\Q(P_{3})$ is cyclic over $\Q$ since it is a subfield of
$\Q(\zeta_{27})$. By Table 1 we see that $G_{E}(3)$ must be contained in the
Borel subgroup of $\GL_{2}(\mathbb{Z}/3\mathbb{Z})$. Applying Theorem [loz,
Theorem 9.3] it follows that $[\Q(P_{3}):\Q]\in\\{1,2\\}$, so we have
$\Q(P_{3})\subseteq\Q(\zeta_{3})$. Assume that $G_{E}(2)\subseteq 2B$. By
[growth, Proposition 4.6] it follows that $C_{2}\oplus C_{4}\subseteq
E(\Q(\zeta_{3}))$. We conclude that $C_{2}\oplus C_{12}\subseteq
E(\Q(\zeta_{3}))$, which is impossible by [najman_torshypN, Theorem 1].
Consider the case when $G_{E}(2)=2Cn$. By [growth, Proposition 4.8] it follows
that the point $P_{4}$ of order $4$ is defined over cubic or sextic subfield
contained in $\mathbb{Q}(\zeta_{27})$. A computation in Magma [3] shows that
if a point of order $4$ is defined over cubic or sextic number field, then
$G_{E}(2)=\GL_{2}(\mathbb{Z}/2\mathbb{Z})$, a contradiction.

$C_{2}\oplus C_{16}$: By [pprimarnatorzija, Corollary 3.5] it follows that
this is impossible.,

$C_{6}\oplus C_{6}$: Since $\Q(E[6])$ is contained in $\Q(\zeta_{27})$, we
have that $|G_{E}(6)|$ divides $[\Q(\zeta_{27}):\Q]=18$. Additionally, the
group $G_{E}(6)$ is cyclic. If $|G_{E}(6)|<6$, then it follows that $E$
obtains a full $6$-torsion over a number field of degree $1,2$ or $3$, but
this is impossible by Theorem 1 and Theorem 1. Assume that
$|G_{E}(6)|\in\\{6,9\\}$. A search in Magma [3] shows that there exists only
one cyclic group $G$ with such property, namely:

$G:=\Big{\langle}\bmatrix 2&3\\\ 31\Big{\rangle}.$

Reducing $G$ modulo $2$ and modulo $3$ we see that $G_{E}(2)=2Cn$ and
$G_{E}(3)=3Cs.1.1$. By [zywina, Theorem 1.1] we have that $j(E)=t^{2}+1728$,
for some $t\in\Q$. Similary from Theorem [zywina, Theorem 1.2] we have that
$j(E)=f(s)^{3}$, for some rational function $f(s)$ and $s\in\Q$. A computation
in Magma [3] shows that the only affine point on the elliptic curve

$t^{2}+1728=x^{3}$

is $(t,x)=(0,12)$. A direct computation shows that the equation $12=f(s)$ does
not have a rational solution in $\Q$. Therefore $E$ cannot have $C_{6}\oplus
C_{6}$ torsion over $\Q(\zeta_{27})$.

##  5 Torsion growth over
$\mathbb{Q}(\zeta_{5}),\mathbb{Q}(\zeta_{7})\;\text{and}\;\mathbb{Q}(\zeta_{11})$Qp

We note that if $p$ is a prime number and $E/\mathbb{Q}$ an elliptic curve
such that $E(\Q(\zeta_{p}))_{tors}$ contains a subgroup isomorphic to
$C_{n}\oplus C_{mn}$, then by the properties of Weil pairing we have
$\Q(\zeta_{n})\subseteq\Q(\zeta_{p})$ which forces $n\leq 2$.

Let $E/\mathbb{Q}$ be an elliptic curve and let $p\in\\{5,7,11\\}$ be a prime
number. Apart from the groups in Mazur’s theorem, the group
$E(\mathbb{Q}(\zeta_{p}))_{tors}$ can only be isomorphic to one of the
following groups:

  * •

If $p=5$,

$C_{5}\oplus
C_{5}\;(\href{https://www.lmfdb.org/EllipticCurve/Q/550k2}{550k2}),\quad
C_{15}\;(\href{https://www.lmfdb.org/EllipticCurve/Q/50/a/2}{50a2})\quad\text{and}\quad
C_{16}\;(\href{https://www.lmfdb.org/EllipticCurve/Q/15/a/7}{15a7}).$

  * •

$Ifp=7$,

$C_{13}\;\href{https://www.lmfdb.org/EllipticCurve/Q/147/c/2}{(147c2)},\quad
C_{14}\;\href{https://www.lmfdb.org/EllipticCurve/Q/49a1}{(49a1)},$
$C_{18}\;\href{https://www.lmfdb.org/EllipticCurve/Q/14a4}{(14a4)},\quad
C_{2}\oplus
C_{14}\;\href{https://www.lmfdb.org/EllipticCurve/Q/49a4}{(49a4)},\quad
C_{2}\oplus
C_{18}\;\href{https://www.lmfdb.org/EllipticCurve/Q/14a5}{(14a5)}.$

  * •

If $p=11$,

$C_{11}\;\href{https://www.lmfdb.org/EllipticCurve/Q/121b2}{(121b2)},\quad
C_{25}\;\href{https://www.lmfdb.org/EllipticCurve/Q/11a3}{(11a3)},\quad
C_{2}\oplus
C_{10}\;\href{https://www.lmfdb.org/EllipticCurve/Q/10230bg2}{(10230bg2)}.$

Assume that $p=5$. By [najman_bruin, Theorem 5] we conclude that the only
possibilities are the ones listed in this Lemma and $C_{17}$. It is easy to
rule out $C_{17}$ by using [growth, Theorem 5.8].

Consider the case when $p=11$. If $E/\mathbb{Q}$ is an elliptic curve and if
$C_{n}\oplus C_{mn}\in E(\mathbb{Q}(\zeta_{11}))$, then by the properties of
the Weil pairing we have
$\mathbb{Q}(\zeta_{n})\subseteq\mathbb{Q}(\zeta_{11})$, which forces
$n\in\\{1,2,11\\}$. Applying the Theorem 2 we eliminate the possibility
$n=11$. By [tomi3, Lemma 6.0.8], it remains to show that the groups
$C_{15},C_{16}$ and $C_{2}\oplus C_{12}$ do not occur. We note that the set
$\Phi_{\mathbb{Q}}(2)$ is described by Theorem 1 and the description of the
set $\Phi_{\mathbb{Q}}(5)$ can be found in [38, Theorem 1].

  * •

$C_{15}$: Assume that $P_{15}\in E(\mathbb{Q}(\zeta_{11}))$ is a point of
order $15$. Obviously we have

$E(\mathbb{Q}(\zeta_{11}))[15]\cong C_{15}.$

By Lemma 2 it follows that $[\mathbb{Q}(P_{15}):\mathbb{Q}]\in\\{1,2\\}$, but
this contradicts the Theorem 1.

  * •

$C_{16}$: If $P_{16}\in E(\mathbb{Q}(\zeta_{11}))$ is a point of order $16$,
then $8P_{16}$ has order $2$ and is defined over at most quadratic extension
contained in $\mathbb{Q}(\zeta_{11})$, which is $\mathbb{Q}(\sqrt{-11})$. By
[growth, Proposition 4.8] it follows that
$[\mathbb{Q}(P_{16}):\mathbb{Q}]\in\\{1,2\\}.$ A computation in Magma [3]
shows that $X_{1}(16)(\mathbb{Q}(\sqrt{-11})$ contains only cusps.

  * •

$C_{2}\oplus C_{12}$: As in the previous case, we show that $C_{2}\oplus
C_{12}\subseteq E(\mathbb{Q}(\sqrt{-11}))$. The modular curve
$X_{1}(2,12)(\mathbb{Q}(\sqrt{-11}))$ has rank $0$ and same torsion as over
$\mathbb{Q}$, which means that there does not exist an elliptic curve with
$C_{2}\oplus C_{12}$ torsion over $\mathbb{Q}(\sqrt{-11})$.

It remains to consider $p=7$. Assume that
$n\in\\{11,15,16,17,21,25,19,37,43,67,163,27\\}$ and that
$E(\Q(\zeta_{7}))\cong C_{n}$.

  * •

$n\in\\{11,15,17,25\\}$: Lemma 2 gives us that if $P_{n}$ is a point of order
$n$, we have $[\Q(P_{n}):\Q]\mid 2$. Now Theorem 1 gives us the contradiction.

  * •

$n\in\\{19,37,43,67,163\\}$: From [growth, Theorem 5.8], we get that the point
of order $n$ cannot be defined over the field $\Q(\zeta_{7})$ (a degree $6$
extension).

  * •

$n=27$: This follows from [tomi1, Theorem 1.1].

  * •

$n=16$: Lemma 2 gives us that if $P_{16}\in E(\Q(\zeta_{7}))$ is a point of
order $16$, we have $[\Q(P_{16}):\Q]\mid 2$. That means that $P_{16}\in
E(\Q(\sqrt{-7}))$. We can use the similar methods as before in Magma to
consider $X_{1}(16)(\Q(\sqrt{-7}))$ and prove that $E$ cannot have a point of
order $16$ defined over $\Q(\sqrt{-7})$.

  * •

$n=21$: For $C_{21}$, we conclude from [chou2, Lemma 2.7] that $E$ has a
$\Q$-rational $21$-isogeny. There are $4$ elliptic curves (up to
$\overline{\Q}$-isomorphism) with a rational $21$-isogeny (see [2, p.78-80]).
Therefore, we can use the division polynomial method since the elliptic curves
with the same $j$-invariant have identical division polynomials, up to scalar.
We will consider the seventh division polynomials. We can use Magma [3] to
factor those polynomials in the field $\Q(\zeta_{7})$ and see that they have
no zeroes there. Hence, this case is impossible.

It remains to eliminate only three non-cyclic groups. For $C_{2}\oplus
C_{16}$, we can use Lemma 2 to show that if $P_{16}\in E(\Q(\zeta_{7}))$ is of
order $16$, then $[\Q(2P_{16}):\Q]\mid 2$. By [growth, Proposition 4.6], we
can conclude that $[\Q(P_{16}):\Q]\in\\{1,2\\}$. Hence, we have a point of
order $16$ defined over $\Q(\sqrt{-7})$. However, we already proved that this
can’t happen when we considered $C_{16}$. For $C_{2}\oplus C_{10}$ and
$C_{2}\oplus C_{12}$, we consider modular curves $X_{1}(2,10)(\Q(\zeta_{7}))$
and $X_{1}(2,12)(\Q(\zeta_{7}))$ and use Magma to show that they don’t have
non-cuspidal points, which completes the proof.

Ideally, one would like to give a useful description of possible isomorphism
classes of $E/\mathbb{Q}(\zeta_{p})$, where $E/\mathbb{Q}$ is an elliptic
curve and $p$ is a prime number. One can start with the following question
that seems to be out of reach for the authors at the time of writing this
paper.
Let $n\in\\{13,16,18,25\\}$. Under what conditions on the prime number $p$
does there exist an elliptic curve $E/\mathbb{Q}$ with a point $P_{n}\in
E(\mathbb{Q}(\zeta_{p}))$ of order $n$?

###  5.1 Appendix: Images of Mod $p$Appendix Galois representations associated
to elliptic curves over $\mathbb{Q}$Appendix

For each possible known subgroup $G_{E}(p)\subsetneq\GL_{2}(\mathbb{F}_{p})$
where $E/\mathbb{Q}$ is a non-CM elliptic curve and $p$ is a prime, Tables 1
and 2 list in the first and second column the corresponding labels in
Sutherland and Zywina notations, and the following data:

  * •

$d_{v}=[G_{E}(p):G_{E}(p)_{v}]=|G_{E}(p).v|$ for $v\in\mathbb{F}_{p}^{2}$,
$v\neq(0,0)$; equivalently, the degrees of the extensions $\mathbb{Q}(P)$ over
$\mathbb{Q}$ for points $P\in E(\overline{\mathbb{Q}})$ of order $p$.

  * •

$d=|G_{E}(p)|$; equivalently, the degree $\mathbb{Q}(E[p])$ over $\mathbb{Q}$.

Note that Tables 1 and 2 are partially extracted from Table 3 of
[sutherlandgalrep]. The difference is that [sutherlandgalrep, Table 3] only
lists the minimum of $d_{v}$, which is denoted by $d_{1}$ therein.

|  | Sutherland | Zywina | $d_{v}$ | $d$
---|---|---|---|---
| 2Cs | $G_{1}$ | 1 | 1
| 2B | $G_{2}$ | 1 , 2 | 2
| 2Cn | $G_{3}$ | 3 | 3
| 3Cs.1.1 | $H_{1,1}$ | 1 , 2 | 2
| 3Cs | $G_{1}$ | 2 , 4 | 4
| 3B.1.1 | $H_{3,1}$ | 1 , 6 | 6
| 3B.1.2 | $H_{3,2}$ | 2 , 3 | 6
| 3Ns | $G_{2}$ | 4 | 8
| 3B | $G_{3}$ | 2 , 6 | 12
| 3Nn | $G_{4}$ | 8 | 16
| 5Cs.1.1 | $H_{1,1}$ | 1 , 4 | 4
| 5Cs.1.3 | $H_{1,2}$ | 2 , 4 | 4
|  | Sutherland | Zywina | $d_{v}$ | $d$
---|---|---|---|---
| 5Cs.4.1 | $G_{1}$ | 2 , 4 , 8 | 8
| 5Ns.2.1 | $G_{3}$ | 8 , 16 | 16
| 5Cs | $G_{2}$ | 4 , 4 | 16
| 5B.1.1 | $H_{6,1}$ | 1 , 20 | 20
| 5B.1.2 | $H_{5,1}$ | 4 , 5 | 20
| 5B.1.4 | $H_{6,2}$ | 2 , 20 | 20
| 5B.1.3 | $H_{5,2}$ | 4 , 10 | 20
| 5Ns | $G_{4}$ | 8 , 16 | 32
| 5B.4.1 | $G_{6}$ | 2 , 20 | 40
| 5B.4.2 | $G_{5}$ | 4 , 10 | 40
| 5Nn | $G_{7}$ | 24 | 48
| 5B | $G_{8}$ | 4 , 20 | 80
| 5S4 | $G_{9}$ | 24 | 96
|
|  | Sutherland | Zywina | $d_{v}$ | $d$
---|---|---|---|---
| 7Ns.2.1 | $H_{1,1}$ | $6\,,\,9\,,\,18$ | 18
| 7Ns.3.1 | $G_{1}$ | $12\,,\,18$ | 36
| 7B.1.1 | $H_{3,1}$ | $1\,,\,42$ | 42
| 7B.1.3 | $H_{4,1}$ | $6\,,\,7$ | 42
| 7B.1.2 | $H_{5,2}$ | $3\,,\,42$ | 42
| 7B.1.5 | $H_{5,1}$ | $6\,,\,21$ | 42
| 7B.1.6 | $H_{3,2}$ | $2\,,\,21$ | 42
| 7B.1.4 | $H_{4,2}$ | $3\,,\,14$ | 42
| 7Ns | $G_{2}$ | $12\,,\,36$ | 72
| 7B.6.1 | $G_{3}$ | $2\,,\,42$ | 84
| 7B.6.3 | $G_{4}$ | $6\,,\,14$ | 84
| 7B.6.2 | $G_{5}$ | $6\,,\,42$ | 84
|  | Sutherland | Zywina | $d_{v}$ | $d$
---|---|---|---|---
| 7Nn | $G_{6}$ | $48$ | 96
| 7B.2.1 | $H_{7,2}$ | $3\,,\,42$ | 126
| 7B.2.3 | $H_{7,1}$ | $6\,,\,21$ | 126
| 7B | $G_{7}$ | $6\,,\,42$ | 252
| 11B.1.4 | $H_{1,1}$ | 5 , 110 | 110
| 11B.1.5 | $H_{2,1}$ | 5 , 110 | 110
| 11B.1.6 | $H_{2,2}$ | 10 , 55 | 110
| 11B.1.7 | $H_{1,2}$ | 10 , 55 | 110
| 11B.10.4 | $G_{1}$ | 10 , 110 | 220
| 11B.10.5 | $G_{2}$ | 10 , 110 | 220
| 11Nn | $G_{3}$ | 120 | 240
|
Table 1: Possible images $G_{E}(p)\neq\GL_{2}(\mathbb{F}_{p})$, for $p\leq 11$, for non-CM elliptic curves $E/\mathbb{Q}$. |  | Sutherland | Zywina | $d_{v}$ | $d$
---|---|---|---|---
| 13S4 | $G_{7}$ | $72\,,\,96$ | 288
| 13B.3.1 | $H_{5,1}$ | $3\,,\,156$ | 468
| 13B.3.2 | $H_{4,1}$ | $12\,,\,39$ | 468
| 13B.3.4 | $H_{5,2}$ | $6\,,\,156$ | 468
| 13B.3.7 | $H_{4,2}$ | $12\,,\,78$ | 468
| 13B.5.1 | $G_{2}$ | $4\,,\,156$ | 624
| 13B.5.2 | $G_{1}$ | $12\,,\,52$ | 624
| 13B.5.4 | $G_{3}$ | $12\,,\,156$ | 624
|  | Sutherland | Zywina | $d_{v}$ | $d$
---|---|---|---|---
| 13B.4.1 | $G_{5}$ | $6\,,\,156$ | 936
| 13B.4.2 | $G_{4}$ | $12\,,\,78$ | 936
| 13B | $G_{6}$ | $12\,,\,156$ | 1872
| 17B.4.2 | $G_{1}$ | $8\,,\,272$ | 1088
| 17B.4.6 | $G_{2}$ | $16\,,\,136$ | 1088
| 37B.8.1 | $G_{1}$ | 12 , 1332 | 15984
| 37B.8.2 | $G_{2}$ | 36 , 444 | 15984
|
Table 2: Known images $G_{E}(p)\neq\GL_{2}(\mathbb{F}_{p})$, for $p=13,17$ or
$37$, for non-CM elliptic curves $E/\mathbb{Q}$.

The authors gratefully acknowledge support from the QuantiXLie Center of
Excellence, a project co-financed by the Croatian Government and European
Union through the European Regional Development Fund - the Competitiveness and
Cohesion Operational Programme (Grant KK.01.1.1.01.0004) and by the Croatian
Science Foundation under the project no. IP-2018-01-1313.

## References

  * namelastnametitle __jtitleetal __journalvolumeISBNISSNurlnumeral
  * [1]
  * [2] B.J. Birch W. Kuyk () Modular functions of one variable IV, 476\.  Springer Berlin, Heidelberg, 1975.
  * [3] W. Bosma, J. Cannon C. Playoust The Magma Algebra System I: The User LanguageThe Magma Algebra System I: The User Language.  J. Symbolic Comput., 24(3-4):235–265, 1997


# Optical Character Recognition, Using K-Nearest Neighbors

Wei Wang

###### Abstract

The problem of optical character recognition, OCR, has been widely discussed
in the literature. Having a hand-written text, the program aims at recognizing
the text. Even though there are several approaches to this issue, it is still
an open problem. In this paper we would like to propose an approach that uses
K-nearest neighbors algorithm, and has the accuracy of more than 90%. The
training and run time is also very short.

##  1 Introduction

Visual pattern recognition has been a point of interest of machine learning
experts for several years. A common approach of recognizing a pattern consists
of two important tasks, feature extraction and classification. As a
preprocessing task, feature extraction must be done in order to find a pattern
for the image. Obviously, there is no way to choose best features for a
specific problem. It is usually done by some heuristics that means some
different features, based on experience and problem especific knowledge, are
tested on the same problem and the most efficient one is chosen.
As a start point for solving the problem, a list of reference objects is
required in order to match any new intance to reference images. The result of
the matching task is used to recognize the new object. Although finding
features of an object is very helpful in recognition process, it usualy
requires lots of computations and has its own complexity. Thus it is a
challenging task to choose features useful and at the same time easy to be
computed.
One of the most important reasons of choosing features instead of using pixels
of an image themselves, is curse of dimensionality [1]. It is considerable
that the number of features must not be so high that increases the number of
dimnesions.
Optical Character Recognition (OCR) is a field of research in artificial
intelligence, computer vision and patern recognition that is the recognition
of characters from optical data. The characters may be machine printed or
handwritten [2]. One of the aspect of OCR concerns digit recognition. Digit
recognition is simpler than OCR, because of less number of possibilities and
simpler format of combination of symbols. In this paper what we explain and
try to solve is handwritten digit recognittion.
Though the problem of digit recognition has been intensively investigated, the
improvement of its performance is still a major issue in a number of
industrial applications, e.g. parcel sorting [3]. We here describe our
approaches that uses simple methods to extract features and classify diferent
digits. We developed two simillar algorithms. The first approach is a simple
feature extraction and classification approach that works based on
partitioning the image and averaging the partitions. The second one is very
similar to the former one but gradient based.
The remainder of this report is organized as follows. In section two we
explain two approaches that we developed and tested. In section III we show
the reaults of both methods according to accuracy and timing of free run
tests. The last section is the conclusion of our work.

##  2 Our Approaches

Often vision systems tend to become quite complex. We here propose a simple
two-stage, and another simple one-stage method of feature extraction. We first
start with two-stage method.
The first step of feature extraction is skeletonization. Morphological
skeleton is a skeleton representation of a shape or binary image, computed by
morphological operators [5]. Applying skeletonization, we get rid of some
shadows that might disturb our approach for recognition. This approach is done
by Matlab using bwmorph recursively until there is no difference between the
last two produced images. In figure 1 the skeletonization of a digit is shown.
The second step of our approach is longer than the first introductory one.
Receiving a thin image from the first step, we try to extract features from
it. The first task is to extract the effective part of image from its whole.
The image might have some useless area. This area does not store any
information. This is usually a totally blank part of the image. So it is
better to remove these parts as soon as possible. We name the area of
meaningful data as the effective part of the image. Extracting this area, we
compute the following equations to get minimum and maximum of X and Y,
positions of pixels. Then based on these values we make a rectangle of
efficient area, and extract features only for this area.

$\displaystyle min_{x}=min(X,:)\mbox{ where value(X,Y) $>$ 0}$ $\displaystyle
min_{y}=min(:,Y)\mbox{ where value(X,Y) $>$ 0}$ $\displaystyle
max_{x}=max(X,:)\mbox{ where value(X,Y) $>$ 0}$ $\displaystyle
max_{y}=max(:,Y)\mbox{ where value(X,Y) $>$ 0}$

Extracting features is the most important part of digit recognition. There are
two considerable points in feature extraction. The first one is to avoid many
number of features. Otherwise solving curse of dimensionality is neccessary.
The second one is computational complexity. If feature extraction needs lots
of computational power, it takes a long time to compute features for all
instances. That is a very good reason for choosing simpler mathematical
formulas for feature extraction. We tried two of the simplest forms of feature
extraction to reach to a good timing and accuracy.

Figure 1: Left: raw image, Right: Skeletonized and partitioned image

In the first approach we split the effective image to several rectangular
sections and compute the average value of pixels for each rectangle as the
representative of that section. Amongst different partitioning types, we found
32 sections more efficient than others. In this case we partitioned an
effective image into 4 vertical and 8 horizontal rectangles. The results of
different methods of partitioning is shown in the next section.
In the second approach, we do not do the skeletonization at all. Because
gradient is done to find the highest intensity of black amongst different
pixels of a section. Applying skeletonization makes the gradient useless. We
partitioned the image exactly the same as the previous one. Instead of
computing the average value, we measured the gradient of the pixels in each
section. Applying gradient to a section results in 2D list of derivatives. We
choose the maximum absolute value of derivatives with respect to X as well as
Y. Thus, having the same number of sections, number of features is twice of
the previous algorithm.

When feature extraction is done, we put the array of vectors of features in a
K-Dimentional tree [4]. The kd-tree that we used was an implemented package,
and we utilized it for our purpose. During the classification task, we extract
first three neighbors of the input image using the kd-tree. Given a kd-tree
the algorithm computes a k-nearest neighbor query (kNN) [7, 8] with a
preprocessing time of $O(d\times N\log{N})$ and an expected query time of
$O(\log{N})$ where N is number of points and d dimensionality of a point in
the set.

##  3 Experiment and Results

We tested our algorithms of digit recognition on a set of 2000 instances of
digit images that can be found [6]. Each image was 2D, gray-scale by the size
of 256X256 pixels. We trained our datamodel on 1500 samples of this dataset
and tested it on the rest. It is notable that the test cases contain 50
intances of each class. The computer that we used, had a Core 2 Duo processor
1.66 GHz associated by 2 GBytes of RAM. We tested all of our cases on the same
computer and showed the results in table 1. Running software was Matlab 2008a
with image processing package enabled. Image processing toolbox is required
for skeletonization purpose. The Kd-tree used in this algorithm, is developed
under C++ languague that needs mex and C++ compiler. We installed Microsft
Visual C 2008 express edition for the compiling reason.

Name | Accuracy(percentage) | Run Time(seconds)
---|---|---
4 vertical, 4 horizontal | 87.8 | 0.795
8 vertical, 4 horizontal | 86 | 1.341
4 vertical, 8 horizontal | 92.6 | 1.092
8 vertical, 8 horizontal | 91.8 | 2.152
Gradient Based |  |
4 vertical, 8 horizontal | 80.6 | 4.72
Table 1: Performances of different digit recognition methods

From the results it is easily understandable that the best performance was
obtained by using 32 features when we partitioned it like figure 1.

##  4 Conclusions

In this report we explained our simple approaches of digit detection based on
partitioning an image into rectangles, and using either the average pixels of
sections or their maximum absolute gradients as representatives of the image,
in the other words features of the image. We used kd-tree to find the nearest
neighbors of an image. We showed that despite its simplicity, averaging pixel
values is very efficient, esp in timing issues.

## References

  * [1] Powell Warren B., ”Approximate Dynamic Programming: Solving the Curses of Dimensionality”, Wiley 2007.
  * [2] Charles C. Tapper, Ching Y. Suen and Toru Wakahara ”The State of the Art in On-line Handwriting Recognition”, IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol 12 No 8, August 1990
  * [3] Kai Labusch , Erhardt Barth and Thomas Martinetz , ”Simple method for High-Performance Digit Recognition Based on Sparse Coding”, IEEE transactions on Neural Networks, Vol 19, No 11, November 2008
  * [4] John L. Bentley: ”Multidimensional binary search trees used for associative searching” Comm. ACM 18, 509-517 (1975)
  * [5] www.wikipedia.org as on October 13 2009
  * [6] http://www.faculty.iu-bremen.de/hjaeger/courses/MLFall09/AML_MP1.zip as on October 13 2009
  * [7] Hamraz, Seyed Hamid, and Seyed Shams Feyzabadi. ”General-Purpose learning machine using k-nearest neighbors algorithm.” RoboCup 2005: Robot Soccer World Cup IX. Springer Berlin Heidelberg, 2006. 529-536.
  * [8] Nuchter, Andreas, Seyedshams Feyzabadi, Deyuan Qiu, and Stefan May. ”SLAM à la carte-GPGPU for Globally Consistent Scan Matching.” In ECMR, pp. 271-276. 2011.