# On the Kepler problem on the Heisenberg group

Sergey Basalaev, Sergei Agapov  The work of the second author was performed
according to the Government research assignment for IM SB RAS, project
FWNF-2022-0004.

###### Abstract

We study the nonholonomic motion of a point particle on the Heisenberg group
around the fixed ‘‘sun’’ whose potential is given by the fundamental solution
of the sub-Laplacian. We find three independent first integrals of the system
and show that its bounded trajectories of the system are wound up around
certain surfaces of the fourth order.

_Keywords:_ Heisenberg group, Kepler problem, nonholonomic dynamics, almost
Poisson bracket, first integral.

_2010 Mathematics Subject Classification:_ 37N05, 53C17, 70F25, 37J60.

##  1 Introduction

How would a planet move around the Sun on the Heisenberg group? While studying
the problem we found that the authors of the paper [MS15] aim to answer this
very question. However, a known feature of nonholonomic mechanics (see e. g.
[B03]) is that the variatonal problem (control, geodesics, how to move from A
to B) and the dynamics problem (how does it move on its own?) are generally
not equivalent. Indeed, for instance, in the geodesic problem on the
Heisenberg group, to any initial velocity there corresponds a one-parametric
family of geodesics. On the other hand the dynamics is uniquely determined by
an initial position and a velocity so it can’t be _any_ geodesic (the actual
solutions in in this case are what is known in nonholonomic geometry as the
‘‘straightest’’ lines). It seems to us that the problem actually studied in
[MS15] is the variatonal one — how to move efficiently on Heisenberg group in
the presence of a gravitational field. Here, we aim to solve the dynamics
problem instead.

We consider the Heisenberg group with the left-invariant sub-Riemannian metric
and a fixed ‘‘sun’’ at the origin. The potential is given by the fundamental
solution of the sub-Laplacian — a generalization of the Laplace–Beltrami
operator to the sub-Riemannian manifolds. Traditionally, to derive the non-
holonomic equations of motion the Lagrange–d’Alembert principle is used. In
Section 2 we remind how the equations of motion can be translated to the form
that uses the intrinsic structure of nonholonomic distribution. This allows
one to use Hamiltonian language best suited for finding integrals of the
system. In Section 3 we apply this to study the Kepler problem on the
Heisenberg group and find its first integrals. In contrast to the
6-dimensional variational problem which is not Liouville integrable (proved in
[SM21]), the dynamics problem is 5-dimensional and turns out to have at least
three independent first integrals. This allows us to rather qualitatively
describe the geometry of trajectories of the system. In particular, a typical
trajectory of the system winds up around the surface of order 4 which we found
explicitly. Due to the nonholonomic constraint the surface in the Heisenberg
group uniquely defines the trajectory by its starting point. We also describe
a few special trajectories.

In relation to our research we note that the Kepler problem on the Riemannian
manifolds was studied extensively starting from works of Lobachevsky [L1835]
in hyperbolic space and Serret [S1860] on the sphere. The survey of related
works in the spaces of constant curvature may be found in [DPM12]. The
aforementioned paper [MS15] has a few followups [DS21, SM21] all of which seem
to address the variational problem.

##  2 Motion on sub-Riemannian manifolds

Here we derive the equations of nonholonomic dynamics in the generalized
Hamiltonian form, simplified for the case considered. The general form may be
found in [B03].

Consider a mechanical system in $\mathbb{R}^{n}$ with $k$ ideal functionally
independent nonintegrable constraints linear in velocities. In Lagrangian
coordinates $q^{i}$, $\dot{q}^{i}$ these can be given by

$\sum_{i=1}^{n}a^{j}_{i}(q)\dot{q}^{i}=0,\quad j=1,\ldots,k.$ (2.1)

Locally the equations (2.1) can be solved to $k$ dependent velocities and
represented in the form

$\dot{q}^{m+j}=\sum_{i=1}^{m}f^{j}_{i}(q)\dot{q}^{i},\quad j=1,\ldots,k,$
(2.2)

where $m=n-k$ and the velocities $\dot{q}^{1},\ldots,\dot{q}^{m}$ are assumed
to be independent.

Recall that for a nonholonomic system with the Lagrangian $L(q,\dot{q},t)$ and
the constraints (2.1) the equations of motion are derived using the
Lagrange–d’Alembert principle (see, e. g. [B03])

$\frac{d}{dt}\frac{\partial L}{\partial\dot{q}^{i}}-\frac{\partial L}{\partial
q^{i}}=\sum_{j=1}^{k}\lambda_{j}a^{j}_{i},\quad i=1,\ldots,n,$ (2.3)

where the Lagrange multipliers $\lambda_{j}$ are determined in such a way that
the trajectory satisfies constraints (2.1).

It may be useful, especially for problems with the constraints of form (2.2),
instead of the Lagrangian coordinates use the ones in the distribution of
admissible velocities. Consider the vector fields

$X_{i}(q)={\partial q_{i}}+\sum_{j=1}^{k}f^{j}_{i}(q){\partial q_{m+j}},\quad
i=1,\ldots,m.$ (2.4)

Then the velocity $\dot{q}$ satisfies the constraints (2.2) iff
$\dot{q}=\sum\limits_{i=1}^{m}\dot{q}^{i}X_{i}$. Introduce the momentum 1-form

$P=\sum_{i=1}^{n}\frac{\partial L}{\partial\dot{q}^{i}}dq^{i}.$ (2.5)

Then we can describe the dynamics by the following generalization of
Euler–Lagrange equations (in what follows we denote the action of 1-form
$\tau$ on the vector field $X$ as $\tau\langle X\rangle$).

######  Proposition 2.1.

The dynamical motion in the system with the Lagrangian $L(q,\dot{q},t)$ and
the constraints (2.2) is described by the system of equations

$\displaystyle\frac{d}{dt}P\langle X_{i}\rangle=X_{i}L,$ $\displaystyle
i=1,\ldots,m,$ (2.6)
$\displaystyle\dot{q}^{m+j}=\sum_{i=1}^{m}f^{j}_{i}(q)\dot{q}^{i},$
$\displaystyle j=1,\ldots,k.$

###### Proof.

Introduce 1-forms of our constraints

$\tau^{j}(q)=dq^{m+j}-\sum_{i=1}^{m}f^{j}_{i}(q)\,dq^{i},\qquad j=1,\ldots,k.$
(2.7)

Then $\tau^{j}\langle X_{i}\rangle=0$ for $i=1,\ldots,m$ and
$\tau^{j}\langle\partial q_{m+l}\rangle=\delta^{j}_{l}$ for $l=1,\ldots,k$.
The Lagrange–d’Alembert equations (2.3) can be rewritten in our terms as

$\frac{d}{dt}P\langle{\partial q_{i}}\rangle-dL\langle{\partial
q_{i}}\rangle=\sum_{j=1}^{k}\lambda_{j}\tau^{j}\langle{\partial
q_{i}}\rangle,\qquad i=1,\ldots,n.$

Then, since the expression is linear w. r. t. the term in the angle brackets

$\frac{d}{dt}P\langle X_{i}\rangle-dL\langle
X_{i}\rangle=\sum_{j=1}^{k}\lambda_{j}\tau^{j}\langle X_{i}\rangle=0,\qquad
i=1,\ldots,m.$

The sufficiency of the equations (2.2), (2.6) follows from the fact that we
can recover Lagrange–d’Alembert equations from them. Indeed, let

$\lambda_{j}=\frac{d}{dt}P\langle{\partial q_{m+j}}\rangle-dL\langle{\partial
q_{m+j}}\rangle=\frac{d}{dt}\frac{\partial
L}{\partial\dot{q}^{m+j}}-\frac{\partial L}{\partial q^{m+j}},\quad
j=1,\ldots,k.$

This gives us the equations (2.3) for $i=m+1,\ldots,n$. Then, since $\partial
q_{i}=X_{i}-\sum\limits_{j=1}^{k}f^{j}_{i}(q){\partial q_{m+j}}$ for
$i=1,\ldots,m$ we have

$\frac{d}{dt}\frac{\partial L}{\partial\dot{q}^{i}}-\frac{\partial L}{\partial
q^{i}}=\frac{d}{dt}P\langle X_{i}\rangle-dL\langle
X_{i}\rangle-\sum_{j=1}^{k}f^{j}_{i}\Big{(}\frac{d}{dt}\frac{\partial
L}{\partial\dot{q}^{m+j}}-\frac{\partial L}{\partial
q^{m+j}}\Big{)}=-\sum_{j=1}^{k}\lambda_{j}f^{j}_{i}.$

These are the equations (2.3) for $i=1,\ldots,m$. Thus, the system of
equations (2.2), (2.6) is equivalent to the one of (2.2), (2.3). ∎

Let $\mathcal{D}$ be the distribution spanned by $X_{1},\ldots,X_{m}$, i. e.
$\mathcal{D}_{q}=\mathrm{span}\,\\{X_{1}(q),\ldots,X_{m}(q)\\}$. One thing to
note is that for deriving equations (2.6) for a particular system it is enough
to know the Lagrangian only on $\mathcal{D}$, not on the whole
$T\mathbb{R}^{n}$, which allows us to immerse the problem in the sub-
Riemannian setting.

Recall that the (regular) sub-Riemannian structure on a smooth manifold $M$ is
given by the constant rank distribution $\mathcal{D}\subset TM$ (i. e.
$\mathcal{D}_{x}\subset T_{x}M$ is a subspace and $\dim\mathcal{D}_{x}$ is
independent of $x$) and the sub-Riemannian metric tensor
$\langle\cdot,\cdot\rangle$ on $\mathcal{D}$, i. e.
$\langle\cdot,\cdot\rangle_{x}$ is a scalar product on $\mathcal{D}_{x}$.

Let us reformulate the problem in Hamiltonian terms. The _energy_
$E(q,\dot{q},t)$ of the system is defined as usual:

$E=P\langle\dot{q}\rangle-L$

and satisfies $\frac{d}{dt}E=-\frac{\partial}{\partial t}L$ on trajectories of
the system. Note, that the dual basis to the one of vector fields
$X_{1},\ldots,X_{m},\partial q_{m+1},\ldots,\partial q_{m+k}$ consists of
1-forms

$dq^{1},\ldots,dq^{m},\tau^{1},\ldots,\tau^{k}.$

In particular, $dq^{1},\ldots,dq^{m}$ form the basis of $\mathcal{D}^{*}$.
Introduce the momenta $p=\sum\limits_{i=1}^{m}p_{i}\,dq^{i}$ on
$\mathcal{D}^{*}$, i. e. $p\langle X_{i}\rangle=p_{i}$, $i=1,\ldots,m$.
Assuming that $\dot{q}(p,q,t)\in\mathcal{D}_{q}$ can be determined uniquely
from the equation $p\langle\dot{q}\rangle=P\langle\dot{q}\rangle$ let us
define the _generalized Hamiltonian_ $H(p,q,t)$ on $\mathcal{D}^{*}$ as

$H(p,q,t)=p\langle\dot{q}\rangle-L(q,\dot{q},t)=\sum_{i=1}^{m}p_{i}\dot{q}^{i}-L(q,\dot{q},t).$

For this assumption to take place it is sufficient to require that the
restriction of the quadratic form
$\frac{\partial^{2}L}{\partial\dot{q}^{i}\partial\dot{q}^{j}}dq^{i}dq^{j}$ on
$\mathcal{D}$ is positive definite.

Reformulating the equations (2.6) in terms of $H$ one obtains

######  Proposition 2.2.

The dynamical motion in the nonholonomic system with the constraints (2.2) and
the generalized Hamiltonian $H(p,q,t)$ on $\mathcal{D}^{*}$ is described by
the system of equations

$\displaystyle\dot{q}^{i}=\frac{\partial H}{\partial
p_{i}},\quad\dot{p}_{i}=-X_{i}H,$ $\displaystyle i=1,\ldots,m,$ (2.8)
$\displaystyle\dot{q}^{m+j}=\sum_{i=1}^{m}f^{j}_{i}(q)\frac{\partial
H}{\partial p_{i}},$ $\displaystyle j=1,\ldots,k.$

###### Proof.

Observe that since the bases $X_{1},\ldots,X_{m},\partial
q^{m+1},\ldots,\partial^{m+k}$, introduced in (2.4), and
$d\dot{q}^{1},\ldots,d\dot{q}^{m},\tau^{1},\ldots,\tau^{k}$, introduced in
(2.7), are dual, for a smooth function $f(q)$ we have

$df=\sum_{i=1}^{m}df\langle
X_{i}\rangle\,dq^{i}+\sum_{j=1}^{k}df\langle\partial
q^{m+j}\rangle\,\tau^{j}=\sum_{i=1}^{m}X_{i}f\,dq^{i}+\sum_{j=1}^{k}\frac{\partial
f}{\partial q^{m+j}}\,\tau^{j}.$

Then

$dL=\sum_{i=1}^{m}X_{i}L\,dq^{i}+\sum_{j=1}^{k}\frac{\partial L}{\partial
q^{m+j}}\,\tau^{j}+\sum_{i=1}^{n}\frac{\partial
L}{\partial\dot{q}^{i}}d\dot{q}^{i}+\frac{\partial L}{\partial t}dt\\\
=\sum_{i=1}^{m}X_{i}L\,dq^{i}+\sum_{j=1}^{k}\frac{\partial L}{\partial
q^{m+j}}\,\tau^{j}+\sum_{i=1}^{n}P\langle\partial q^{i}\rangle
d\dot{q}^{i}+\frac{\partial L}{\partial t}dt.$

Further, restricting $L$ on $\mathcal{D}\times\mathbb{R}$ with coordinates
$q^{1},\ldots,q^{n}$, $\dot{q}^{1},\ldots,\dot{q}^{m}$, $t$, i. e. setting
$\dot{q}=\dot{q}^{1}X_{1}+\ldots+\dot{q}^{m}X_{m}$, we obtain

$dL|_{T\mathcal{D}}=\sum_{i=1}^{m}X_{i}L\,dq^{i}+\sum_{j=1}^{k}\frac{\partial
L}{\partial q^{m+j}}\,\tau^{j}+\sum_{i=1}^{m}P\langle X_{i}\rangle
d\dot{q}^{i}+\frac{\partial L}{\partial t}dt.$

Now, for $H=p\langle\dot{q}\rangle-L(q,\dot{q},t)|_{\mathcal{D}}$ with
$\dot{q}=\dot{q}(p,q,t)$ we have

$dH=\sum_{i=1}^{m}\big{(}p_{i}-P\langle
X_{i}\rangle\big{)}\,d\dot{q}^{i}+\sum_{i=1}^{m}\dot{q}^{i}\,dp_{i}-\sum_{i=1}^{m}X_{i}L\,dq^{i}-\sum_{j=1}^{k}\frac{\partial
L}{\partial q^{m+j}}\,\tau^{j}-\frac{\partial L}{\partial t}dt.$

If $\dot{q}$ satisfies $p\langle\dot{q}\rangle=P\langle\dot{q}\rangle$ then
the first sum vanishes. Finally, since for the function $H(p,q,t)$

$dH=\sum_{i=1}^{m}\frac{\partial H}{\partial
p_{i}}dp_{i}+\sum_{i=1}^{m}X_{i}H\,dq^{i}+\sum_{j=1}^{k}\frac{\partial
H}{\partial q^{m+j}}\tau^{j}+\frac{\partial H}{\partial t}dt,$

the equations (2.8) follow from (2.6). ∎

Observe that for time independent system $H$ is a first integral of (2.8),
(2.2). We can define _almost Poisson bracket_ (see e. g. [B03, Section 3.1])
on $\mathcal{D}^{*}$ as

$\\{F,H\\}=\sum_{i=1}^{m}\Big{(}\frac{\partial H}{\partial
p_{i}}X_{i}F-\frac{\partial F}{\partial p_{i}}X_{i}H\Big{)}.$

It retains all the properties of Poisson bracket but the Jacobi identity.
Nevertheless, since $\dot{F}=\\{F,H\\}$ it follows that $F$ is an integral of
(2.8), (2.2) iff $\\{F,H\\}\equiv 0$.

##  3 Motion in a potential field on the Heisenberg group

Recall that the Heisenberg group
$\mathbb{H}^{1}=(\mathbb{R}^{3},\cdot,\delta_{\lambda})$ is a homogeneous
group with the group operation

$(x,y,z)\cdot(x^{\prime},y^{\prime},z^{\prime})=\Big{(}x+x^{\prime},y+y^{\prime},z+z^{\prime}+\frac{xy^{\prime}-x^{\prime}y}{2}\Big{)},$

and the one-parametric family of anisotropic dilatations

$\delta_{\lambda}(x,y,z)=(\lambda x,\lambda y,\lambda^{2}z),\quad\lambda>0.$

Its Lie algebra $\mathfrak{h}^{1}$ of left-invariant vector fields has the
basis

$X=\partial_{x}-\frac{y}{2}\partial_{z},\quad
Y=\partial_{y}+\frac{x}{2}\partial_{z},\quad Z=[X,Y]=\partial_{z}.$

The dual basis of left-invariant 1-forms is

$dx,\quad dy,\quad\tau=dz+\frac{y\,dx-x\,dy}{2}.$

The _horizontal distribution_
$\mathcal{D}=\mathop{\mathrm{span}}\\{X,Y\\}\subset T\mathbb{H}^{1}$ is
totally nonholonomic. The form $\tau$ is its annihilator. The sub-Riemannian
structure on $\mathbb{H}^{1}$ is given by the quadratic form
$\langle\cdot,\cdot\rangle$ on $\mathcal{D}$. We choose the one such that
$X,Y$ form the orthonormal basis:

$ds^{2}=dx^{2}+dy^{2}.$

While this quadratic form is degenerate on $T\mathbb{H}^{1}$ it is positive
definite on $\mathcal{D}$. For the mechanical motion with the kinetic energy
$T=\frac{1}{2}ds^{2}\langle\dot{q}\rangle=\frac{1}{2}(\dot{x}^{2}+\dot{y}^{2})$
and the potential energy $U=U(x,y,z)$ one has, as usual, the Lagrangian
$L=T-U$. By Proposition 2.2 we can translate equations to the Hamiltonian form
where the Hamiltonian $H$ on $\mathcal{D}^{*}$ takes the form $H=2T-L=T+U$, i.
e.

$H(x,y,z,p_{X},p_{Y})=\frac{p_{X}^{2}+p_{Y}^{2}}{2}+U(x,y,z).$

We are interested in the gravitational potential which in $\mathbb{R}^{n}$ is
given by a fundamental solution of the Laplacian. The analogue of
Laplace–Beltrami operator on the Heisenberg group is the operator
$\Delta_{H}=X^{2}+Y^{2}$. Its fundamental solution111In the cited paper [MS15]
the potential $U$ has the term $\frac{z^{2}}{16}$ instead of $16z^{2}$. One
can check that $16$ is the correct coefficient since only in this case
$\Delta_{H}U=0$ away from the origin. is found in [F79]:

$U=-\frac{k}{\rho^{2}},\quad\text{where
}\rho(x,y,z)=((x^{2}+y^{2})^{2}+16\\*z^{2})^{\frac{1}{4}}$

and $k>0$ is some constant. Since both the distribution and the potential have
a rotational symmetry around $Oz$ it is natural to make the cylindrical
coordinate change $x=r\cos\theta,y=r\sin\theta$. The basis of $\mathcal{D}$
may be given by vector fields

$\displaystyle R$
$\displaystyle=\phantom{-r}\cos\theta\,X+\phantom{r}\sin\theta\,Y=\partial_{r},$
$\displaystyle S$
$\displaystyle=-r\sin\theta\,X+r\cos\theta\,Y=\partial_{\theta}+\tfrac{r^{2}}{2}\partial_{z}.$

Duals to the basis $R,S,\partial_{z}$ are
$dr,d\theta,\tau=dz-\frac{r^{2}}{2}d\theta$ and for the momenta we have

$p_{X}\,dx+p_{Y}\,dy=(p_{X}\cos\theta+p_{Y}\sin\theta)dr+r(p_{Y}\,\cos\theta-
p_{X}\,\sin\theta)d\theta=p_{R}\,dr+p_{S}\,d\theta.$

It follows that $T=\frac{p_{R}^{2}+p_{S}^{2}/r^{2}}{2}$ and the Hamiltonian
becomes

$H(r,\theta,z,p_{R},p_{S})=\frac{p_{R}^{2}+\frac{1}{r^{2}}p_{S}^{2}}{2}-\frac{k}{(r^{4}+16z^{2})^{\frac{1}{2}}}.$

Since constraints in the new coordinates still have the form (2.2) we may
apply Proposition 2.2 to derive the equations of motion:

$\displaystyle\dot{r}$ $\displaystyle=\frac{\partial H}{\partial
p_{R}}=p_{R},$ $\displaystyle\dot{p}_{R}$
$\displaystyle=-RH=\frac{p_{S}^{2}}{r^{3}}-\frac{2kr^{3}}{(r^{4}+16z^{2})^{\frac{3}{2}}},$
$\displaystyle\dot{\theta}$ $\displaystyle=\frac{\partial H}{\partial
p_{S}}=\frac{p_{S}}{r^{2}},$ $\displaystyle\dot{p}_{S}$
$\displaystyle=-SH=-\frac{8kr^{2}z}{(r^{4}+16z^{2})^{\frac{3}{2}}},$
$\displaystyle\dot{z}$ $\displaystyle=\frac{r^{2}}{2}\frac{\partial
H}{\partial p_{S}}=\frac{p_{S}}{2}.$ (3.1)

######  Theorem 3.1.

If $H<0$ the solutions of (3.1) are bounded with
$\sqrt{r^{4}+16z^{2}}\leq\frac{k}{|H|}$.

###### Proof.

This easily follows from the inequality
$H+k/\sqrt{r^{4}+16z^{2}}=\frac{p_{R}^{2}+\frac{1}{r^{2}}p_{S}^{2}}{2}\geq 0$.
∎

In what follows we search for the additional first integrals of the system.

######  Proposition 3.2.

The system (3.1) does not admit any linear in momenta first integrals.

This statement can be checked by straightforward calculations. We skip the
details. However, it turns out that there are a few quadratic integrals in
addition to the Hamiltonian $H$.

######  Theorem 3.3.

The system (3.1) admits quadratic in momenta first integrals

$\displaystyle F_{1}$
$\displaystyle=\phantom{-}\Big{(}p_{R}p_{S}r-2p_{R}^{2}z+\frac{2p_{S}^{2}z}{r^{2}}\Big{)}\cos(2\theta)+\Big{(}\frac{4p_{R}p_{S}z}{r}-p_{S}^{2}+\frac{kr^{2}}{\sqrt{r^{4}+16z^{2}}}\Big{)}\sin(2\theta),$
$\displaystyle F_{2}$
$\displaystyle=-\Big{(}p_{R}p_{S}r-2p_{R}^{2}z+\frac{2p_{S}^{2}z}{r^{2}}\Big{)}\sin(2\theta)+\Big{(}\frac{4p_{R}p_{S}z}{r}-p_{S}^{2}+\frac{kr^{2}}{\sqrt{r^{4}+16z^{2}}}\Big{)}\cos(2\theta),$
$\displaystyle F_{3}$
$\displaystyle=(2zp_{R}-rp_{S})^{2}+4z^{2}\Big{(}\frac{p_{S}^{2}}{r^{2}}+\frac{2k}{\sqrt{r^{4}+16z^{2}}}\Big{)}.$

Any three of $H,F_{1},F_{2},F_{3}$ are functionally independent a. e. wherein
all of them satisfy the relation

$F_{1}^{2}+F_{2}^{2}=2HF_{3}+k^{2}.$ (3.2)

This theorem can be verified by straightforward calculations. The method we
used to construct these integrals is described in Appendix A.

Knowing three independent first integrals allows us to derive the equation of
the surface (in coordinates $(r,\theta,z)$) in which the trajectories lie. To
do that we introduce two more conserved quantities $J\geq 0$ and in the case
$J>0$ also $\theta_{0}\in[0,\pi)$ such that

$F_{1}=J\sin(2\theta_{0}),\quad F_{2}=J\cos(2\theta_{0}).$

As main parameters we choose $H,F_{3}$ that do not depend on the angle and the
phase offset $\theta_{0}$ that captures the rotational symmetry of the
problem. Note, that from the definition $F_{3}\geq 0$ since $k>0$, and
$2HF_{3}+k^{2}=J^{2}\geq 0$. Therefore we have one general case and two cases
that might require special handling:

  * •

The general case $F_{3}>0$, $J>0$. Then $H>-\frac{k^{2}}{2F_{3}}$ and
$\theta_{0}\in[0,\pi)$ is defined.

  * •

The minimum energy case $F_{3}>0$, $J=0$. Then
$H=H_{\min}=-\frac{k^{2}}{2F_{3}}$, $\theta_{0}$ is undefined.

  * •

The degenerate case $F_{3}=0$. In this case $J=k$, $\theta_{0}\in[0,\pi)$ is
defined and $H$ is unbounded.

######  Theorem 3.4.

All trajectories of the system (3.1) with the fixed values of the first
integrals $H,F_{3},\theta_{0}$ lie on the surface which in the general case
$F_{3}>0$, $J>0$ satisfies the equation

$F_{3}=8z^{2}H+k\sqrt{r^{4}+16z^{2}}-\sqrt{k^{2}+2HF_{3}}\,r^{2}\cos\big{(}2(\theta-\theta_{0})\big{)}.$
(3.3)

In the minimum energy case $F_{3}>0$, $J=0$ the surface becomes an ellipsoid
of revolution

$4k^{2}z^{2}+kF_{3}r^{2}=F_{3}^{2}.$ (3.4)

In the degenerate case $F_{3}=0$ the surface degenerates to the straight
horizontal line passing through the origin

$z=0,\quad\theta=\theta_{0}\mod\pi.$ (3.5)

###### Proof.

Observe, that we can write $F_{3}$ as

$\displaystyle
F_{3}=8z^{2}H+k\sqrt{r^{4}+16z^{2}}-r^{2}\Big{(}-p_{S}^{2}+\frac{4p_{R}p_{S}z}{r}+\frac{kr^{2}}{\sqrt{r^{4}+16z^{2}}}\Big{)}.$

From the expressions of $F_{1}$ and $F_{2}$ we have

$\Big{(}-p_{S}^{2}+\frac{4p_{R}p_{S}z}{r}+\frac{kr^{2}}{\sqrt{r^{4}+16z^{2}}}\Big{)}=F_{1}\sin(2\theta)+F_{2}\cos(2\theta).$

Therefore,

$F_{3}=8z^{2}H+k\sqrt{r^{4}+16z^{2}}-r^{2}\big{(}F_{1}\sin(2\theta)+F_{2}\cos(2\theta)\big{)}.$
(3.6)

Let $F_{3}>0$ and $J>0$. We have
$F_{1}\sin(2\theta)+F_{2}\cos(2\theta)=J\cos(2(\theta-\theta_{0}))$ and (3.3)
follows.

Now, let $F_{3}>0$ and $J=0$. In this case $F_{1}=F_{2}=0$ and the last term
in (3.6) vanishes. Since $H=-\frac{k^{2}}{2F_{3}}$ in this case, (3.6) becomes

$F_{3}^{2}+4k^{2}z^{2}=kF_{3}\sqrt{r^{4}+16z^{2}}.$

Squaring and simplifying it we obtain

$(F_{3}^{2}-4k^{2}z^{2})^{2}=k^{2}F_{3}^{2}r^{4}.$ (3.7)

By Theorem 3.1 for the trajectories of the system
$\sqrt{r^{4}+16z^{2}}\leq\frac{k}{|H|}=\frac{2F_{3}}{k}$. Therefore,

$F_{3}\geq\tfrac{k}{2}\sqrt{r^{4}+16z^{2}}\geq\tfrac{k}{2}\sqrt{16z^{2}}=2kz.$

Now (3.4) follows if we take the square root of (3.7).

Lastly, $F_{3}=0$ implies $z=0$ and $J=k>0$. The restriction of (3.6) on $z=0$
becomes

$0=kr^{2}-kr^{2}\cos(2(\theta-\theta_{0})).$

This gives us either $r=0$ (the origin) or $\theta=\theta_{0}\mod\pi$. ∎

######  Remark 3.5.

Solving the equation (3.3) for the square root $\sqrt{r^{4}+z^{2}}$ and then
squaring it we obtain the following equation

$k^{2}(r^{4}+16z^{2})=\big{(}F_{3}-8z^{2}H+\sqrt{k^{2}+2HF_{3}}\,r^{2}\cos\big{(}2(\theta-\theta_{0})\big{)}\big{)}^{2},$

or in Cartesian coordinates

$k^{2}((x^{2}+y^{2})^{2}+16z^{2})=\big{(}F_{3}-8z^{2}H+\sqrt{k^{2}+2HF_{3}}\,(\cos(2\theta_{0})(x^{2}-y^{2})+2\sin(2\theta_{0})xy)\big{)}^{2}.$

We see that this is an equation of the fourth order. However, its solution is
a branched surface and only one of its branches is the solution to the
original equation, i. e. the equation is quadratic in $z^{2}$ but only one of
its two roots solves (3.3).

Examples of the surfaces corresponding to the cases $H=0$ and $H<0$ may be
seen on Fig. 1. Next we note a few properties of the surfaces obtained.

Figure 1: Surfaces corresponding to the cases $H=0$ and $H<0$.

######  Corollary 3.6.

In the non-degenerate case $F_{3}>0$ the surfaces described in Theorem 3.4
have the following properties.

  1. 1.

The surface is topologically

    1. (a)

a sphere in the case $H<0$;

    2. (b)

a cylinder in the case $H\geq 0$.

  2. 2.

The surface has reflection symmetry in the plane $z=0$ and

    1. (a)

in the case $J>0$ it has two more planes of symmetry $\theta=\theta_{0}$ and
$\theta=\theta_{0}+\frac{\pi}{2}$;

    2. (b)

in the case $J=0$ it is the surface of revolution around $r=0$.

  3. 3.

The trace of the surface on the plane $z=0$ is a quadratic curve:

    1. (a)

in the case $H<0$ it is an ellipse with the semiaxes
$\sqrt{\frac{F_{3}}{k-J}}$ and $\sqrt{\frac{F_{3}}{k+J}}$;

    2. (b)

in the case $H=0$ it is two parallel lines at the distance
$\sqrt{\frac{F_{3}}{k}}$ from the origin;

    3. (c)

in the case $H>0$ it is a hyperbola with the semiaxis
$\sqrt{\frac{F_{3}}{k+J}}$.

The properties are straightforward and easy to check.

A smooth surface $S\subset\mathbb{H}^{1}$ is transversal to the horizontal
distribution $\mathcal{D}$ at almost all points, i. e.
$T_{x}S\cap\mathcal{D}_{x}$ is one-dimensional for a. e. $x\in S$. Therefore,
the trajectory of the system is rather uniquely defined by a starting point on
a surface, with the only possible exception being when the solution arrives at
the point tangent to $\mathcal{D}$ with zero velocity. An example of a bounded
trajectory ($H<0$) which we belive to be a typical one is presented in Fig. 2.
Next we describe special solutions corresponding to the degenerate cases.

Figure 2: A trajectory and its projection on the plane $Oxy$ with $k=1$ and
initial values $x(0)=1$, $y(0)=z(0)=p_{X}(0)=0$, $p_{Y}(0)=\frac{1}{10}$.

######  Theorem 3.7.

The only trajectories of (3.1) passing through the origin are straight lines
in the plane $z=0$. In this case $\theta=const$ and $r(t)$ satisfies the
equation $\frac{\dot{r}^{2}}{2}=H+\frac{k}{r^{2}}$. These solutions correspond
to the degenerate case $F_{3}=0$.

###### Proof.

From Corollary 3.6 the trajectories may pass through the origin only in the
case $F_{3}=0$, i. e. only if $z\equiv 0$ and $p_{S}\equiv 0$. Then
$\dot{\theta}=0$ and the Hamiltonian becomes
$H=\frac{\dot{r}^{2}}{2}-\frac{k}{r^{2}}$. ∎

Therefore, the conserved quantity $F_{3}$ serves as a kind of angular/vertical
momentum. Two more special solutions appear in the minimal energy case $J=0$.

######  Theorem 3.8.

The only stationary solutions of (3.1) are points on $Oz$:

$r=0,\quad z=\pm\frac{k}{4H}.$

These solutions correspond to the minimal energy case $J=0$.

###### Proof.

Indeed, outside of the axis $r=0$ the stationary solution must satisfy
$p_{S}\equiv 0$. But in this case $\dot{p}_{R}$ is negative and the solution
is non-stationary. For the stationary solution on $Oz$ we have
$H=-\frac{k}{\sqrt{16z^{2}}}$ and $F_{3}=\frac{8kz^{2}}{\sqrt{16z^{2}}}$.
Therefore $J^{2}=k^{2}+2HF_{3}=0$. ∎

######  Theorem 3.9.

Let $J=0$ and $z_{0}=\frac{k}{4|H|}$. The trajectories of non-stationary
solutions to (3.1) are the curves monotone in $z$ and $\theta$ such that being
parameterized by $z$ they have the form

$r(z)=\Big{(}2\frac{z_{0}^{2}-z^{2}}{z_{0}}\Big{)}^{\frac{1}{2}},\qquad\theta(z)=\frac{1}{2}\log\frac{z_{0}+z}{z_{0}-z}+\theta(0),\qquad|z|\leq
z_{0}.$

These solutions connect the stationary points $(0,0,\pm z_{0})$ and take
infinite time to approach them, i. e. $t(z)\to\pm\infty$ as $z\to\pm z_{0}$.

###### Proof.

Note that $F_{3}\geq 0$. Therefore, $J=\sqrt{k^{2}+2HF_{3}}=0$ implies $H<0$.
Hence, $H=-\frac{k}{4z_{0}}$ and $F_{3}=-\frac{k^{2}}{2H}=2kz_{0}$. The
surface (3.4) in this case is an ellipsoid of revolution

$\frac{z_{0}r^{2}}{2}+\frac{z^{2}}{z_{0}^{2}}=1.$

From this equation we find the dependence $r(z)$. Next, from
$\frac{d\theta}{dz}=\frac{\dot{\theta}}{\dot{z}}=\frac{2}{r^{2}}$ we also find
the expression of $\theta(z)$ in a closed form. This gives us a family of
curves described in the statement of the theorem. Take any such curve. From
the expression of $H$ we have

$\frac{\dot{r}^{2}+\frac{4\dot{z}^{2}}{r^{2}}}{2}=\frac{p_{R}^{2}+\frac{p_{S}^{2}}{r^{2}}}{2}=H+\frac{k}{\sqrt{r^{4}+16z^{2}}}=\frac{k(z_{0}^{2}-z^{2})}{4z_{0}(z_{0}^{2}+z^{2})}>0$
(3.8)

for all points except the poles. Therefore, the velocity along the curve is
nonzero except on the endpoints. Hence the solution restricted to a curve is
monotone in $z$. From the equation of the surface we have

$0=z_{0}r\dot{r}+\frac{2z\dot{z}}{z_{0}^{2}}.$

This together with (3.8) yields the equation on $z$:

$\dot{z}^{2}=\frac{kz_{0}^{4}(z_{0}^{2}-z^{2})^{2}}{4(z^{2}+z_{0}^{2})(z^{2}+z_{0}^{6})}.$

Choosing the solution increasing in $z$ we obtain

$t(z_{1})=\int_{0}^{z_{1}}\frac{dz}{\dot{z}}=\int_{0}^{z_{1}}\frac{2\sqrt{(z^{2}+z_{0}^{2})(z^{2}+z_{0}^{6})}}{\sqrt{k}z_{0}^{2}(z_{0}^{2}-z^{2})}\,dz$

which diverges as $z_{1}\to\pm z_{0}$. The theorem is proved. ∎

##  4 Conclusion

In conclusion we see that the variational problem and the dynamics problem on
the Heisenberg group are vastly different. While the first is Hamiltonian but
non-integrable in Liouville sense, the second one, being non-Hamiltonian has
at least three first integrals. Both problems are interesting but provide a
different insight into the nonholonomic world.

Acknowledgements: While the problem was mainly considered by the first author
he is very grateful to second author for deriving the first integrals
(Proposition 3.2, Theorem 3.3 and Appendix) which significantly progressed the
research and for overall critical comments.

The images were prepared using GNUPlot and Maxima free software.

##  Appendix A Derivation of quadratic first integrals

By definition any first integral $F$ of (3.1) must satisfy the following
relation:

$\frac{dF}{dt}=\frac{\partial F}{\partial r}p_{R}+\frac{\partial
F}{\partial\theta}\frac{p_{S}}{r^{2}}+\frac{\partial F}{\partial
z}\frac{p_{S}}{2}+\frac{\partial F}{\partial
p_{R}}\left(\frac{p_{S}^{2}}{r^{3}}-\frac{2kr^{3}}{(r^{4}+16z^{2})^{3/2}}\right)-\frac{\partial
F}{\partial p_{S}}\frac{8kr^{2}z}{(r^{4}+16z^{2})^{3/2}}=0.$ (A.1)

It is quite natural to search for the first integrals of (3.1) having the form
of non-homogeneous polynomials in momenta.

We shall search for the quadratic integral of (3.1) in the form:

$F=a\,p_{R}^{2}+d\,p_{R}p_{S}+b\,p_{S}^{2}+f\,p_{R}+g\,p_{S}+h,$

where all the coefficients are unknown functions which depend on $r$,
$\theta$, $z$. Writing down the condition (A.1) for such an integral $F$, we
obtain the system of PDEs which splits into two parts: the first one contains
relations between the unknown functions $a$, $b$, $d$, $h$ only, the second
one is between $f$ and $g$. As in Proposition 3.2, it is easy to check that if
$F$ is the first integral, then both functions $f$ and $g$ must vanish
identically. So we start our analysis with an integral of the form

$F=a(r,\theta,z)p_{R}^{2}+d(r,\theta,z)p_{R}p_{S}+b(r,\theta,z)p_{S}^{2}+h(r,\theta,z).$

The condition (A.1) implies:

$\displaystyle a_{r}=0,$ (A.2) $\displaystyle
2a_{\theta}+r^{2}(a_{z}+2d_{r})=0,$ (A.3) $\displaystyle
4a+r^{3}d_{z}+2r(d_{\theta}+r^{2}b_{r})=0,$ (A.4) $\displaystyle
2d+r^{3}b_{z}+2rb_{\theta}=0,$ (A.5) $\displaystyle
2r^{3}(r^{4}+16z^{2})^{3/2}h_{r}-8kr^{6}a-16kr^{5}zd=0,$ (A.6) $\displaystyle
r(r^{4}+16z^{2})^{3/2}(r^{2}h_{z}+2h_{\theta})-4kr^{6}d-32kr^{5}zb=0.$ (A.7)

Integrating the equations (A.2)–(A.4) successively, we obtain

$\displaystyle a(r,\theta,z)=\alpha(\theta,z),\qquad
d(r,\theta,z)=\gamma(\theta,z)-\frac{1}{2}r\alpha_{z}+\frac{\alpha_{\theta}}{r},$
$\displaystyle
b(r,\theta,z)=\frac{\alpha}{r^{2}}+\frac{\alpha_{\theta\theta}}{2r^{2}}+\frac{r^{2}\alpha_{zz}}{8}+\frac{\gamma_{\theta}}{r}-\frac{r\gamma_{z}}{2}+\omega(\theta,z),$

where $\alpha(\theta,z)$, $\gamma(\theta,z)$, $\omega(\theta,z)$ are arbitrary
functions. Then (A.5) takes the form

$\alpha_{zzz}r^{6}-4\gamma_{zz}r^{5}+(2\alpha_{\theta
zz}+8\omega_{z})r^{4}+4(\alpha_{\theta\theta
z}+4\omega_{\theta})r^{2}+16(\gamma_{\theta\theta}+\gamma)r+8(\alpha_{\theta\theta\theta}+\alpha_{\theta})=0.$

This is a polynomial in $r$ with coefficients depending on $\theta$, $z$ only.
Since this polynomial must vanish, all its coefficients must vanish as well.
This allows one to find $\alpha(\theta,z)$, $\gamma(\theta,z)$,
$\omega(\theta,z)$ and, consequently, the coefficients $a(r,\theta,z)$,
$b(r,\theta,z)$, $d(r,\theta,z)$ explicitly. We omit these long but simple
calculations and skip the final form of these coefficients since they are
quite cumbersome.

After that we are left with two equations (A.6), (A.7) on the unknown function
$h(r,\theta,z)$ which take the form:

$\displaystyle(r^{4}+16z^{2})^{3/2}h_{r}+2kr(c_{8}r^{2}-z(4c_{9}+z(4c_{3}+c_{5}r^{2}+4c_{6}z)))\cos(2\theta)$
$\displaystyle-2kr(c_{9}r^{2}+z(4c_{8}+z(4c_{2}-c_{6}r^{2}+4c_{5}z)))\sin(2\theta)-4kr^{3}(c_{7}-c_{4}z^{2})$
$\displaystyle-8kr^{2}z((s_{2}+s_{4}z)\cos\theta+(s_{3}+s_{5}z)\sin\theta)=0,$
(A.8)
$\displaystyle(r^{4}+16z^{2})^{3/2}(r^{2}h_{z}+2h_{\theta})+2c_{1}kr^{2}(r^{4}-16z^{2})$
$\displaystyle-
kr^{2}(4c_{9}r^{2}+c_{2}(r^{4}+16z^{2})-2z(-8c_{8}+2c_{3}r^{2}+c_{5}r^{4}+6c_{6}r^{2}z-8c_{5}z^{2}))\cos(2\theta)$
$\displaystyle+kr^{2}(-4c_{8}r^{2}+c_{3}(r^{4}+16z^{2})+2z(8c_{9}+2c_{2}r^{2}-c_{6}r^{4}+6c_{5}r^{2}z+8c_{6}z^{2}))\sin(2\theta)$
$\displaystyle-4kr^{3}(s_{2}r^{2}+z(8s_{3}-3s_{4}r^{2}+8s_{5}z))\cos\theta+4kr^{3}(-s_{3}r^{2}+z(8s_{2}+3s_{5}r^{2}+8s_{4}z))\sin\theta$
$\displaystyle-4kr^{2}z(8c_{7}+8s_{1}r^{2}+c_{4}(r^{4}+8z^{2}))=0.$ (A.9)

Here $c_{k}$, $s_{k}$ are arbitrary constants. The equation (A.8) can be
integrated. However, the general solution $h(r,\theta,z)$ to (A.8) is
expressed in terms of elliptic integrals. We consider the simplest case

$s_{2}=s_{3}=s_{4}=s_{5}=0.$

In this case $h(x,y,z)$ can be found from (A.8) in terms of elementary
functions as follows:

$h(r,\theta,z)=\psi(\theta,z)+\frac{k}{4z\sqrt{r^{4}+16z^{2}}}\phi(r,\theta,z),$

where $\psi(\theta,z)$ is an arbitrary function and

$\phi(r,\theta,z)=(c_{9}r^{2}+z(4c_{8}+c_{3}r^{2}+c_{6}r^{2}z-4c_{5}z^{2}))\cos(2\theta)\\\
+(c_{8}r^{2}+z(-4c_{9}+c_{2}r^{2}+c_{5}r^{2}z+4c_{6}z^{2}))\sin(2\theta)+8z(c_{4}z^{2}-c_{7}).$

The unknown function $\psi(\theta,z)$ should be chosen such that the relation
(A.9) holds identically. It seems that the only possible way to satisfy this
requirement is to put

$\psi(\theta,z)\equiv 0,\qquad c_{1}=c_{5}=c_{6}=c_{8}=c_{9}=s_{1}=0.$

In this case (A.9) is satisfied. Thus we found all the coefficients of $F.$
Notice that $\widetilde{F}=4F-8c_{7}H$ is also the first integral of (3.1)
having the simpler form:

$\widetilde{F}=\widetilde{a}\,p_{R}^{2}+\widetilde{d}\,p_{R}p_{S}+\widetilde{b}\,p_{S}^{2}+\widetilde{h},$

where

$\displaystyle\widetilde{a}$
$\displaystyle=2z(2c_{4}z-c_{2}\cos(2\theta)+c_{3}\sin(2\theta)),$
$\displaystyle\widetilde{d}$
$\displaystyle=\Big{(}c_{2}r+\frac{4c_{3}z}{r}\Big{)}\cos(2\theta)-\Big{(}c_{3}x-\frac{4c_{2}z}{r}\Big{)}\sin(2\theta)-4c_{4}rz,$
$\displaystyle\widetilde{b}$
$\displaystyle=\frac{1}{r^{2}}\big{(}(-c_{3}r^{2}+2c_{2}z)\cos(2\theta)-(c_{2}r^{2}+2c_{3}z)\sin(2\theta)+c_{4}(r^{4}+4z^{2})\big{)},$
$\displaystyle\widetilde{h}$
$\displaystyle=\frac{k}{\sqrt{r^{4}+16z^{2}}}\big{(}r^{2}(c_{3}\cos(2\theta)+c_{2}\sin(2\theta))+8c_{4}z^{2}\big{)}.$

Here $c_{2}$, $c_{3}$, $c_{4}$ are arbitrary constants. It is left to notice
that $\widetilde{F}$ is linear in these constants, i. e. it has the form
$\widetilde{F}=c_{2}F_{1}+c_{3}F_{2}+c_{4}F_{3}$. This implies that the
functions

$\displaystyle F_{1}$
$\displaystyle=\phantom{-}\Big{(}p_{R}p_{S}r-2p_{R}^{2}z+\frac{2p_{S}^{2}z}{r^{2}}\Big{)}\cos(2\theta)+\Big{(}-p_{S}^{2}+\frac{4p_{R}p_{S}z}{r}+\frac{kr^{2}}{\sqrt{r^{4}+16z^{2}}}\Big{)}\sin(2\theta),$
$\displaystyle F_{2}$
$\displaystyle=-\Big{(}p_{R}p_{S}r-2p_{R}^{2}z+\frac{2p_{S}^{2}z}{r^{2}}\Big{)}\sin(2\theta)+\Big{(}-p_{S}^{2}+\frac{4p_{R}p_{S}z}{r}+\frac{kr^{2}}{\sqrt{r^{4}+16z^{2}}}\Big{)}\cos(2\theta),$
$\displaystyle F_{3}$
$\displaystyle=4z^{2}p_{R}^{2}-4rzp_{R}p_{S}+\frac{r^{4}+4z^{2}}{r^{2}}p_{S}^{2}+\frac{8kz^{2}}{\sqrt{r^{4}+16z^{2}}}.$

are also integrals of (3.1).

## References

  * [B03] A. M. Bloch et al., _Nonholonomic Mechanics and Control_ / Interdisciplinary Applied Math. 24 (2003), 501 pp.
  * [DPM12] F. Diacu, E. Pérez-Chavela, M. Santoprete, _The $n$-body Problem in Spaces of Constant Curvature. Part I: Relative Equilibria_ // J. Nonlinear Sci. 22 (2012), 247–266, DOI: 10.1007/s00332-011-9116-z
  * [DS21] V. Dods, C. Shanbrom, _Self-similarity in the Kepler–Heisenberg Problem_ // J. Nonlinear Sci. 31 (2021), Article 49, DOI: 10.1007/s00332-021-09709-1
  * [F79] G. Folland, _A fundamental solution for a subellipic operator_ // Bulletin of the AMS 79:2 (1973), 373–376.
  * [L1835] N. I. Lobachevsky, _The new foundations of geometry with full theory of parallels_ / Collected Works 1835–1838 2, GITTL, Moscow, 1949, 159 pp. [in Russian]
  * [MS15] R. Montgomery, C. Shanbrom, _Keplerian Dynamics on the Heisenberg Group and Elsewhere_ // Fields Inst. Commun. 73 (2015), 319–342, arXiv:1212.2713 [math.DS]
  * [S1860] P. J. Serret _Théorie nouvelle géométrique et mécanique des lignes a double courbre_ / Librave de Mallet-Bachelier, Paris, 1860.
  * [SM21] T. Stachowiak, A. J. Maciejewski, _Non-Integrability of the Kepler and the Two-Body Problems on the Heisenberg Group_ // SIGMA 17 (2021), Article 074, arXiv:2103.10495 [math-ph], DOI: 10.3842/SIGMA.2021.074

Sergey Basalaev
Novosibirsk State University, 1 Pirogova st., 630090 Novosibirsk Russia
e-mail<EMAIL_ADDRESS>
Sergei Agapov
Sobolev Institute of Mathematics, 4 Acad. Koptyug avenue, 630090 Novosibirsk
Russia
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>


points $x\in\mathrm{int}(K)$ and $x^{\prime}\in\mathrm{int}(K)$ with non-
vanishing components. Note that
$\hat{L}_{\mathrm{lex}}(q)(\uplambda,x)=a_{\alpha}x^{\alpha}\uplambda^{\langle\alpha,z\rangle}$,
$\hat{L}_{\mathrm{lex}}(q)(\uplambda,x^{\prime})=a_{\alpha}x^{\prime\alpha}\uplambda^{\langle\alpha,z\rangle}$,
$\hat{L}_{\mathrm{lex}}(p)(\uplambda,x)=b_{\beta}x^{\beta}\uplambda^{\langle\beta,z\rangle}$
and
$\hat{L}_{\mathrm{lex}}(p)(\uplambda,x^{\prime})=b_{\beta}x^{\prime\beta}\uplambda^{\langle\beta,z\rangle}$.
It is obvious that there are two positive real numbers $\lambda_{1}$ and
$\lambda_{2}$ such that

$\hat{q}(\lambda,x),\hat{p}(\lambda,x),\hat{L}_{\mathrm{lex}}(q)(\lambda,x),\hat{L}_{\mathrm{lex}}(p)(\lambda,x)>0$

holds for all $\lambda\geq\lambda_{1}$ resp.

$\hat{q}(\lambda,x^{\prime}),\hat{p}(\lambda,x^{\prime}),\hat{L}_{\mathrm{lex}}(q)(\lambda,x^{\prime}),\hat{L}_{\mathrm{lex}}(p)(\lambda,x^{\prime})<0$

holds for all $\lambda\geq\lambda_{2}$.
Set
$M(\alpha)=\left\\{i:\mathrm{sgn}\left(x_{i}^{\alpha_{i}}\right)\neq\mathrm{sgn}\left(x_{i}^{\prime\alpha_{i}}\right)\right\\}$
and
$M(\beta)=\left\\{i:\mathrm{sgn}\left(x_{i}^{\beta_{i}}\right)\neq\mathrm{sgn}\left(x_{i}^{\prime\beta_{i}}\right)\right\\}$.
Then both sets are not empty, because otherwise we would get

$\mathrm{sgn}\left(\hat{L}_{\mathrm{lex}}(q)(\lambda,x)\right)=\mathrm{sgn}\left(\hat{L}_{\mathrm{lex}}(q)(\lambda,x^{\prime})\right)$

or

$\mathrm{sgn}\left(\hat{L}_{\mathrm{lex}}(p)(\lambda,x)\right)=\mathrm{sgn}\left(\hat{L}_{\mathrm{lex}}(p)(\lambda,x^{\prime})\right)$

for all $\lambda\geq 1$, which would result in a contradiction. The
intersection $M(\alpha)\cap M(\beta)$ is not empty, since $S(q)\subseteq
S(p)$:

Suppose that the intersection would be empty. Take an element $i\in M(\beta)$.
Set $y_{1}=(1,\ldots,1,x_{i},1,\ldots,1)$ and
$y_{2}=(1,\ldots,1,x_{i}^{\prime},1,\ldots,1)$. Then
$\hat{L}_{\mathrm{lex}}(p)(\lambda,y_{1})$ and
$\hat{L}_{\mathrm{lex}}(p)(\lambda,y_{2})$ have different signs for all
$\lambda\geq 1$, while $\hat{L}_{\mathrm{lex}}(q)(\lambda,y_{1})$ and
$\hat{L}_{\mathrm{lex}}(q)(\lambda,y_{2})$ have the same sign for all
$\lambda\geq 1$. Thus $\hat{q}(\lambda,y_{1})$ and $\hat{q}(\lambda,y_{2})$
are both negative or positive, while $\hat{p}(\lambda,y_{1})$ and
$\hat{p}(\lambda,y_{2})$ have different signs for all $\lambda$ large enough.
This can only work if $\hat{q}(\lambda,y_{1})$ and $\hat{q}(\lambda,y_{2})$
are negative for all sufficiently large $\lambda$. Thus $a_{\alpha}<0$. On the
other side, we get $b_{\beta}>0$ by repeating the same arguments with $j\in
M(\alpha)$. It is not hard to see that this cannot work. Set

$\tilde{y}=\begin{cases}(1,\ldots,1,x_{j}^{\prime},1,\ldots,1,x_{i}^{\prime},1,\ldots,1)\,\text{if}\,x_{j}^{\prime}<0\,\text{and}\,x_{i}^{\prime}<0\\\
(1,\ldots,1,x_{j},1,\ldots,1,x_{i}^{\prime},1,\ldots,1)\,\text{if}\,x_{j}<0\,\text{and}\,x_{i}^{\prime}<0\\\
(1,\ldots,1,x_{j},1,\ldots,1,x_{i},1,\ldots,1)\,\text{if}\,x_{j}<0\,\text{and}\,x_{i}<0\\\
(1,\ldots,1,x_{j}^{\prime},1,\ldots,1,x_{i},1,\ldots,1)\,\text{if}\,x_{j}^{\prime}<0\,\text{and}\,x_{i}<0\end{cases}$

In fact, $\hat{q}(\lambda,\tilde{y})$ is positive for all $\lambda>1$ large
enough, while $\hat{p}(\lambda,\tilde{y})$ is negative, since
$a_{\alpha}\tilde{y}^{\alpha}>0$ and $b_{\beta}\tilde{y}^{\beta}<0$. But this
contradicts $S(q)\subseteq S(p)$. Take a natural number $k$ out of the set
$M(\alpha)\cap M(\beta)$ and let $\pi_{k}$ denote the map
$\mathbb{R}^{n}\rightarrow\mathbb{R}^{n},(x_{1},\ldots,x_{k},\ldots,x_{n})\mapsto(x_{1},\ldots,-x_{k},\ldots,x_{n})$.
Then we can find a two positive real numbers $c_{1}$ and $c_{2}$ such that the
following two equations
$c_{1}\hat{L}_{\mathrm{lex}}(q)(\uplambda,\pi_{k}(x^{\prime}))=\hat{L}_{\mathrm{lex}}(q)(\uplambda,x)$
and
$c_{2}\hat{L}_{\mathrm{lex}}(p)(\uplambda,\pi_{k}(x^{\prime}))=\hat{L}_{\mathrm{lex}}(p)(\uplambda,x)$
hold. Hence it is easy to see that there is a positive real number
$\lambda_{1}^{\prime}$ such that

$c_{1}\hat{q}(\lambda,\pi_{k}(x^{\prime})),c_{2}\hat{p}(\lambda,\pi_{k}(x^{\prime}))>0$

for all $\lambda\geq\lambda_{1}^{\prime}$. On the other side, we still have

$c_{1}\hat{q}(\lambda,x^{\prime}),c_{2}\hat{p}(\lambda,x^{\prime})<0$

for all $\lambda\geq\lambda_{2}$. We are now ready to construct the rational
tentacles needed for the second part of Theorem 1.6. For each $i=1,\ldots,n$
we define
$\varphi_{i}=(\uplambda+\uplambda_{3}-1)^{z_{i}}\in\mathbb{R}[\uplambda]$,
where $\lambda_{3}=\max\\{\lambda_{1}^{\prime},\lambda_{2}\\}$. As usual set
$\varphi=(\varphi_{1},\ldots,\varphi_{n})$. By what we have done so far and by
using the fact that, if we substitute $x$ by some nearby other point $y$, none
of those inequalities used in this proof will be affected (see Proposition
2.3), we see that there is an open neighborhood $U$ of $x$ such that
$T_{\overline{U},\varphi}\in\mathcal{T}(S(q))$ and
$\pi_{k}(T_{\overline{U},\varphi})\notin\mathcal{T}(S(p))$. Repeating part
(ii)-(iv) of Theorem 1.6 concludes the proof. $\boxempty$

######  Remark 3.3.

Let $n$ be a natural number greater than $2$. Suppose $f$ and $g$ are two
different irreducible homogeneous polynomials in
$\mathbb{R}[\mathrm{x}_{1},\ldots,\mathrm{x}_{n}]$ of odd degree. If $f$ and
$g$ satisfy the condition (a) of Theorem 3.1, then it has some interesting
geometric consequences for $V_{1}=\mathcal{V}(f)$ and $V_{2}=\mathcal{V}(g)$.
Let $\Lambda_{1}$ resp. $\Lambda_{2}$ denote all singular points of
$V_{1}(\mathbb{R})$ resp. $V_{2}(\mathbb{R})$. According to [Sha I, Theorem 1,
p. 239] the intersection $V_{1}(\mathbb{R})\cap V_{2}(\mathbb{R})$ is not
empty. Consider a point $x\in V_{1}(\mathbb{R})\cap V_{2}(\mathbb{R})$. The
following cases may occur:

  * •

The point $x$ is in $\Lambda_{1}$: If $x$ is a local minimum of $g$, then is
must also a local minimum of $f$. Thus $x$ is a point in $\Lambda_{2}$.

  * •

The point $x$ is not in $\Lambda_{1}$: It is easy to see that $x$ is a
boundary point of $S(g)$. If $x$ is in $\Lambda_{2}$, then $x$ is either a
local minimum of $f$ or it is a saddle point of $f$. If $x$ is not in
$\Lambda_{2}$, then $V_{1}$ and $V_{2}$ intersect non-transversely at $x$.

[3, r, ,]

######  Example 3.4.

Consider the polynomials
$g=\mathrm{x}_{1}+\mathrm{x}_{2}+\mathrm{x}_{1}\mathrm{x}_{2}^{3}$ and
$f=\mathrm{x}_{1}^{5}+\mathrm{x}_{1}^{5}\mathrm{x}_{2}+\mathrm{x}_{2}^{2}$.
Let us check if $f$ and $g$ satisfy all conditions of Theorem 3.1.
By using [Mathematica] as we did in 4.1, we see that $S(g)\subseteq S(f)$.
(a): Obviously, $L_{\mathrm{lex}}(g)=\mathrm{x}_{1}\mathrm{x}_{2}^{3}$ and
$L_{\mathrm{lex}}(f)=\mathrm{x}_{1}^{5}\mathrm{x}_{2}$. Since
$\deg_{\mathrm{x}_{2}}\left(L_{\mathrm{lex}}(g)\right)>\deg_{\mathrm{x}_{2}}\left(L_{\mathrm{lex}}(f)\right)$,
the polynomial $L_{\mathrm{lex}}(f)$ cannot lie in
$L_{\mathrm{lex}}(g)\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2}]$.
(b): Set $y^{\prime}=(5,5)$, $y=(5,-5)$, and consider
$\hat{g}(\uplambda,y^{\prime})=g\left(\uplambda^{z_{1}}y_{1}^{\prime},\uplambda^{z_{2}}y_{2}^{\prime}\right)=5\uplambda^{z_{1}}+5\uplambda^{z_{2}}+625\uplambda^{z_{1}+3z_{2}}$,
$\hat{f}(\uplambda,y^{\prime})=f\left(\uplambda^{z_{1}}y_{1}^{\prime},\uplambda^{z_{2}}y_{2}^{\prime}\right)=3125\uplambda^{5z_{1}}+25\uplambda^{2z_{2}}+15625\uplambda^{5z_{1}+z_{2}}$,
$\hat{g}(\uplambda,y)=g\left(\uplambda^{z_{1}}y_{1},\uplambda^{z_{2}}y_{2}\right)=5\uplambda^{z_{1}}-5\uplambda^{z_{2}}-625\uplambda^{z_{1}+3z_{2}}$,
$\hat{f}(\uplambda,y)=f\left(\uplambda^{z_{1}}y_{1},\uplambda^{z_{2}}y_{2}\right)=3125\uplambda^{5z_{1}}+25\uplambda^{2z_{2}}-15625\uplambda^{5z_{1}+z_{2}}$,
where $(z_{1},z_{2})\in\mathbb{Z}^{2}$. It is easy to see that if we take
$(z_{1},z_{2})\in\mathbb{N}_{1}$, then
$\hat{g}(\lambda,y^{\prime}),\hat{f}(\lambda,y^{\prime})$ are positive for all
$\lambda\geq 1$, while
$\hat{g}(\lambda,y^{\prime}),\hat{f}(\lambda,y^{\prime})$ are negative for all
$\lambda\geq 1$. By taking a small compact neighborhood $U^{\prime}$ of
$y^{\prime}$ resp. $U$ of $y$, we get a tentacle $T_{U^{\prime},z}$ belonging
to $\mathcal{T}(S(g))$ and another one $T_{U,z}$ belonging to
$\mathcal{T}(S(-f))$. Finally, Theorem 3.1 implies that there is no non-
negative polynomial $t\in\mathbb{R}[x_{1},x_{2}]$ such that $f(x)-t(x)g(x)\geq
0$ for all $x\in\mathbb{R}^{2}$.

[8, r, ,]

######  Example 3.5.

Let us revisit the counterexample
$f=\mathrm{x}_{1}^{3}+\mathrm{x}_{1}^{3}\mathrm{x}_{2}+\mathrm{x}_{2}^{2}$ and
$g=\mathrm{x}_{1}+\mathrm{x}_{2}+\mathrm{x}_{1}\mathrm{x}_{2}$. Unfortunately,
condition (a) of Theorem 3.1 is violated. We have
$L_{\mathrm{lex}}(g)=\mathrm{x}_{1}\mathrm{x}_{2}$ and
$L_{\mathrm{lex}}(f)=\mathrm{x}_{1}^{3}\mathrm{x}_{2}$. Thus
$L_{\mathrm{lex}}(f)=\mathrm{x}_{1}^{2}L_{\mathrm{lex}}(g)$. On the other
side, it is easy to see that the condition (b) of Theorem 3.1 is satisfied.
Can we still make use of Theorem 3.1?
Let $t\in\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2}]$ be a non-negative
polynomial. If $L_{\mathrm{lex}}(t)\neq\mathrm{x}_{1}^{2}$ then we can just
repeat the arguments of Theorem 3.1. However, the problematic case is if
$L_{\mathrm{lex}}(t)=\mathrm{x}_{1}^{2}$ holds. Interchange the variables
$\mathrm{x}_{1}$ and $\mathrm{x}_{2}$ in $f$ and $g$, giving
$f=\mathrm{x}_{2}^{3}+\mathrm{x}_{2}^{3}\mathrm{x}_{1}+\mathrm{x}_{1}^{2}$ and
$g=\mathrm{x}_{1}+\mathrm{x}_{2}+\mathrm{x}_{1}\mathrm{x}_{2}$. Now, we are
out for some suitable $z$-grading such that $L_{z}(f)\notin
L_{z}(g)\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2}]$. Choose
$z\in\mathbb{N}_{1}^{2}$ with $3z_{2}=z_{1}$. Then we have
$L_{z}(g)=\mathrm{x}_{1}\mathrm{x}_{2}$,
$L_{z}(f)=\mathrm{x}_{2}^{3}\mathrm{x}_{1}+\mathrm{x}_{1}^{2}$ resp.
$\hat{L}_{z}(g)(\uplambda,\mathrm{x})=\mathrm{x}_{1}\mathrm{x}_{2}\uplambda^{z_{1}+z_{2}}$,
$\hat{L}_{z}(f)(\uplambda,\mathrm{x})=\left(\mathrm{x}_{2}^{3}\mathrm{x}_{1}+\mathrm{x}_{1}^{2}\right)\uplambda^{2z_{1}}$.
Obviously, $L_{z}(f)$ does not lie in
$L_{z}(g)\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2}]$. Let us check condition
(b). Set $x=(5,5)$. Then
$\hat{L}_{z}(g)(\uplambda,x)=25\uplambda^{z_{1}+z_{2}}$ and
$\hat{L}_{z}(f)(\uplambda,x)=5^{4}\uplambda^{2z_{1}}+5^{2}\uplambda^{2z_{1}}$.
On the other hand, if we set $x=(-5,5)$ we get
$\hat{L}_{z}(g)(\uplambda,x)=-25\uplambda^{z_{1}+z_{2}}$ and
$\hat{L}_{z}(f)(\uplambda,x)=-5^{4}\uplambda^{2z_{1}}+5^{2}\uplambda^{2z_{1}}$.
This means that if we take an appropriate neighborhood $U$ of $x$ and define
$\pi:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2},(x_{1},x_{2})\mapsto(-x_{1},x_{2})$
we get $T_{U,z}\in\mathcal{T}(S(g))$ and $\pi(T_{U,z})\in\mathcal{T}(S(-f))$.
By using the here defined $z$-grading, instead of the one defined in the proof
of Theorem 3.1, and repeating the arguments in part (ii)-(iv) we get that
there is no non-negative polynomial
$t\in\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2}]$ such that $f(y)-t(y)g(y)\geq
0$ for all $y\in\mathbb{R}^{2}$. In fact we can state:

######  Theorem 3.6.

Let $S_{1}=S(q)$ and $S_{2}=S(p)$ be two semi-algebraic sets in
$\mathbb{R}^{n}$. Suppose the following condition is satisfied:

There is a $z\in\mathbb{N}_{1}^{n}$ such that $L_{z}(p)\notin L_{z}(q)A$, a
rational tentacle $T\in\mathcal{T}(S_{1})$ of degree $z$ and an element
$\pi\in G$ such that $\pi(T)\notin\mathcal{T}(S_{2})$. Furthermore, all
unbounded $T^{\prime}\in\mathcal{T}(S_{1})$ with
$\pi(T^{\prime})\subseteq\pi(T)$ satisfy
$\pi(T^{\prime})\notin\mathcal{T}(S_{2})$.

Then there is no non-negative polynomial $t\in A$ such that $p(y)-t(y)q(y)\geq
0$ for all $y\in\mathbb{R}^{n}$.

Proof: Simply repeat the same arguments in part (ii)-(iv) of Theorem 3.1.
$\boxempty$

In the same manner we can modify Proposition 3.2:

######  Proposition 3.7.

Let $S_{1}=S(q)$ and $S_{2}=S(p)$ be two semi-algebraic sets in
$\mathbb{R}^{n}$, where neither $q=0$ nor $p=0$. Suppose the following
condition is satisfied:

There is a $z\in\mathbb{N}_{1}^{n}$ such that $L_{z}(p)\notin L_{z}(q)A$ and
the quadratic modules $\mathrm{QM}(q)$ and $\mathrm{QM}(-p)$ are totally
stable with respect to the $z$-grading.

Then there is no non-negative polynomial $t\in A$ such that $p(y)-t(y)q(y)\geq
0$ for all $y\in\mathbb{R}^{n}$.

######  Remark 3.8.

The difference between Theorem 3.1 and Theorem 3.6 is simple. In Theorem 3.1
we demanded that $L_{\mathrm{lex}}(p)\notin L_{\mathrm{lex}}(q)A$. Then we
constructed a special $z$-grading, where $z\in\mathbb{N}_{1}^{n}$. Under this
$z$-grading we had $L_{z}(p)=L_{\mathrm{lex}}(p)$ resp.
$L_{z}(q)=L_{\mathrm{lex}}(q)$ and therefore $L_{z}(p)\notin L_{z}(q)A$. Let
us refer to the $z$-gradings satisfy $z\in\mathbb{N}_{1}^{n}$ and
$L_{z}(p)\notin L_{z}(q)A$ as special gradings. The difference between Theorem
3.6 and Theorem 3.1 is, that the special $z$-grading constructed in part (i)
of the proof of Theorem 3.1, is already given in the prerequisites of Theorem
3.6. The disadvantage of Theorem 3.6 compared to Theorem 3.1 is, that one must
find such a special $z$-grading for Theorem 3.6 to work, while Theorem 3.1
does not require such a procedure. The advantage of Theorem 3.6 is, that it
allows a wider range of gradings as Example 3.5 illustrates it. In fact, we
have an two explanations why we used
$g=\mathrm{x}_{1}\mathrm{x}_{3}+\mathrm{x}_{2}\mathrm{x}_{3}+\mathrm{x}_{1}\mathrm{x}_{2}$.
The first explanation is a geometric one. According to Example 3.5 the
polynomial $g$ has all the necessary properties for Theorem 3.6 to work. The
second explanation is an algebraic one. According to Theorem 2.5 quadratic
forms $q$ that have a vanishing diagonal part give rise to quadratic module
$\mathrm{QM}(q)$ that is totally stably with respect to any $z$-grading. Thus
it is (was) convenient to choose
$g=\mathrm{x}_{1}\mathrm{x}_{3}+\mathrm{x}_{2}\mathrm{x}_{3}+\mathrm{x}_{1}\mathrm{x}_{2}$.

###  4 A non-geometric counterexample

[2, r, ,] In the last chapter we saw, that Counterexample 6 proved the
S4-conjecture wrong because of geometric reasons. The straightforward question
is obvious: Can we find for all counterexamples a geometric reason? At
least,we will give a counterexample that does not work because of a
arithmetical reason. From now on, we set
$g=(-3+\mathrm{x}_{1}-\mathrm{x}_{2})(3+\mathrm{x}_{1}-\mathrm{x}_{2})$,
$p=-\mathrm{x}_{1}^{3}+\mathrm{x}_{2}^{3}+2\mathrm{x}_{1}+1$,
$l_{1}=-3+\mathrm{x}_{1}-\mathrm{x}_{2}$,
$l_{2}=3+\mathrm{x}_{1}-\mathrm{x}_{2}$ and $f=-l_{2}p$. By using
[Mathematica] (see 6) it is easy to verify that $S(g)\subseteq S(f)$. The next
obvious step is:

######  Proposition 4.1.

There is no non-negative polynomial
$t\in\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2}]_{2}$ such that
$f(y)-t(y)g(y)\geq 0$ for all $y\in\mathbb{R}^{2}$.

Proof: Suppose there is a non-negative polynomial
$t\in\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2}]_{2}$ contradicting the
statement of this proposition. In the following we fix a real number
$x_{2}\in\mathbb{R}$ and consider
$f(\mathrm{x}_{1},x_{2}),g(\mathrm{x}_{1},x_{2}),t(\mathrm{x}_{1},x_{2})$ as
polynomials in $\mathrm{x}_{1}$. The polynomial
$f(\mathrm{x}_{1},x_{2})-t(\mathrm{x}_{1},x_{2})g(\mathrm{x}_{1},x_{2})\in\mathbb{R}[\mathrm{x}_{1}]$
has a root at $x_{1}=x_{2}-3$. Since
$f(\mathrm{x}_{1},x_{2})-t(\mathrm{x}_{1},x_{2})g(\mathrm{x}_{1},x_{2})\in\mathbb{R}[\mathrm{x}_{1}]$
is non-negative for every $x_{2}\in\mathbb{R}$, it must be divided by
$l_{2}^{2}(\mathrm{x}_{1},x_{2})\in\mathbb{R}[\mathrm{x}_{1}]$.

The remainder of the polynomial division
$f(\mathrm{x}_{1},x_{2}):l_{2}^{2}(\mathrm{x}_{1},x_{2})$, as polynomials in
$\mathbb{R}[\mathrm{x}_{1}]$, is
$r_{1}(\mathrm{x}_{1},x_{2})=9x_{2}^{3}-9\mathrm{x}_{1}x_{2}^{2}-52x_{2}^{2}-25\mathrm{x}_{1}x_{2}+97x_{2}^{2}-22\mathrm{x}_{1}-66\in\mathbb{R}[\mathrm{x}_{1}]$.
And for $g(\mathrm{x}_{1},x_{2}):l_{2}^{2}(\mathrm{x}_{1},x_{2})$ we have
$r_{2}(\mathrm{x}_{1},x_{2})=6x_{2}-6\mathrm{x}_{1}-18\in\mathbb{R}[\mathrm{x}_{1}]$.
Finally, let $r_{3}(\mathrm{x}_{1},x_{2})\in\mathbb{R}[\mathrm{x}_{1}]$ denote
the remainder of $t(\mathrm{x}_{1},x_{2}):l_{2}^{2}(\mathrm{x}_{1},x_{2})$.
Since $f(\mathrm{x}_{1},x_{2})-t(\mathrm{x}_{1},x_{2})g(\mathrm{x}_{1},x_{2})$
is divided by $l_{2}^{2}(\mathrm{x}_{1},x_{2})$, we get the identity
$r_{1}(\mathrm{x}_{1},x_{2})-r_{3}(\mathrm{x}_{1},x_{2})r_{2}(\mathrm{x}_{1},x_{2})=0$.
This leads to
$r_{3}(\mathrm{x}_{1},x_{2})=\frac{r_{1}(\mathrm{x}_{1},x_{2})}{r_{2}(\mathrm{x}_{1},x_{2})}=\frac{1}{6}\left(22-25x_{2}+9x_{2}^{2}\right)$.
Set
$\tilde{t}=a_{x_{2}}l_{2}^{2}+\frac{1}{6}\left(22-25\mathrm{x}_{2}+9\mathrm{x}_{2}^{2}\right)$
and choose $a_{x_{2}}\in\mathbb{R}$ such that the equality
$\tilde{t}(\mathrm{x}_{1},x_{2})=t(\mathrm{x}_{1},x_{2})$ holds for
$x_{2}\in\mathbb{R}$. It is easy to see that the leading term of
$f-\tilde{t}g$ in $\mathrm{x}_{2}$ is
$\left(-\frac{1}{2}-a_{x_{2}}\right)\mathrm{x}_{2}^{4}$ and that
$f(0,0)-\tilde{t}(0,0)g(0,0)=30+81a_{x_{2}}$. For large $x_{2}\in\mathbb{R}$
we see that $a_{x_{2}}$ must satisfy $a_{x_{2}}\leq-\frac{1}{2}$ and
$a_{x_{2}}\geq-\frac{30}{81}$, which is impossible. $\boxempty$

######  Remark 4.2.

Under an appropriate change of coordinates the homogenization of $g$ can be
written as $\overline{g}=-9\mathrm{x}_{1}^{2}+2\mathrm{x}_{2}^{2}$. Thus the
signature of $\overline{g}$ is $0$. Note, that the signature in Counterexample
4.1 was $-1$. Furthermore, it is easy to see that neither Theorem 3.1 nor
Theorem 3.6 can be applied. In fact, $g$ violates condition (b) in Theorem 3.1
and the condition in Theorem 3.6.

###  5 Final thoughts

Finally, we have come to the end of this article. Hence let us summarize what
we have learned so far. First, we learned that the S4-conjecture is not true.
Second, we learned that there are geometric reasons why the S4-conjecture
cannot be true. Finally, we learned that there are arithmetic reasons why the
S4-conjecture cannot work. Still there are many questions left. The signature
of $g$ in 4.1 resp. 4.1 was $-1$ resp. $0$. So it is quite naturally to ask,
if there is an counterexample, where $g$ has signature $1$. There is another
obvious question: Under which conditions does the S4-conjecture work? Can
these conditions be expressed in geometric or algebraic terms? Results in this
direction can be found in [Sch, Corollary 4.5]:

######  Proposition 5.1.

[Sch, Corollary 4.5]: Let
$h_{1},\ldots,h_{r}\in\mathbb{R}[\mathrm{x}_{0},\ldots,\mathrm{x}_{n}]$ be
homogeneous polynomials of even degree, and let

$S=S(h_{1},\ldots,h_{r}).$

Assume there is $\xi\in\mathbb{R}^{n+1}$ with $h_{i}(\xi)>0$ for
$i=1,\ldots,r$. If $p,q\in\mathbb{R}[\mathrm{x}_{0},\ldots,\mathrm{x}_{n}]$
are homogeneous of even degree and positive on $S\backslash\\{0\\}$, then
$pq^{m}$ lies in the preordering generated by $h_{1},\ldots,h_{r}$, for all
sufficiently large $m\geq 0$.

In contrast to the S4-conjecture, we need a homogeneous polynomial $p$ of even
degree that is positive on the set $S\backslash\\{0\\}$. And even then we can
only conclude that there is a natural number $m\geq 1$ such that $p^{m}\in
T(h_{1},\ldots,h_{r})$. If we could show that $p\in T(h_{1},\ldots,h_{r})$, we
would still have considerable obstacles. For example, we do not know what kind
of degree bounds the various representations of $p$ in $T(h_{1},\ldots,h_{r})$
have. Nevertheless, let us consider the polynomials
$g=\mathrm{x}_{1}^{2}+\mathrm{x}_{2}^{2}-4\mathrm{x}_{3}^{2}$ and
$f=\mathrm{x}_{1}^{4}+\mathrm{x}_{2}^{4}-\mathrm{x}_{3}^{4}$. Then $f$ and $g$
satisfy the following conditions:

  * •

There is a point $x^{\prime}\in\mathbb{R}^{3}$ such that $g(x^{\prime})>0$.

  * •

The polynomial $f$ is positive on the set $S(g)\backslash\\{0\\}$.

  * •

The polynomial $f-\mathrm{x}_{3}^{2}g$ is non-negative.

In other words, the two polynomials satisfy the S4-conjecture. It is easy to
see that $V_{1}$ and $V_{2}$ are both non-singular curves. According to [BC,
Proposition 11.6.2, p. 286] the set
$V_{1}(\mathbb{R})\subseteq\mathbb{P}^{2}(\mathbb{R})$ decomposes into 1 or 0
ovals, and the set $V_{2}(\mathbb{R})\subseteq\mathbb{P}^{2}(\mathbb{R})$ into
at most $4$ different ovals. For the definition of an oval see [BC, p. 286].
In our case $V_{1}(\mathbb{R})$ and $V_{2}(\mathbb{R})$ decompose into one
oval that does not intersect the plane at infinity. Since
$\mathrm{x}_{1}^{2}+\mathrm{x}_{2}^{2}$ and
$\mathrm{x}_{1}^{4}+\mathrm{x}_{2}^{4}$ are positive on the set
$\mathbb{R}^{2}\backslash\\{0\\}$, the two sets $V_{1}(\mathbb{R})$,
$V_{2}(\mathbb{R})$ do not intersect the hyperplane at infinity. Another
consequence is that we can tell something about the geometry of $S(g)$ resp.
$S(f)$. For $y\in\mathbb{R}\backslash\\{0\\}$ set
$H_{y}=\mathbb{R}^{2}\times\\{y\\}$ and interpret $V_{1}$ and $V_{2}$ as
affine varieties in $\mathbb{A}^{3}$. Then we have $\partial S(g)\cap
H_{y}=V_{1}(\mathbb{R})\cap H_{y}$ and $\partial S(f)\cap
H_{y}=V_{2}(\mathbb{R})\cap H_{y}$ for all $y\in\mathbb{R}\backslash\\{0\\}$.
The two inclusions $\partial S(g)\cap H_{y}\subseteq V_{1}(\mathbb{R})\cap
H_{y}$ and $\partial S(f)\cap H_{y}\subseteq V_{2}(\mathbb{R})\cap H_{y}$ are
obvious. The other two inclusions $\partial S(g)\cap H_{y}\supseteq
V_{1}(\mathbb{R})\cap H_{y}$ and $\partial S(f)\cap H_{y}\supseteq
V_{2}(\mathbb{R})\cap H_{y}$ hold, because $V_{1}(\mathbb{R})$ and
$V_{2}(\mathbb{R})$ do not have singular points in $H_{y}$. Thus each slice
$\partial S(f)\cap H_{y}$, $\partial S(g)\cap H_{y}$ ’looks’ like a circle. If
we replace $g$ by an arbitrary quadratic form
$q\in\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2},\mathrm{x}_{3}]$ and $f$ by an
arbitrary ternary quartic $p$, then all geometric statements111Of course with
an adjusted number of ovals made so far in this matter remain true, if the
ovals of $\mathcal{V}(q)(\mathbb{R})$ and $\mathcal{V}(p)(\mathbb{R})$ do not
intersect the plane at infinity. Interestingly, the ovals that do not
intersect the plane at infinity have different topological properties than
their counterparts that intersect the plane at infinity: It is a well known
fact that the fundamental group of $\mathbb{P}^{2}(\mathbb{R})$ is exactly
$\mathbb{Z}/2\mathbb{Z}$. By interpreting an oval as a loop, it turns out that
all ovals that do not intersect the hyperplane at infinity, represent the
identity element of the fundamental group. If all ovals of
$\mathcal{V}(q)(\mathbb{R})$ and $\mathcal{V}(p)(\mathbb{R})$ do not intersect
the hyper plane at infinity, then it is obvious that the topological situation
compared to $\mathcal{V}(g)(\mathbb{R})$ and $\mathcal{V}(f)(\mathbb{R})$ has
not changed much.
Hence it is convenient to ask this final question:

######  Question 5.2.

Let $q\in\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2},\mathrm{x}_{3}]$ be a
quadratic form and $p$ a ternary quartic. Set
$V_{1}=\mathcal{V}(q)\subseteq\mathbb{P}^{2}$ and
$V_{2}=\mathcal{V}(p)\subseteq\mathbb{P}^{2}$. Suppose that the following
conditions are satisfied:

  * •

There is a point $x^{\prime}\in\mathbb{R}^{3}$ such that $q(x^{\prime})>0$.

  * •

The ternary quartic $p$ is positive on the set $S(q)\backslash\\{0\\}$.

  * •

The projective varieties $V_{1}$ and $V_{2}$ are non-singular.

  * •

The set $V_{1}(\mathbb{R})$ is an oval and $V_{2}(\mathbb{R})$ decomposes into
at least one oval.

  * •

All ovals of $V_{1}(\mathbb{R})$ and $V_{2}(\mathbb{R})$ do not intersect the
hyperplane at infinity.

Can we find a non-negative homogeneous polynomial
$t\in\mathbb{R}[\mathrm{x}_{1},\mathrm{x}_{2},\mathrm{x}_{3}]_{2}$ such that
$p(y)-t(y)q(y)\geq 0$ for all $y\in\mathbb{R}^{3}$?

## References

  * [Bar] Alexander Barvinok: A Course in Convexity, Graduate Studies in Mathematics, AMS Volume 54, ISSN 1065-7339.
  * [Sha I] Igor Shavarevich: Basic Algebraic Geometry I, Springer, Second Edition, 1994.
  * [Sha II] Igor Shavarevich: Basic Algebraic Geometry II, Springer, Second Edition, 1997.
  * [BC] Jacek Bochnak, Michel Coste, Marie-Francois Roy: Real Algebraic Geometry, Springer, 1991.
  * [PT] Imre Polik, Tamas Terlaky: A Survey of the S-lemma, SIAM journals, http://epubs.siam.org/doi/abs/10.1137/S003614450444614X?journalCode=siread.
  * [ZS] Kuize Zhang, Lijun Zhang, Fuchun Sun: High-Order S-Lemma With Application To Stability Of A Class Of Switched Nonlinear Systems, http://arxiv.org/abs/1403.1016.
  * [Ne] Tim Netzer, Stability of Quadratic Modules, http://arxiv.org/abs/0807.4403.
  * [RS] Remmert, Schumacher, Funktionentheorie 1, Springer, 5.Auflage, 2001.
  * [Hi] D.Hilbert: Über die Darstellung definiter Fromen als Summe von Formenquadraten, Math. Ann. 32 342-350, 1888.
  * [Mathematica] Wolfram Research, Inc., Mathematica, Version 10.0, Champaign, IL (2014)
  * [MathOver] Mathoverflow, http://mathoverflow.net/questions/112230/an-s-lemma-for-polynomials-of-degree-4-in-three-variables
  * [Sch] C.Scheiderer: A Positivstellensatz for projective real varieties, http://arxiv.org/abs/1104.1772
  * [Bo] S.Bosch: Algebra, 7.Auflage, Springer, 2008


11institutetext: Chinese University of Hong Kong

# Cascaded Robust Learning at Imperfect Labels for Chest X-ray Segmentation

Cheng Xue  Qiao Deng  Qi Dou  Pheng-Ann Heng

###### Abstract

The superior performance of CNN on medical image analysis heavily depends on
the annotation quality, such as the number of labeled image, the source of
image, and the expert experience. The annotation requires great expertise and
labour. To deal with the high inter-rater variability, the study of imperfect
label has great significance in medical image segmentation tasks. In this
paper, we present a novel cascaded robust learning framework for chest X-ray
segmentation with imperfect annotation. Our model consists of three
independent network, which can effectively learn useful information from the
peer networks. The framework includes two stages. In the first stage, we
select the clean annotated samples via a model committee setting, the networks
are trained by minimizing a segmentation loss using the selected clean
samples. In the second stage, we design a joint optimization framework with
label correction to gradually correct the wrong annotation and improve the
network performance. We conduct experiments on the public chest X-ray image
datasets collected by Shenzhen Hospital. The results show that our methods
could achieve a significant improvement on the accuracy in segmentation tasks
compared to the previous methods.

###### Keywords:

Robust learning imperfect label Segmentation.

##  1 Introduction

Deep neural networks (DNNs) have achieved human-level performance on many
medical image analysis tasks, such as melanoma diagnosis [5], pulmonary
nodules detection [13], retinal disease [4], and lumpy node metastases
detection [1]. These outstanding performances heavily rely on massive training
data with high-quality annotations. Annotation of medical images, especially
for pixel-level annotation for segmentation tasks, is costly and time-
consuming. The process is experience-prone, while the annotations from
different clinical experts may have disagreements that are usually inevitable
for the blurred boundary of lesions and organs.

Previous studies show that the DNNs trained by noisy labeled datasets can
cause performance degradation. That is because the huge memory capacity and
strong learning ability of DNNs can remember the noisy labels and easily
overfit to them [17, 18, 14]. Hence, it is important to develop DNNs with
strong robustness to noisy labels. Many studies have addressed the issue of
noisy label in medical analysis community. Goldberger et al.[6] added an
additional softmax layer to estimate the correct labels. Xue et al. [16]
proposed to consider the noisy sample and hard sample by an on-line sample
selection module and re-weighting module. Zhu et al. [18] proposed the
automatic quality evaluation module and overfitting control module to update
the network parameters. Shu et al.[14] presented a LVC-Net losses function by
combining noisy labels with image local visual cues to generate better
semantic segmentation. Le et al. [10] utilized a small set of clean training
samples to assign weights to training samples.

To tackle the challenging problem of noisy labeled segmentation masks, we
present a cascaded learning framework for lung segmentation using the X-ray
images with imperfectly annotated ground truth. In the first stage, our
framework selects clean annotated samples according to the prediction
confidence and uncertainty of samples, that is inspired by the ideas of Co-
teaching [7]. Specifically, our model consists of three independent networks
being trained simultaneously, each network is real-time updated according to
the prediction results of the other two networks. For a clean annotated
sample, the three networks tend to produce high confidence prediction with
smaller inter-rater variance. Thus, the samples with close prediction and high
confidence are selected as the high-quality sample, which will be only used to
contribute in the weight backpropagation process. Since the selection stage
leads to a low utilization efficiency of the valuable training data, we
propose a label correction module in the second stage, which can correct the
imperfect label. Furthermore, a joint optimization scheme is designed to
cooperatively supervise the three networks with the original label and the
corrected one. Our method was extensively evaluated on the dataset Shenzhen
chest x-ray [8, 3, 15]. The results demonstrate a good capability of our
method to the issue of noisy label, that the cascaded robust learning
framework can more accurately perform the lung segmentation comparing to other
methods.

##  2 Method

Figure 1: Illustration of the pipeline of our cascaded robust learning
framework. (a) shows the first sample selection stage, where three networks
trained independently, but updated according to the prediction of the other
two peer networks. (b) and (c) are the second stage. (b) shows our proposed
label correction module, using the average prediction of two peer networks
followed by a sharpening function to produce corrected label $\bar{y}$. (c)
shows the joint optimization scheme, the network is supervised by the original
label $\hat{y}$ and the corrected label $\bar{y}$. The final output is given
by the average of the three networks.

Fig 1. illustrates the framework of our cascaded robust learning method. In
the first stage the sample selection module filters the clean samples and
update the three networks with the selected clean samples. In the second
stage, our method start to correct the imperfect labels, then use the
corrected label and original label to jointly optimize the three networks.

###  2.1 Sample Selection Stage

We study the task of chest x-ray segmentation, where the training set contains
images $x$ and noisy labeled ground truth $\hat{y}$, while the clean ground
truth $y$ is unknown. The goal for this fully supervised segmentation task is
to minimize the following object function:

$\min_{\theta}\sum_{i=1}^{N}\mathcal{L}(f(x_{i};\theta),\hat{y}_{i})$ (1)

where $\mathcal{L}$ denotes the loss function (e.g., cross-entropy loss) to
evaluate the quality of the network output on inputs. $f(\theta)$ denotes the
segmentation neural network with weights $\theta$.

Recent studies show that by updating the network with high confidence samples
can improve the robustness to noisy labels [11, 9, 7]. Therefore, we propose a
novel sample selection framework (SS) to select high confidence samples as the
useful training instances. Our framework consisted by three independent
networks, where they have identical architecture. We adopt the vanilla
U-Net[12] as the classifier in our experiment. In the training process, we
select the high confidence samples with small uncertainty to update each
network, because those samples are more likely to be clean labeled instances.
In our experiment, we select half batch data as useful information.
Concretely, the three networks feed forward and predict the same mini-batch of
data. Then for each network, the useful samples for weight updating is
obtained by the other two networks as shown in Fig. 1(a). Taking network A as
an example, the useful sample for network A is obtained from network B and C,
where we first filter out the high uncertainty samples by exclude the ones
showing disagreed prediction, then among the low uncertainty samples, the
small loss samples was further selected as useful samples for network A. We
calculate the uncertainty according to equation 2.

$\mu=|\mathcal{L}(f_{B}(x_{i};\theta_{B}),\hat{y}_{i})-\mathcal{L}(f_{C}(x_{i};\theta_{C}),\hat{y}_{i})|$
(2)

where $\mathcal{L}$ denotes the cross-entropy loss. $f_{2}$ and $f_{3}$ denote
the network B and network C. $\theta_{B}$ and $\theta_{C}$ represent the
weight of network B and C.

###  2.2 Joint Optimization with Label Correction

In the stage of sample selection, only partial samples can be used for
training, where it does not take fully advantage of the imperfect training
data. Therefore, we design a joint optimization (JO) framework to train the
network with the original label and corrected label, so that the utilization
efficiency of training data can be maintained. In order to correct the noisy
label, we design an label correction module to work together with the joint
optimization scheme.

#### Label correction

The sample selection stage first trains an initial network by using image x
with noisy label $\hat{y}$. Then we proceed to the label correction phase, as
shown in Figure 1 (b). We compute the average of three model’s prediction in
each iteration, that is followed by an entropy minimization step widely
adopted in semi-supervised learning [2]. Specifically, for the average
prediction of the three models, we apply a sharpening function to reduce the
entropy of the per pixel label distribution through adjusting the temperature.
The sharpening function is the equation 3.

$\displaystyle\begin{split}q=\frac{1}{2}(f_{B}((x_{i};\theta_{B}),\hat{y}_{i})+f_{C}((x_{i};\theta_{C}),\hat{y}_{i}))\\\
sharpen(q,T)_{i}=q_{i}^{\frac{1}{T}}/\sum_{j=1}^{L}q_{j}^{\frac{1}{T}}\end{split}$
(3)

where $q$ is the average prediction feature map over two models, $T$ is a
hyperparameter that adjusts the temperature. As $T$ closes to zero, the output
of $Sharpen(q,T)$ will approach a one-hot distribution. Since we will use
$q=Sharpen(q,T)$ as a corrected target for the model’s prediction later,
$T=0.5$ is chosen to encourage the model to produce lower-entropy prediction.

#### Joint optimization

We start the joint optimization stage after $k$ epochs of sample selection.
For each uncertain sample, we produce a corrected label for the imperfect
input by the label correction module. The corrected label is used in the
training process together with the original label as a complementary
supervision to jointly supervise the network, as shown in Equation 4.

$\mathcal{L}_{total}=\alpha\times\mathcal{L}(f(x_{i};\theta),\hat{y}_{i})+(1-\alpha)\times\mathcal{L}(f(x_{i};\theta),\bar{y}_{i})$
(4)

where $\mathcal{L}$ is the cross entropy loss, $\hat{y}$ is the original noisy
label, and $\bar{y}$ is the corrected label produced by the label correction
phase. The weight factor $\alpha$ controls the important weight of the two
terms, we set $\alpha=0.5$ in our study.

##  3 Experiments

###  3.1 Dataset and Pre-processing.

We evaluated our method on the public Shenzhen chest x-ray datasets [8, 3,
15], the segmentation mask were prepared manually by Computer Engineering
Department, Faculty of Informatics and Computer Engineering, National
Technical University of Ukraine. The dataset contains 566 chest x-ray images
and each image has the left and the right lungs. We split the 566 chest x-ray
images into 396 images for training and 170 for evaluation. All the images
were resized to $256\times 256$, and normalized as zero mean and unit
variance.

###  3.2  Implementation

The framework was implemented in PyTorch, using a TITAN Xp GPU. We used the
SGD optimizer to update the network parameters with weight decay of 0.001 and
momentum of 0.9. We adopt an exponential learning rate with initial learning
rate set as 0.001. We totally trained 100 epochs, the batch size was 32. We
adopted the data augmentation including randomly rotation and randomly
horizontal flipping. In order to produce noisy label for the training data, we
randomly selected 25$\%$, 50$\%$, 75$\%$ samples from the training set to
erode or dilate with the number of iterations $n$ between $5-15$ ($5\leq n\leq
15$). We adopted the dice coefficient as evaluation criteria for segmentation
accuracy evaluation. Fig.2 shows the example of some noisy annotation of the
segmentation mask.

|  |  |
---|---|---|---
(a) chest X-ray | (b) ground truth | (c)dilation | (d)erosion
Figure 2: Examples of the noise annotation  |  |  |
---|---|---|---
|  |  |
(a) Chest X-ray | (b) ground truth | (c)noise 75 $\%$ baseline | (d)noise 75 $\%$ ours
Figure 3: Two examples of the segmentation results in the test data by
different methods. (a) is the input image, (b) is the ground truth. (c) and
(d) shows the results of U-Net and our method under 75$\%$ noise ratio.

###  3.3 Quantitative Evaluation

Table 1: Comparison between our method and various methods. Noise ratio | Noise level | Strategy | Dice | k
---|---|---|---|---
|  |  |  | 20 | 50 | 80
No noise | - | Vanilla U-Net | 0.8989 | - | - | -
No noise | - | Co-teaching [7] | 0.9146 | - | - | -
No noise | - | SS + JO | - | 0.9252 | 0.9250 | 0.9236
25$\%$ | $5\leq n\leq 15$ | Vanilla U-Net | 0.8758 | - | - | -
25$\%$ | $5\leq n\leq 15$ | Co-teaching [7] | 0.8906 | - | - | -
25$\%$ | $5\leq n\leq 15$ | SS | 0.9142 | - | - | -
25$\%$ | $5\leq n\leq 15$ | SS + JO | - | 0.9214 | 0.9253 | 0.9248
50$\%$ | $5\leq n\leq 15$ | Vanilla U-Net | 0.8478 | - | - | -
50$\%$ | $5\leq n\leq 15$ | Co-teaching [7] | 0.8896 | - | - | -
50$\%$ | $5\leq n\leq 15$ | SS | 0.8956 | - | - | -
50$\%$ | $5\leq n\leq 15$ | SS + JO | - | 0.8790 | 0.9003 | 0.8862
75$\%$ | $5\leq n\leq 15$ | Vanilla U-Net | 0.8496 | - | - | -
75$\%$ | $5\leq n\leq 15$ | Co-teaching [7] | 0.9023 | - | - | -
75$\%$ | $5\leq n\leq 15$ | SS | 0.9041 | - | - | -
75$\%$ | $5\leq n\leq 15$ | SS + JO | - | 0.9065 | 0.9108 | 0.9067

The experiments were conducted on the Chest X-ray dataset. We trained the
network on the samples with different levels of noisy labels and tested it by
the clean labels. Table 1 presents the segmentation performance of vanilla
U-Net (baseline) and our cascaded robust learning framework that were all
trained by noisy labels. We first trained the fully supervised vanilla U-Net
with the noisy level set to zero, which can be regarded as the upper-line
performance. Compared with the vanilla U-Net, our framework improves the
segmentation performance and achieves an average Dice of $0.925$ on the clean
annotated dataset, indicating that the sample selection stage and joint-
optimization stage can actually encourage the model to learn more distinguish
features.

For the training dataset with different level of noisy labels, we observed
that as the noise level increases, the segmentation performance of the vanilla
U-Net decreases dramatically. Compared with vanilla U-Net, the sample
selection stage (SS) can consistently improves the performance by encouraging
the model to be trained by the selected data. Through the joint optimization
(JO) stage supervised by the corrected label and original ones, the
segmentation accuracy are further improved, suggesting that our method can
effectively eliminates the effect of the noisy and gain performance by
producing correct label. In Fig. 2, we show some segmentation results under
75$\%$ noise, in which our results have higher Dice score than the baseline
method. At all the noise level, we compared our method with the state-of-the-
art noise robust method [7], which select the small loss samples according to
the prediction of peer network. The results show that our method outperform
the state-of-the-art method in all the noise level setting.

In our experiment, we also investigated the impact of the starting epoch $k$
on the performance of our method. As shown in Table 1, the joint optimization
(JO) with label correction stage is started at 20, 50, and 80 epochs,
respectively. The experimental results show that the segmentation has the best
accuracy at the intermediate epochs.

###  3.4 Analysis of Our Method

|  |
---|---|---
(a) | (b) | (c)
Figure 4: (a) The segmentation accuracy of different sample selection
criteria. (b) The segmentation accuracy of the U-Net and U-Net with only
sample selection stage on 100$\%$ noise setting. (c) The label accuracy of
labels in the original dataset , and labels corrected by the model at the end
of training.

#### Sample selection

Compared with the vanilla U-Net, our sample selection stage (SS) shows higher
segmentation accuracy under different noisy level, as shown in Table 1. To
validate the criteria of our sample selection, we conducted another experiment
that only selected the small loss sample. Fig.3(a) shows the test accuracy
with different sample selection criteria. It reveals that the test accuracy
significantly improved when considering the uncertainty in the selection
stage. To further validate the effectiveness of our method at the sample
selection stage, we applied our method on training dataset with 100$\%$ noise
and noise level $n=5,20$. Under this setting, the sample selection stage shows
worse segmentation accuracy than vanilla U-Net, because no clean sample can be
selected. The results decreased due to the low sample utilization efficiency.

#### Joint optimization

To analyze the contribution of the joint optimization stage, we explore the
label accuracy with and without the stage of joint optimization and label
correction. We calculated the Dice coefficient of the initial noisy label
($\hat{y}$) and the corrected label ($\bar{y}$) of the final model at the end
of training. Fig.3(c) shows the overall accuracy. We see that the label
quality is improved by the scheme of joint optimization and label correction.

We also investigated the impact of the weight factor $\alpha$ by adjusting
$\alpha$ between 0 and 1. If $\alpha=0$, the network is trained only using the
noisy labels without correction, and $\alpha=1$ represents the network
discards the original label and only use the corrected labels. The test
accuracy with different $\alpha$ under 50$\%$ noise ratio is shown in Fig
3.(d). Results show that model trained jointly by the corrected label
$\bar{y}$ and original label $\hat{y}$ achieves the best performance when
$\alpha=0.5$. $\alpha=1$ leads to sub-optimal performance as it might correct
some hard samples, and thus reduced the network generalization ability.

##  4 Conclusion

In this paper, we present a novel Cascaded Robust Learning framework for the
segmentation of noisy labeled chest x-ray images. Our method consists of two
stages: sample selection stage, and the stage of joint optimization with label
correction. In the first stage, the clean annotated samples are selected for
network updating, so that the influence of noisy sample can be interactively
eliminated in the three networks. In the second stage, the label correction
module work together with the joint optimization scheme to revise the
imperfect labels. Thus the training of whole network is supervised by the
corrected labels and the original ones. Compared with other state-of-the-art
models, our cascaded robust learning framework keeps high robustness when the
training data contains various noisy labels. Experimental results on the
benchmark dataset demonstrate that our network outperforms other methods on
segmentation tasks and achieves very competitive results on the noisy-labels
dataset.

## References

  * [1] Bejnordi, B.E., Veta, M., Van Diest, P.J., Van Ginneken, B., Karssemeijer, N., Litjens, G., Van Der Laak, J.A., Hermsen, M., Manson, Q.F., Balkenhol, M., et al.: Diagnostic assessment of deep learning algorithms for detection of lymph node metastases in women with breast cancer. Jama 318(22), 2199–2210 (2017)
  * [2] Berthelot, D., Carlini, N., Goodfellow, I., Papernot, N., Oliver, A., Raffel, C.A.: Mixmatch: A holistic approach to semi-supervised learning. In: Advances in Neural Information Processing Systems. pp. 5050–5060 (2019)
  * [3] Candemir, S., Jaeger, S., Palaniappan, K., Musco, J.P., Singh, R.K., Xue, Z., Karargyris, A., Antani, S., Thoma, G., McDonald, C.J.: Lung segmentation in chest radiographs using anatomical atlases with nonrigid registration. IEEE transactions on medical imaging 33(2), 577–590 (2013)
  * [4] De Fauw, J., Ledsam, J.R., Romera-Paredes, B., Nikolov, S., Tomasev, N., Blackwell, S., Askham, H., Glorot, X., O’Donoghue, B., Visentin, D., et al.: Clinically applicable deep learning for diagnosis and referral in retinal disease. Nature medicine 24(9), 1342–1350 (2018)
  * [5] Esteva, A., Kuprel, B., Novoa, R.A., Ko, J., Swetter, S.M., Blau, H.M., Thrun, S.: Dermatologist-level classification of skin cancer with deep neural networks. Nature 542(7639), 115 (2017)
  * [6] Goldberger, J., Ben-Reuven, E.: Training deep neural-networks using a noise adaptation layer. In: ICLR (2017)
  * [7] Han, B., Yao, Q., Yu, X., Niu, G., Xu, M., Hu, W., Tsang, I., Sugiyama, M.: Co-teaching: Robust training of deep neural networks with extremely noisy labels. In: Advances in neural information processing systems. pp. 8527–8537 (2018)
  * [8] Jaeger, S., Karargyris, A., Candemir, S., Folio, L., Siegelman, J., Callaghan, F., Xue, Z., Palaniappan, K., Singh, R.K., Antani, S., et al.: Automatic tuberculosis screening using chest radiographs. IEEE transactions on medical imaging 33(2), 233–245 (2013)
  * [9] Jiang, L., Zhou, Z., Leung, T., Li, L.J., Fei-Fei, L.: MentorNet: Regularizing very deep neural networks on corrupted labels. In: ICML (2018)
  * [10] Le, H., Samaras, D., Kurc, T., Gupta, R., Shroyer, K., Saltz, J.: Pancreatic cancer detection in whole slide images using noisy label annotations. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 541–549. Springer (2019)
  * [11] Ren, M., Zeng, W., Yang, B., Urtasun, R.: Learning to reweight examples for robust deep learning. In: ICML (2018)
  * [12] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [13] Setio, A.A.A., Traverso, A., De Bel, T., Berens, M.S., van den Bogaard, C., Cerello, P., Chen, H., Dou, Q., Fantacci, M.E., Geurts, B., et al.: Validation, comparison, and combination of algorithms for automatic detection of pulmonary nodules in computed tomography images: the luna16 challenge. Medical image analysis 42, 1–13 (2017)
  * [14] Shu, Y., Wu, X., Li, W.: Lvc-net: Medical image segmentation with noisy label based on local visual cues. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 558–566. Springer (2019)
  * [15] Stirenko, S., Kochura, Y., Alienin, O., Rokovyi, O., Gordienko, Y., Gang, P., Zeng, W.: Chest x-ray analysis of tuberculosis by deep learning with segmentation and augmentation. In: 2018 IEEE 38th International Conference on Electronics and Nanotechnology (ELNANO). pp. 422–428. IEEE (2018)
  * [16] Xue, C., Dou, Q., Shi, X., Chen, H., Heng, P.A.: Robust learning at noisy labeled medical images: Applied to skin lesion classification. In: ISBI (2019)
  * [17] Zhang, C., Bengio, S., Hardt, M., Recht, B., Vinyals, O.: Understanding deep learning requires rethinking generalization. In: ICLR (2017)
  * [18] Zhu, H., Shi, J., Wu, J.: Pick-and-learn: Automatic quality evaluation for noisy-labeled image segmentation. arXiv preprint arXiv:1907.11835 (2019)