# Active queue management
for alleviating Internet congestion
via a nonlinear differential equation
with a variable delay

###### Abstract

Active Queue Management (AQM) for mitigating Internet congestion has been
addressed via various feedback control syntheses, especially P, PI, and PID
regulators, by using a linear approximation where the “round trip time”, i.e.,
the delay, is assumed to be constant. This constraint is lifted here by using
a nonlinear modeling with a variable delay, introduced more than 20 years ago.
This delay, intimately linked to the congestion phenomenon, may be viewed as a
“flat output.” All other system variables, especially the control variable,
i.e., the packet loss ratio, are expressed as a function of the delay and its
derivatives: they are frozen if the delay is kept constant. This flatness-like
property, which demonstrates the mathematical discrepancy of the linear
approximation adopted until today, yields also our control strategy in two
steps: Firstly, designing an open-loop control, thanks to straightforward
flatness-based control techniques, and secondly, closing the loop via Model-
Free Control (MFC) in order to take into account severe model mismatches,
like, here, the number of TCP sessions. Several convincing computer
simulations, which are easily implementable, are presented and discussed.

> Hugues Mounier1, Cédric Join2,5, Emmanuel Delaleau3 and Michel Fliess4,5,∗
>
>
> 1 L2S (CNRS, UMR 8506), Université Paris-Saclay,
> Centrale Supélec, 3 rue Joliot-Curie, 91192 Gif-sur-Yvette, France
><EMAIL_ADDRESS>
>  2 CRAN (CNRS, UMR 7039), Université de Lorraine,
> Campus Aiguillettes, B.P. 70239, Vandœuvre-lès-Nancy, 54506, France
><EMAIL_ADDRESS><EMAIL_ADDRESS>
>  3 ENI Brest, UMR CNRS 6027, IRDL, F–29 200, France
><EMAIL_ADDRESS>
>  4 LIX (CNRS, UMR 7161), École polytechnique, 91128 Palaiseau, France
><EMAIL_ADDRESS><EMAIL_ADDRESS>><EMAIL_ADDRESS>
>  5 AL.I.E.N., 7 rue Maurice Barrès, 54330 Vézelise, France
>
>  ∗ Corresponding author.

Highlights:

  * •

In order to mitigate Internet congestion, this work is among the first ones to
use a $20$ years old modeling via a nonlinear differential equation with a
variable delay, where a new flatness-like property is encountered: the delay
is a flat output. Combining flatness-based open-loop control and closed-loop
control via the intelligent proportional controller deduced from model-free
control yields easily implementable and convincing computer experiments which
display a remarkable robustness with respect to large uncertainties on the
number of TCP connections.

  * •

The above nonlinear modeling has mainly been employed until today to derive
time-invariant linear delay approximate systems, which are quite popular, not
only for investigating control-theoretic questions but also for computer
experiments. The flatness-like property of the nonlinear model shows that
freezing the delay implies that all other system variables, including the
control one, are kept constant. The validity of the linear approximations is
therefore questioned.

Keywords: Internet congestion, active queue management, flatness-based
control, model-free control, intelligent proportional control, delay,
nonstandard analysis, time series, prediction.

##  1 Introduction

In order to alleviate Internet congestion, an _active queue management_ is a
dropping packets policy inside a router buffer yielding a corresponding queue
length management (see, e.g., [Adams(2013)], [Varma(2015)], [Grazia et
al.(2017)], [Hotchi(2021)] for surveys and comparisons). It is often related
to various control techniques and should perhaps be viewed, according to
[Varma(2015)], as “the largest human-made feedback-controlled system in the
world.” A modeling of the most popular _transmission control protocol_ (_TCP_)
has been derived more than twenty years ago in [Misra et al.(2000)] and
[Hollot et al.(2002)] via some relationship with fluid mechanics. It is a
nonlinear system of differential equations with a time-dependent delay, where
the control variable is the packet loss ratio. Although this work is much
cited, it seems, to the best of our knowledge, that almost only linear
approximations with constant delays have been exploited to propose various
applicable AQM techniques (see, however, [Barbera et al.(2010), Belamfedel
Alaoui et al.(2018), Li and Peng(2022)]). Let us restrict our short review to
a few examples where this approximation has been employed:

  * •

The familiar _random early detection_ (_RED_) algorithm, which was invented by
[Floyd and Jacobson(1993)], has been commented by [Hollot et al.(2001)], [Ryu
et al.(2005)].

  * •

The well-known _proportional-integral enhanced_ (PIE) controller is introduced
by [Pan et al.(2013)].

  * •

New algorithms are initiated by [Bisoy and Pattnaik(2017)], [Hotchi et
al.(2020)], [Bisoy et al.(2021)], [Hotchi(2021)], [Hotchi and Kubo(2022)].

Our work relies on a remarkable attribute of the above mentioned nonlinear
system: it should be called _flat_ in a more or less analogous sense of
[Fliess et al.(1995), Fliess et al.(1999)]. The delay is a flat output, the
queue length another one. This means that the knowledge of the delay time
variation, or of the queue length, determines all the other system variables,
including the control one.

In particular, freezing the delay implies at once that all the other system
variables are constant. This property thus questions the frequent use that
seemed until today so self-evident, both for simulation and control purposes,
of the linear approximations, where the delay is assumed to be constant and
not the other system variables (see, e.g., [Alli-Oke(2022)], and references
therein).

This paper shows that appropriate control-theoretic tools do exist for
handling the nonlinear modeling:

  * •

Open-loop control strategies are deduced at once by exploiting this flatness-
like property of the nonlinear modeling. We here choose to regulate the delay:

    * –

It is obviously related to congestion.

    * –

_Controlling queue Delay_ (_CoDel_) ([Nichols and Jacobson(2012)]), which has
become a popular setting, also puts delay control on the forefront but via a
completely different viewpoint.

  * •

In order to counteract the unavoidable model mismatches (see, e.g., [Xu et
al.(2015)] for a summary of the shortcomings of the above nonlinear modeling)
and disturbances, we follow [Villagra and Herrero-Perez(2012)], [Fliess et
al.(2021)], [Join et al.(2022a)] by closing the loop via _model-free_ control
in the sense of [Fliess and Join(2013), Fliess and Join(2022)]. The presence
of a delay requires a predictor which cannot, here, be the celebrated Smith’s
predictor ([Smith(1957)]), because the latter is model-based (see, e.g., [Deng
et al.(2022)] for a recent survey). We thus adapt here the viewpoint of
[Hamiche et al.(2019)] for studying _supply chain management_ (see also [Join
et al.(2022b)]). This is achieved by removing the unpredictable _quick
fluctuations_ via a theorem due to [Cartier and Perrin(1995)] which is
expressed in the language of _nonstandard analysis_. Note that this result has
led to a new understanding of time series (see, e.g, [Fliess et al.(2018)],
and references therein).

Our paper is organized as follows. After showing that the delay or the queue
length can be viewed as a flat output, Section 2 explains the inherent
weakness of linear approximations with a constant delay. Section 3 recalls the
basic facts of model-free control, which has already been used successfully
many times. In order to take the delay into account in the model-free
approach, Section 4 exploits techniques stemming from nonstandard analysis.
Numerical simulations are presented in Section 5: they show that the model
mismatch on the number of TCP sessions is well compensated by the closed-loop
control without any clear-cut superiority of the techniques developed in
Section 4. Various questions are raised in Section 6.

##  2 Some consequences of the nonlinear modeling

###  2.1 Flatness

The nonlinear TCP/AQM network model ([Misra et al.(2000), Hollot et
al.(2002)]) reads

$\displaystyle\dot{W}(t)$
$\displaystyle=\frac{1}{R(t)}-\frac{W(t)W(t-R(t))}{2[R(t-R(t))]}u(t-R(t))$
(1a) $\displaystyle\dot{Q}(t)$ $\displaystyle=\frac{W(t)}{R(t)}N(t)-C(t)$ (1b)

where

  * •

$W(t)>0$ is the length of the TCP window;

  * •

$R(t)>0$ is the _round trip time_ (_RTT_) which appears as a time-dependent
delay in Equation (1);

  * •

$Q(t)>0$ is the queue length;

  * •

the control variable is the dropping packet policy $u(t)$, $0\leqslant
u(t)\leqslant 1$: it is called the _packet loss ratio_ ([Alli-Oke(2022)]), or,
as often in the literature, the _packet drop probability_ ;

  * •

$C(t)>0$ is the bottleneck link capacity;

  * •

$N(t)>0$ is the number of TCP sessions. It plays the role of external
disturbance.

In practice $C(t)$ and $N(t)$ are piecewise constant. Thus
$\dot{C}(t)=\dot{N}(t)=0$, with the exception of a finite number of points on
any finite-time interval. The RTT $R(t)$ and the queue length $Q(t)$ are
related by a simple affine relation

$R(t)=T+\frac{Q(t)}{C(t)}$ (2)

where $T$ is the round trip propagation time. Assume, in the following
computations, $\dot{C}(t)=\dot{N}(t)=0$ on some time interval, Set $C(t)=C$,
$N(t)=N$, where $C$ and $N$ are constant. Then Equation (2) yields
$\dot{Q}(t)=C\dot{R}(t)$. Equation (1b) becomes

$\dot{R}(t)=\frac{NW(t)}{CR(t)}-1$ (3)

It yields

$W(t)=\frac{CR(t)(\dot{R}(t)+1)}{N}$ (4)

Thus Equations (4) and (1a) show that $W(t)$ and $u(t-R(t))$ depend on $R(t)$
and its first and second order derivatives. In other words, we may call System
(1) _flat_ and the RTT $R(t)$ a _flat output_ (compare with [Mounier et
al.(2003)]). Equation (2) shows that the queue length $Q(t)$ is another flat
output.

###### Remark 1

Classic flatness has been formally defined in [Fliess et al.(1995)] via
differential algebra and in [Fliess et al.(1999)] via differential geometry of
infinite jets and prolongations. Combining differential and difference
algebras (see, e.g., [Cohn(1970)]) permits a precise definition of flatness
for nonlinear systems with constant delays ([Mounier and Rudolph(1998),
Mounier and Rudolph(2008)]). Such a setting does not however work with a
variable delay such as $R(t)$, since the time derivation $\frac{d}{dt}$ and
the time shift $t\mapsto t-R(t)$ with a time-varying quantity do not commute.
Before developing an adequate general mathematical formalism for our example,
it may be wise to wait for other concrete case-studies. It would open a path
to a new understanding of nonlinear systems with variable delays.

###  2.2 Critical appraisal of the linear approximation

See [Alli-Oke(2022)] for a nice survey on linear approximations. Let $u_{0}$,
$Q_{0}$, $W_{0}$, $R_{0}$ be the numerical values of $u(t)$, $Q(t)$, $W(t)$,
$R(t)$ at an operating, or equilibrium, point. Contrarily to the variables
$u$, $Q$, $W$, the delay $R$ is kept frozen at the value $R_{0}$: the delay in
the linear approximation is constant. It is obvious that such an assumption
contradicts Section 2.1, where Equations (2), (4) and (1a) show that $Q(t)$,
$W(t)$ and $u(t)$ become constant when $R(t)$ is constant. This fact is
casting some doubt not only about AQM via such approximations, but also on the
computer simulations, which rely on it (see [Alli-Oke(2021), Alli-Oke(2022)],
and the references therein).

###### Remark 2

Define the control variable $\delta u(t)=u(t)-u_{0}$ and the output variable
$\delta Q(t)=Q(t)-Q_{0}$. They are often related in the literature (see, e.g.,
[Alli-Oke(2022)]) by the time-invariant linear delay system defined by the
transfer function

$-\frac{(2N\frac{W_{0}}{2})^{3}e^{-R_{0}s}}{(R_{0}s+1)(\frac{W_{0}R_{0}}{2}s+1)}$

where the number $N$ of sessions is assumed to be constant. A system defined
by such a transfer function is sometimes called _quasi-finite_ ([Fliess et
al.(2002)]). The output $\delta Q$ is said to be flat, or _basic_ ([Fliess et
al.(2002)]).

##  3 Closed-loop control via model-free control without delay: a short review

###  3.1 Ultra-local model

Consider a single-input single-output (SISO) nonlinear system ($\Sigma$).
Denote by $\mathfrak{u}(t)$ (resp. $\mathfrak{y}(t)$) the control (resp.
output) variable. It has been demonstrated ([Fliess and Join(2013)]) via
elementary techniques from functional analysis and differential algebra that
the often poorly known modeling of ($\Sigma$) may be replaced, if some quite
weak assumptions are satisfied, by an _ultra-local model_ :

$\mathfrak{y}^{(\nu)}=F+\alpha\mathfrak{u}$ (5)

where $\alpha\in\mathbb{R}$ is chosen by the practitioner such that
$\alpha\mathfrak{u}$ and $\mathfrak{y}^{(\nu)}$ are of the same order of
magnitude: it does not need to be precisely known. Numerous successful
applications (see, e.g., references in [Fliess and Join(2013), Fliess and
Join(2022)]) have shown that $\nu=1$ in Equation (5) yields most often a
convenient ultra-local model:

$\dot{\mathfrak{y}}=F+\alpha\mathfrak{u}$ (6)

The following comments are useful:

  * •

Equation (6) is only valid during a short time lapse: it must be continuously
updated.

  * •

$F$ is _data-driven_ , i.e., it is estimated via the knowledge $\mathfrak{u}$
and $\mathfrak{y}$ ([Fliess and Join(2013)]):

$F_{\text{est}}(t)=-\frac{6}{\tau^{3}}\int_{t-\tau}^{t}\left[(\tau-2\sigma){\mathfrak{y}}(\sigma)+\alpha\sigma(\tau-\sigma){\mathfrak{u}}(\sigma)\right]d\sigma$
(7)

The quantity $\tau>0$ may be chosen to be quite “small.” The above integral,
which is a low pass filter, may, in practice, be replaced by a classic digital
filter.

  * •

$F$ subsumes not only the unknown structure of the system, which most of the
time is nonlinear, but also any external disturbance.

###  3.2 Intelligent controllers and local stability

The loop is closed with the following _intelligent proportional controller_
([Fliess and Join(2013)]), or _iP_ ,

$\mathfrak{u}=-\frac{F_{\text{est}}-\dot{\mathfrak{y}}^{\star}+K_{P}e}{\alpha}$
(8)

where:

  * •

$\mathfrak{y}^{\star}$ is the reference trajectory of the output,

  * •

$e=\mathfrak{y}-\mathfrak{y}^{\star}$ is the tracking error,

  * •

$K_{P}\in\mathbb{R}$ is a tuning gain.

Combining Equations (6) and (8) yields $\dot{e}+K_{P}e=F-F_{\text{est}}$. If
the estimate $F_{\text{est}}$ is “good,” i.e., if $F-F_{\text{est}}\approx 0$,
then $\lim_{t\to+\infty}e(t)\approx 0$, if, and only if, $K_{P}>0$.

##  4 Closed-loop control via model-free control with delay

###  4.1 Ultra-local model with delay

Set $v(t)=u(t-R(t))$ in Equation (1):

$\displaystyle\dot{W}(t)$
$\displaystyle=\frac{1}{R(t)}-\frac{W(t)W(t-R(t))}{2[R(t-R(t))]}v(t)$
$\displaystyle\dot{Q}(t)$ $\displaystyle=\frac{W(t)}{R(t)}N(t)-C(t)$

This trivial change of variable shows that the techniques from [Fliess and
Join(2013)] remain valid for introducing the ultra-local model with a time-
varying delay

$\dot{\mathfrak{y}}(t)=\mathfrak{F}+\alpha\mathfrak{u}(t-R(t))$ (10)

where $\mathfrak{F}$ plays the same rôle as $F$ in Equation (6). Equation (7)
becomes

$\mathfrak{F}_{\text{est}}(t)=-\frac{6}{\tau^{3}}\int_{t-\tau}^{t}\left[(\tau-2\sigma)\mathfrak{y}(\sigma)+\alpha\sigma(\tau-\sigma)\mathfrak{u}(\sigma-R(\sigma)\right]d\sigma$
(11)

###  4.2 Prediction via time series

####  4.2.1 Time series and the Cartier-Perrin theorem

Consider the time interval $[0,1]\subset\mathbb{R}$. Introduce as often in
_nonstandard analysis_ (see [Robinson(1974)], [Diener and Diener(1995)],
[Lobry and Sari(2008)]) the infinitesimal sampling of $[0,1]$:
${\mathfrak{T}}=\\{0=t_{0}<t_{1}<\dots<t_{\nu}=1\\}$ where $t_{i+1}-t_{i}$,
$0\leqslant i<\nu$, is infinitesimal, i.e., “very small”. A _time series_
$X(t)$ is a function $X:{\mathfrak{T}}\rightarrow\mathbb{R}$.

A time series ${\mathcal{X}}:{\mathfrak{T}}\rightarrow\mathbb{R}$ is said to
be quickly fluctuating, or oscillating, if the integral
$\int_{A}{\mathcal{X}}dm$ is infinitesimal, i.e., very small, for any
_appreciable_ interval, i.e., an interval which is neither “very small” nor
“very large”.

According to a theorem due to [Cartier and Perrin(1995)], the following
additive decomposition holds for any time series $X$, which satisfies a weak
integrability condition,

$X(t)=E(X)(t)+X_{\tiny{\rm fluctuation}}(t)$ (12)

where

  * •

the _mean_ , or _trend_ , $E(X)$ is “quite smooth”;

  * •

$X_{\tiny{\rm fluctuation}}$ is quickly fluctuating.

The decomposition (12) is unique up to an additive infinitesimal: It means
that the two terms on the right handside of Equation (12) are unique up to a
“very small” additive quantitity.

####  4.2.2 Derivative estimate

Let us start with a polynomial time function of degree $1$

$p_{1}(\tau)=a_{0}+a_{1}\tau$

where $\tau\geqslant 0$, $a_{0},a_{1}\in\mathbb{R}$. Operational calculus
(see, e.g., [Yosida(1984)]) with respect to the variable $\tau$, permits to
express $p_{1}$ as

$P_{1}={a_{0}}/{s}+{a_{1}}/{s^{2}}$

Multiply both sides by $s^{2}$:

$s^{2}P_{1}=a_{0}s+a_{1}$ (13)

Take the derivative of both sides with respect to $s$, which corresponds in
the time domain to the multiplication by $-\tau$:

$s^{2}\frac{dP_{1}}{ds}+2sP_{1}=a_{0}$ (14)

The coefficients $a_{0},a_{1}$ are obtained via the triangular system of
linear equations (13)-(14). We get rid of the time derivatives, i.e., of
$sP_{1}$, $s^{2}P_{1}$, and $s^{2}\frac{dP_{1}}{ds}$, by multiplying both
sides of Equations (13)-(14) by $s^{-n}$, $n\geq 2$. The corresponding
iterated time integrals are low pass filters which attenuate the corrupting
noises. A quite short time window is sufficient for obtaining accurate values
of $a_{0}$, $a_{1}$.

###### Remark 3

See [Mboup et al.(2009)] and [Othmane et al.(2021), Othmane et al.(2022)] for
more details. Note also that estimating derivatives via integrals seems to
have been first introduced by [Lanczos(1956)].

####  4.2.3 Prediction

Set the following forecast $X_{\text{forecast}}(t+\Delta T)$, where $\Delta
T>0$ is not too “large”,

$X_{\text{forecast}}(t+\Delta
T)=E(X)(t)+\left[\frac{dE(X)(t)}{dt}\right]_{e}\Delta T$ (15)

where $E(X)(t)$ and $\left[\frac{dE(X)(t)}{dt}\right]_{e}$ are estimated like
$a_{0}$ and $a_{1}$ above. Let us stress that what we predict is the mean and
not the quick fluctuations.

###### Remark 4

The above construction is obviously reminiscent of the _sliding window_
techniques in the applied literature on time series (see, e.g.,
[Mélard(2008)]).

Note that estimating $a_{0}$ and $a_{1}$ yields respectively the mean and the
derivative.

####  4.2.4 Local closed-loop stability

Equation (10) may be rewritten as

$\dot{\mathfrak{y}}(t+S(t))=\mathfrak{F}_{\text{forecast}}(t+S(t))+\alpha\mathfrak{u}(t)$
(16)

where the advance $S(t)>0$ is defined by

$S(t)=\min\left\\{\,\tau\ |\ \tau-R(\tau)=t\,\right\\}-t$ (17)

If $R(t)$ is “slowly” varying, it is clear that $R(t)$ and $S(t)$ remain
close. Evaluating $S(t)$ however requires a prediction of $R(t)$. Replace
therefore in Equation (17) $R(t)$ by the reference trajectory $R^{\star}(t)$.
It yields

$S^{\star}(t)=\min\\{\tau|\tau-R^{\star}(\tau)=t\\}-t$ (18)

Equation (8) then becomes

$\mathfrak{u}(t)=-\frac{\mathfrak{F}_{\text{forecast}}(t+S^{\star}(t))-\dot{\mathfrak{y}}^{\star}(t+S^{\star}(t))+K_{P}e{(t+S^{\star}(t))}}{\alpha}$
(19)

where:

  * •

the forecast of $\mathfrak{F}$ is obtained via Formulae (11) and (15);

  * •

$e(t+S^{\star}(t))=\mathfrak{y}_{\text{forecast}}(t+S^{\star}(t))-\mathfrak{y}^{\star}(t+S^{\star}(t))$,
where $\mathfrak{y}^{\star}$ is the reference trajectory and $e$ is the
tracking error;

  * •

$\mathfrak{y}_{\text{forecast}}(t+S^{\star}(t))=\mathfrak{z}(t+S^{\star}(t))$
is obtained via the linear differential equation

$\dot{\mathfrak{z}}(\tau)=\mathfrak{F}_{\text{forecast}}(\tau)+\alpha\mathfrak{u}(\tau-R^{\star}(\tau))\quad
t\leqslant\tau\leqslant t+S(t)$

  * •

$K_{P}\in\mathbb{R}$ is the tuning gain.

It yields

$\displaystyle\dot{e}(t+S^{\star}(t))+K_{P}e(t+S^{\star}(t))=\mathfrak{F}(t+S^{\star}(t))-\mathfrak{F}_{\text{forecast}}(t+S^{\star}(t))$
(20)

Local stability is ensured, i.e.,
$\lim_{t\to+\infty}\mathfrak{y}(t+S^{\star}(t))\approx\mathfrak{y}_{\text{forecast}}(t+S^{\star}(t))$
, if

  * •

$K_{P}>0$,

  * •

the forecast is “good,” i.e.,
${\mathfrak{F}(t+S^{\star}(t))-\mathfrak{F}_{\text{forecast}}(t+S^{\star}(t))}\approx
0$.

##  5 Computer simulations111Contact C. Join<EMAIL_ADDRESS>for the simulation codes.

###  5.1 Various situations

Introduce the following control settings:

  1. 1.

Reference trajectory and nominal control: The choice of a reference trajectory
$R^{\star}(t)$ for the delay $R(t)$ yields at once via Section 2.1 an open-
loop nominal control $u^{\star}(t)$ for $u(t)$, i.e.,

$u^{\star}(t-R^{\star}(t))=2\left(\frac{1}{R^{\star}(t)}-\dot{W}^{\star}(t)\right)\left(\frac{R^{\star}(t-R^{\star}(t))}{W^{\star}(t)W^{\star}(t-R^{\star}(t))}\right)$

or

$u^{\star}(t)=2\left(\frac{1}{R^{\star}(t+S^{\star}(t))}-\dot{W}^{\star}(t+S^{\star}(t))\right)\left(\frac{R^{\star}(t)}{W^{\star}(t+S^{\star}(t))W^{\star}(t)}\right)$

and

$W^{\star}(t)=\frac{C(t)R^{\star}(t)(\dot{R^{\star}}(t)+1)}{N_{0}}$

where $N_{0}=60$ and $C=3000$.

  2. 2.

Open-loop control (OL): Inject $u^{\star}$ in Equation (1) with $N(t)=N_{0}$.

  3. 3.

Closing the loop via an iP without delay (iP): In Equations (6), (7), (8), set
$\mathfrak{y}(t)=e(t)=R(t)-R^{\star}(t)$, $\mathfrak{y}^{\star}(t)=0$,
$\mathfrak{u}(t)=\Delta u(t)=u(t)-u^{\star}(t)$. The iP (8) becomes

$\Delta u=-\frac{F_{\text{est}}+K_{P}e}{\alpha}$

where $\alpha=-1000$, $K_{P}=1$. Consider the estimation of $F_{\text{est}}$
via Formula (7) as a classic _finite impulse response_ (_FIR_) (see, e.g.,
[Rabiner and Gold(1975)]). We thus apply a control $u(t)$ of the form

$u(t)=u^{\star}(t)-\frac{F_{\text{est}}+K_{P}e}{\alpha}$

where the first part incorporates our knowledge of the system, and the second
one deals with perturbations, model imperfections and unknown dynamics. It is
clear that the practical implementation, which has already been achieved
successfully a number of times, is straightforward.

  4. 4.

Closing the loop via an iP with delay (iPWD): In Equations (10), (11), (19),
set as above $\mathfrak{y}(t)=e(t)=R(t)-R^{\star}(t)$,
$\mathfrak{y}^{\star}(t)=0$, $\mathfrak{u}(t)=\Delta u(t)=u(t)-u^{\star}(t)$.
The iP (19) becomes

$\Delta
u(t)=-\frac{\mathfrak{F}_{\text{forecast}}(t+S^{\star}(t))+K_{P}e(t+S^{\star}(t))}{\alpha}$

where, as above, $\alpha=-10$, $K_{P}=1$. The calculations related to
predictions are detailed in [Fliess et al.(2018)]. The sequel is similar to
the case without delay.

###  5.2 Scenarios

We illustrate our control laws through two different scenarios, where the
mismatch is the number of TCP sessions. The first scenario corresponds to
normal operation: small variation of R. The second scenario represents an exit
from a congestion situation and corresponds to a large variation of R. The
command is designed with a known value of the number of connections N but
operates with a N constant piecewise: in our simulation N first goes from 60
to 70 then goes down to 50. The number of connection N plays the role of an
external disturbance. It is moreover this effect that the open loop curves
Figs.1 and 4 show where the trajectory deviates from the reference one when
the number of connections changes.

  1. 1.

Scenario $1$: $0.25$s $\leqslant R(t)\leqslant 0.3$s, $50\leqslant
N(t)\leqslant 70$.

  2. 2.

Scenario $2$: $0.3$s $\leqslant R(t)\leqslant 0.7$s, $50\leqslant
N(t)\leqslant 70$.

Table 1: Simulations and Figures Scenarios | OL | iP | iPWD
---|---|---|---
1 | Figure 1 | Figure 2 | Figure 3
2 | Figure 4 | Figure 5 | Figure 6

Those Figures tell us that:

  * •

the mismatch is well compensated by iPs with or without delay,

  * •

iPs with or without delay exhibit very similar behaviors.

Our results are especially well displayed in Figure 7. Indeed, we can see that
the behavior of the tracking errors are totally comparable. In other words the
iP with delay seems to be useless!

##  6 Conclusion

It has been shown that

  * •

the linear constant delay approximations, which also play a key rôle in
computer simulations, contradict the more complete nonlinear modeling;

  * •

control-theoretic tools are available for an active queue management via this
nonlinear modeling.

Many points remain of course to be addressed:

  * •

Other mismatches and external disturbances ought to be examined: noisy
measurements, abrupt changes of the round trip time, … Would the _intelligent
proportional-derivative controller_ (_iPD_) advocated by [Fliess and
Join(2022)] be helpful? See, e.g., [Sun et al.(2003)] and [Ryu et al.(2005)]
for results with classic PD controllers.

  * •

The simulations in Section 5.2 indicate the futility of an iP with delay in
order to compensate a model mismatch. Without a precise mathematical analysis,
it is not clear whether this property is always valid. Let us suggest
nevertheless that the open-loop nonlinear control, where the delay is taken
into account, is doing the job!

  * •

The coefficient $\alpha$ in Equations (6) and (10), which does not need to be
determined precisely, is obtained via trials and errors. A more subtle
estimation would be welcome. Let us also add that important variations of some
quantities like the number $N(t)$ of TCP connections might necessitate the
introduction of a time-varying $\alpha$ (see [Gédoin et al.(2011), Moreno-
Gonzales et al.(2022)] for first results in other engineering domains).

  * •

Many network simulation for investigating Internet congestion (see, e.g.,
[Riley and Henderson(2010)], [Alli-Oke(2021), Alli-Oke(2022)], and references
therein) seem to have employed time-invariant linear delay systems (see
Section 2.2). It should therefore be most rewarding to develop and integrate
the tools of this paper.

(a) Control (–) and nominal control (- -)  (b) $R$ (–) and reference
trajectory (- -)

(c) TCP Window (d) Queue length (e) Number of connections

Figure 1: Scenario 1 – OL

(a) Control (–) and nominal control (- -)  (b) $R$ (–) and reference
trajectory (- -)

(c) TCP Window (d) Queue length (e) Number of connections

Figure 2: Scenario 1 – iP

(a) Control (–) and nominal control (- -)  (b) $R$ (–) and reference
trajectory (- -)

(c) TCP Window (d) Queue length (e) Number of connections

Figure 3: Scenario 1 – iPWD

(a) Control (–) and nominal control (- -)  (b) $R$ (–) and reference
trajectory (- -)

(c) TCP Window (d) Queue length (e) Number of connections

Figure 4: Scenario 2 – OL

(a) Control (–) and nominal control (- -)  (b) $R$ (–) and reference
trajectory (- -)

(c) TCP Window (d) Queue length (e) Number of connections

Figure 5: Scenario 2 – iP

(a) Control (–) and nominal control (- -)  (b) $R$ (–) and reference
trajectory (- -)

(c) TCP Window (d) Queue length (e) Number of connections

Figure 6: Scenario 2 – iPWD

(a) Scenario 1 iP (–), iPWD (- .) and zero line (- -)  (b) Scenario 2 iP (–),
iPWD (- .) and zero line (- -)

Figure 7: Tracking errors

## References

  * [Adams(2013)] Adams, R., 2013.  Active queue management: A survey.  IEEE Communications Surveys & Tutorials 15, 1425–1476.
  * [Alli-Oke(2021)] Alli-Oke R.O., 2021.  Repository.  https://github.com/droa28/AQM-Simulations
  * [Alli-Oke(2022)] Alli-Oke, R.O., 2022.  On the validity of numerical simulations for control-theoretic AQM schemes in computer networks.  Mathematics and Computers in Simulation 193, 466–480.
  * [Barbera et al.(2010)] Barbera M., Lombardo A., Panarello C., Schembra G., 2010.  Queue stability analysis and performance evaluation of a TCP-compliant window management mechanism.  IEEE/ACM Transactions on Networking 18, 1275–1288.
  * [Belamfedel Alaoui et al.(2018)] Belamfedel Alaoui S., Tissir E.H., Chaibi N., 2018.  Active queue management based feedback control for TCP with successive delays in single and multiple bottleneck topology.  Comput. Communicat. 117, 58–70.
  * [Bisoy and Pattnaik(2017)] Bisoy, S.K., Pattnaik, P.K., 2017\.  Design of feedback controller for TCP/AQM networks.  Engineering Science and Technology 20, 116–132.
  * [Bisoy et al.(2021)] Bisoy, S.K., Pattnaik, P.K., Sain, M., Jeong, D.U., 2021\.  A self-tuning congestion tracking control for TCP/AQM network for single and multiple bottleneck topology.  IEEE Access 9, 27723–27735.
  * [Cartier and Perrin(1995)] Cartier, P., Perrin, Y., 1995\.  Integration over finite sets. In Diener, F., Diener, M. (Eds.), Nonstandard Analysis in Practice. Springer, Berlin, Heidelberg, pp. 185–204.
  * [Cohn(1970)] Cohn, R.M., 1970.  A difference-differential basis theorem.  Canadian Journal of Mathematics 22, 1224–1237.  doi:10.4153/CJM-1970-141-3.
  * [Deng et al.(2022)] Deng, Y., Léchappé, V., Moulay, E., Chen, Z., Liang, B., Plestan, F., Han, Q.L., 2022.  Predictor-based control of time-delay systems: a survey.  International Journal of Systems Science 53, 2496–2534.
  * [Diener and Diener(1995)] Diener, F., Diener, M., 1995\.  Tutorial. In Diener, F., Diener, M. (Eds.), Nonstandard Analysis in Practice. Springer, Berlin, Heidelberg, pp. 1–21.
  * [Fliess and Join(2013)] Fliess, M., Join, C., 2013.  Model-free control.  International Journal of Control 86, 2228–2252.
  * [Fliess and Join(2022)] Fliess, M., Join, C., 2022.  An alternative to proportional-integral and proportional-integral-derivative regulators: Intelligent proportional-derivative regulators.  International Journal of Robust and Nonlinear Control 32, 9512–9524.
  * [Fliess et al.(2021)] Fliess, M., Join, C., Moussa, K., Djouadi, S.M., Alsager, M.W., 2021.  Toward simple in silico experiments for drugs administration in some cancer treatments.  IFAC-PapersOnLine 54, 245–250.
  * [Fliess et al.(2018)] Fliess, M., Join, C., Voyant, C., 2018.  Prediction bands for solar energy: New short-term time series forecasting techniques.  Solar Energy 166, 519–528.
  * [Fliess et al.(1995)] Fliess, M., Lévine, J., Martin, P., Rouchon, P., 1995\.  Flatness and defect of non-linear systems: introductory theory and examples.  International Journal of Control 61, 1327–1361.
  * [Fliess et al.(1999)] Fliess, M., Lévine, J., Martin, P., Rouchon, P., 1999\.  A Lie-Bäcklund approach to equivalence and flatness of nonlinear systems.  IEEE Transactions on Automatic Control 44, 922–937.
  * [Fliess et al.(2002)] Fliess, M., Marquez, R., Mounier, H., 2002.  An extension of predictive control, PID regulators and Smith predictors to some linear delay systems.  International Journal of Control 75, 728–743.
  * [Floyd and Jacobson(1993)] Floyd, S., Jacobson, V., 1993\.  Random early detection gateways for congestion avoidance.  IEEE/ACM Transactions on Networking 1, 397–413.
  * [Gédoin et al.(2011)] Gédouin P.-A., Delaleau E., Bourgeot J.-M., Join C., Arbab Chirani S., Calloch S., 2011.  Experimental comparison of classical PID and model-free control: Position control of a shape memory alloy active spring.  Control Engineering Practice 19, 433–441.
  * [Grazia et al.(2017) ] Grazia C.A., Patriciello N., Klapez M., M. Casoni M., 2017.  A cross-comparison between TCP and AQM algorithms: Which is the best couple for congestion control?. 2017 14th International Conference on Telecommunications (ConTEL), Zagreb, pp. 75–82.
  * [Hamiche et al.(2019)] Hamiche, K., Fliess, M., Join, C., Abouaïssa, H., 2019\.  Bullwhip effect attenuation in supply chain management via control-theoretic tools and short-term forecasts: A preliminary study with an application to perishable inventories, in: 2019 6th International Conference on Control, Decision and Information Technologies (CoDIT), Paris, pp. 1492–1497.
  * [Hollot et al.(2001)] Hollot, C., Misra, V., Towsley, D., Gong, W.B., 2001\.  A control theoretic analysis of RED. In: Proc. IEEE INFOCOM 2001, pp. 1510–1519, vol. 3.
  * [Hollot et al.(2002)] Hollot, C., Misra, V., Towsley, D., Gong, W., 2002\.  Analysis and design of controllers for AQM routers supporting TCP flows.  IEEE Transactions on Automatic Control 47, 945–959.
  * [Hotchi(2021)] Hotchi, R.A., 2021.  Active queue management based on control-theoretic approaches for diversified communication services.  PhD Thesis, Keio University, Tokyo.
  * [Hotchi et al.(2020)] Hotchi, R., Chibana, H., Iwai, T., Kubo, R., 2020.  Active queue management supporting TCP flows using disturbance observer and smith predictor.  IEEE Access 8, 173401–173413.
  * [Hotchi and Kubo(2022)] Hotchi, R., Kubo, R., 2022.  Quality of service aware adaptive target queue length generation for active queue management.  IET Control Theory & Applications 16, 398–413.
  * [Join et al.(2022a)] Join, C., d’Onofrio, A., Fliess, M., 2022a.  Toward more realistic social distancing policies via advanced feedback control.  Automation 3, 286–301.
  * [Join et al.(2022b)] Join, C., Mounier, H., Delaleau, E., Fliess, M., 2022b.  Active queue management: First steps toward a new control-theoretic viewpoint.  International Conference on Systems and Control (ICSC’2022), Marseille, 448–453.
  * [Lanczos(1956)] Lanczos, C., 1956.  Applied Analysis.  Prentice-Hall, Englewood Cliffs.
  * [Li and Peng(2022)] Li, C., Deng, D., 2022.  Uniform stability of nonlinear systems with state-dependent delay.  Automatica 137, 110098.
  * [Lobry and Sari(2008)] Lobry, C., Sari, T., 2008.  Non-standard analysis and representation of reality.  International Journal of Control 81, 519–536.
  * [Mboup et al.(2009)] Mboup, M., Join, C., Fliess, M., 2009.  Numerical differentiation with annihilators in noisy environment.  Numerical Algorithms 50, 439–467.
  * [Mélard(2008)] Mélard, G., 2008.  Méthodes de prévision à court terme ($2^{e}$ éd.).  Ellipses & Presses Universitaires de Bruxelles, Paris & Bruxelles.
  * [Misra et al.(2000)] Misra, V., Gong, W.B., Towsley, D., 2000.  Fluid-based analysis of a network of AQM routers supporting TCP flows with an application to RED.  SIGCOMM Comput. Commun. Rev. , 151–160.
  * [Moreno-Gonzales et al.(2022)] Moreno-Gonzalez M., Artu$\tilde{\rm n}$edo A., Villagra J., Join C., Fliess M., 2022.  Speed-adaptive model-free lateral control for automated cars.  IFAC-PapersOnLine 55, 84–89.
  * [Mounier et al.(2003)] Mounier, H., Bastin, G., Vèque, V., Zitoune, L., 2003\.  Round trip time TCP tracking: a first step towards QoS pricing.  International Journal of Systems Science 34, 607–614.
  * [Mounier and Rudolph(1998)] Mounier, H., Rudolph, J., 1998\.  Flatness-based control of nonlinear delay systems: A chemical reactor example.  International Journal of Control 71, 871–890.
  * [Mounier and Rudolph(2008)] Mounier, H., Rudolph, J., 2008\.  Flatness and quasi-static state feedback in non-linear delay systems.  International Journal of Control 81, 445–456.
  * [Nichols and Jacobson(2012)] Nichols, K., Jacobson, V., 2012\.  Controlling queue delay: A modern AQM is just one piece of the solution to bufferbloat.  ACMQueue 5, 1–20.
  * [Othmane et al.(2022)] Othmane, A., Kiltz, L., Rudolph, J., 2022\.  Survey on algebraic numerical differentiation: historical developments, parametrization, examples, and applications.  International Journal of Systems Science 53, 1848–1887.
  * [Othmane et al.(2021)] Othmane, A., Rudolph, J., Mounier, H., 2021.  Systematic comparison of numerical differentiators and an application to model-free control.  European Journal of Control 62, 113–119.
  * [Pan et al.(2013)] Pan, R., Natarajan, P., Piglione, C., Prabhu, M.S., Subramanian, V., Baker, F., VerSteeg, B., 2013.  PIE: A lightweight control scheme to address the bufferbloat problem, in: IEEE 14th International Conference on High Performance Switching and Routing (HPSR), pp. 148–155.
  * [Rabiner and Gold(1975)] Rabiner, L.R., Gold, B., 1975\.  Theory and Application of Digital Signal Processing.  Prentice-Hall, Englewood Cliffs.
  * [Riley and Henderson(2010)] Riley, G.F., Henderson, T.R., 2010\.  The NS–3 network simulator.  In: Wehrle, K., Güneş, M., Gross, J. (Eds): Modeling and Tools for Network Simulation, pp. 15–34, Springer, Berlin.
  * [Robinson(1974)] Robinson, A., 1974.  Non-standard Analysis (2nd ed.).  North-Holland, Amsterdam.
  * [Ryu et al.(2005)] Ryu S., Ryu B., Jeong, Park S., 2005.  PI-PD controller for adaptive and robust active queue management for Internet congestion control.  Simulation 81, 437–459.
  * [Smith(1957)] Smith, O.J.M., 1957.  Closer control of loops with dead time.  Chemistry Engineering Progress 53, 217–219.
  * [Sun et al.(2003)] Sun J., Chen G., Ko K.-T., Chan S., Zukerman M., 2003.  PD-controller: a new active queue management scheme, GLOBECOM’03, San Francisco, 3103–3107.
  * [Varma(2015)] Varma, S., 2015.  Internet Congestion Control.  Morgan Kaufmann, Waltham.
  * [Villagra and Herrero-Perez(2012)] Villagra, J., Herrero-Perez, D., 2012\.  A comparison of control techniques for robust docking maneuvers of an AGV.  IEEE Trans. Control Systems Technology 20, 1116–1123.
  * [Xu et al.(2015)] Xu Q., Li F., Sun J., Zukerman M., 2015.  A new TCP/AQM system analysis.  J. Network Computer Applications 57, 43–60.
  * [Yosida(1984)] Yosida, K., 1984.  Operational Calculus (translated from the Japanese).  Springer, New York.


# Localized geometric moves to compute hyperbolic structures on triangulated
3-manifolds

Clément Maria  Owen Rouillé

###### Abstract

A fundamental way to study 3-manifolds is through the geometric lens, one of
the most prominent geometries being the hyperbolic one. We focus on the
computation of a complete hyperbolic structure on a connected orientable
hyperbolic 3-manifold with torus boundaries. This family of 3-manifolds
includes the knot complements.

This computation of a hyperbolic structure requires the resolution of gluing
equations on a triangulation of the space, but not all triangulations admit a
solution to the equations.

In this paper, we propose a new method to find a triangulation that admits a
solution to the gluing equations, using convex optimization and localized
combinatorial modifications. It is based on Casson and Rivin’s reformulation
of the equations. We provide a novel approach to modify a triangulation and
update its geometry, along with experimental results to support the new
method.

##  1 Introduction

A main problem of knot theory is to tell whether two knots are equivalent or
distinct. Equivalence between knots is defined by the existence of an isotopy
of the ambient space that would turn one knot into the other, i.e., a
continuous deformation of the space that preserves the entanglement.

Isotopies are too difficult to compute in practice, and practitioners use
invariants to tackle the knot equivalence problem. A topological invariant is
a quantity assigned to a presentation of a knot, that is invariant by isotopy.

An important family of knots are the hyperbolic knots, which are the knots
whose complements admit a complete hyperbolic metric. They are the subjects of
active mathematical research, which motivates the introduction of efficient
algorithmic tools to study their geometric properties, and most notably their
hyperbolic volume. The hyperbolic volume of a hyperbolic knot is a topological
invariant which is powerful at distinguishing between non-equivalent knots, is
non-trivial to compute, and is at the heart of several deep conjectures in
topology [10].

If it exists, the complete hyperbolic metric on a 3-manifold is unique[9], and
may be combinatorially represented by a _complete hyperbolic structure_ (CHS).
In order to compute geometric properties (such as volume) of a hyperbolic
knot, one triangulates the knot complement, and try to assign hyperbolic
shapes to its tetrahedra. If these shapes verify a set of non-linear
constraints called the gluing equations, they form a CHS and encode the
complete hyperbolic metric of the space. A major issue with this approach is
that a solution to the constraints may not exist on all triangulations of a
manifold, even if, as topological object, the manifold can carry a complete
hyperbolic metric. Worst, it is not known whether every hyperbolic 3-manifold
admits a triangulation on which a solution as CHS exists.

In practice, given an input knot, software can construct a triangulation of
the knot complement, and simplify it to have a CHS. But if this fails, the
only implemented practical solution, proposed by SnapPy[3], is to randomly
modify and simplify the triangulation before trying again until a CHS is
found.

| Alternating | Non-Alternating
---|---|---
#crossings | 12 | 13 | 14 | 15 | 16 | 17 | 12 | 13 | 14 | 15 | 16 | 17
% Failure on first try | 0.9 | 1.6 | 2.4 | 3.3 | 4.4 | 5.7 | 0.8 | 0.7 | 1.2 | 1.7 | 2.4 | 3.3
Expected nb of retriang. | 2.9 | 3.9 | 6.2 | 7.1 | 9.5 | 15.9 | 2.0 | 3.1 | 6.1 | 11.8 | 9.7 | 13.5
Table 1: On all $\sim 9.7$ millions prime knots with crossing numbers ranging
from 12 to 17, alternating and non-alternating, we indicate (% Failure on
first try) the percentage of knot complements (after triangulation and
simplification) on which SnapPy fails to compute a CHS on first try. We also
indicate (Highest expected nb of retriang.) the highest, over all knots,
expected number of random re-triangulations necessary for SnapPy to succeed
finding a CHS.  Figure 1: Distribution of the number of randomizations
required to find a CHS for the complement of the knot “17nh_2654001” of the
census, for 10000 tries. The mean is 10.3 randomizations, and the standard
deviation is 8.7. The minimum number of re-triangulations is 1, and the
maximum 92.

Table 1 provides data on the search for a CHS with SnapPy, on the $\sim 9.7$
millions prime knots with crossing numbers up to 17. As observed, SnapPy has a
high rate of success in finding a CHS after a standard triangulation and
simplification of the knot complement. However, this standard construction of
a triangulation fails to admit a CHS on more than $350$ thousand knots in the
census, and the percentage of problematic triangulations needing re-
triangulations tends to increase with the number of crossings. Additionally,
we observe that some knot complements may require in expectation a high number
of random re-triangulations (up to 15.9), and the number of re-triangulations
may itself suffer a high variance, as illustrated in Figure 1.

Checking for the existence of a CHS requires the resolution of the non-linear
gluing equations (with, e.g., Newton optimization method). Reducing the number
of re-triangulations is consequently critical for performance of computation
in knot theory, most notably for knots on which the state-of-the-art methods
implemented in SnapPy require large numbers of re-triangulations, and even
more so when proceeding to very large scale experiments such as the
computation of the knot censuses [1] that are of great use to practitioners,
where “difficult” knots are many.

Contribution: This paper introduces a new heuristic based algorithm to improve
on the random approach for re-triangulating. The method is inspired by Casson
and Rivin’s reformulation of the gluing equations [14, 4]: the gluing
equations are split into a linear part and a non-linear part, and the
resolution reduces to a convex optimization problem on a polytope domain. If
the triangulation does not admit a complete structure, the optimization
problem will converge on the boundary of the polytope and we exploit this
information to introduce new localized moves to modify combinatorially the
triangulation while reusing the partially computed geometry. We introduce
necessary background on triangulations and geometry in Section 2, and the
computation of hyperbolic structures with optimization in Section 3. We
analyze precisely the behavior of the optimization phase on triangulations not
admitting a CHS in Section 4, and introduce a re-triangulation algorithm in
Section 5 guided by the partially computed geometry of the optimization phase.
We illustrate experimentally the interest of the approach in Section 6 and
propose a hybrid method with SnapPy in Section 6.3, that outperforms the
state-of-the-art.

Note that, in this article we focus on complements of hyperbolic knots.
However, the techniques introduced extend to more general hyperbolic
3-manifolds with torus boundaries [4].

##  2 Background

In this article, we focus on knot complements, i.e., non-compact 3-manifolds
obtained by removing a closed regular neighborhood of a knot $K$ in the sphere
$S^{3}$.

###  2.1 Generalized and ideal triangulations

Figure 2: Left: ideal triangulation of the complement of the 8-knot with two
tetrahedra. The ideal vertex is truncated, and the red surface gives the torus
link of the ideal vertex after gluing of the tetrahedra following the edge
identifications. Middle: ideal tetrahedron in the upper half-space model, the
dihedral angles are denoted $\alpha$, $\beta$ and $\gamma$, the complex shape
parameter $z$ is associated to the edge between $0$ and $\infty$. Right:
example of shearing singularity, each triangle corresponds to a ideal
tetrahedron seen from above (gray arrow).

A _generalized triangulation_ $T$ is a collection of $n$ abstract tetrahedra
whose triangular facets are identified (e.g. Figure 2 (left)), or _glued_ , in
pairs. Note that the facets of the same tetrahedron may be glued together, and
generalized triangulations are more general than simplicial complexes. The
link of a vertex in $T$ is the frontier of a closed regular neighborhood, and
is itself a closed triangulated surface embedded in the triangulation. If the
link of a vertex $v$ is a 2-sphere, we call $v$ an _internal_ vertex,
otherwise (e.g., if the link if a torus), we call $v$ an ideal vertex.

A $1$-vertex ideal triangulation is a triangulation with exactly 1 ideal
vertex, and no internal vertices. They represent non-compact 3-manifolds, that
can be recovered from the 1-vertex ideal triangulation by considering their
realization where the vertex has been removed. Every knot complement can be
represented by a 1-vertex ideal triangulation, where the link of the ideal
vertex gives the frontier of a closed regular neighborhood of the knot.
Intuitively, this is a triangulation of the sphere $S^{3}$ where the knot has
been shrunk into a single point, distorting its neighborhood. Such 1-vertex
ideal triangulations of a knot complement can be computed in polynomial time
[6, 5] from a planar drawing of a knot.

Any two ideal triangulations of the same 3-manifold can be connected by a
sequence of _Pachner moves_ [12]. For 1-vertex ideal triangulations, they
consist of the moves 2-3 and 3-2, inverse of each other, pictured in Figure 3.

Figure 3: Illustration of the Pachner moves 2-3 and 3-2.

###  2.2 Combinatorial description of hyperbolic geometry

Certain topological 3-manifolds can be equipped uniformly with a complete
hyperbolic metric, which is unique up to isometry. They are called _hyperbolic
manifolds_. They include the vast and important family of complements of
_hyperbolic knots_.

We use the upper half-space model to represent the hyperbolic space
${\mathbb{H}}^{3}$. This representation corresponds to
$\\{(z,r)|z\in\mathbb{C},r\in\mathbb{R}_{+}^{*}\\}$ with
$\partial\overline{{\mathbb{H}}^{3}}=\mathbb{C}\cup\\{\infty\\}$ (from now on
denoted $\partial{\mathbb{H}}^{3}$) consisting of the bottom plane
$\mathbb{C}\times\\{0\\}$, together with the point at infinity, where
geodesics are arcs of circles orthogonal to $\partial{\mathbb{H}}^{3}$. This
model is _conformal_ , i.e., Euclidean angles in the upper half space give the
values of the angles in the hyperbolic space.

An _ideal hyperbolic tetrahedron_ is the hyperbolic convex hull of four
distinct points of $\partial{\mathbb{H}}^{3}$. These points are the (ideal)
vertices of the tetrahedron.

A vertex on $\partial{\mathbb{H}}^{3}$ is a vertex at infinity, thus it is not
part of the tetrahedron. An example of ideal tetrahedron is shown in Figure 2
(middle). Up to isometry, the geometric shape of an ideal tetrahedron can be
represented by a single complex number:

######  Definition 2.1 (Shape parameter).

Given an ideal hyperbolic tetrahedron, there exists an isometry sending three
of its vertices to $0$, $1$, and $\infty$, and ensuring that the fourth vertex
has positive imaginary part in the complex plane in
$\partial{\mathbb{H}}^{3}=\mathbb{C}\cup\\{\infty\\}$. The coordinate $z$ of
this fourth vertex is the _shape parameter_ of the tetrahedron.

The shape parameter defines and illustrates the shape of an ideal tetrahedron.
It depends on which vertices are sent to $0$, $1$, and $\infty$. Other
permutations of the vertices give other equivalent shape parameters
$z^{\prime}=\frac{z-1}{z}$ and $z^{\prime\prime}=\frac{1}{1-z}$ but the
underlying tetrahedron is the same. The construction is well defined as
isometries of ${\mathbb{H}}^{3}$ are determined by their action on three
vertices of $\partial{\mathbb{H}}^{3}$.

Another way of characterizing the shape of an ideal hyperbolic tetrahedron is
to consider its dihedral angles, i.e., the angle formed by two faces meeting
on a common edge; see Figure 2 (middle). In an ideal tetrahedron, opposite
angles are equal, and the sum of the six dihedral angles is $2\pi$. We denote
in the following the dihedral angles of a tetrahedron by a triplet
$(\alpha,\beta,\gamma)$ with $(\alpha+\beta+\gamma)=\pi$.

##  3 Angle structures and hyperbolic volume

In this section, we introduce notions connected to the computation of complete
hyperbolic structures on triangulations, via optimization methods. The
approach given in this section is another formulation of Thurston’s gluing
equations, which are non-linear in the complex shape parameters mentioned
above, we refer the reader to [4, 14] for more details.

###  3.1 Linear equations and angle structures

Let $T$ be a 1-vertex ideal triangulation of a knot complement $M$ with $n$
tetrahedra. Since opposite edges have the same dihedral angles, all possible
shapes of the tetrahedra can be represented by a vertex in
${\mathbb{R}}^{3n}$. We define an _angle structure_ :

######  Definition 3.1 (Angle structures).

Given an ideal triangulation $T$, an _angle structure_ is a value assignment
to the dihedral angles of the tetrahedra of $T$ such that:

  1. 1.

all the angles are strictly positive;

  2. 2.

the three dihedral angles $(\alpha,\beta,\gamma)$ of a tetrahedron sum to
$\pi$;

  3. 3.

the dihedral angles around each edge of $T$ sum to $2\pi$.

The set of angle structures on $T$ is denoted $\mathcal{A}(T)$.

Conditions 1 and 2 ensure the angles are in $(0,\pi)$ and that the tetrahedra
have the same orientation. Condition 3 is necessary for points on the interior
of edges to have a neighborhood isometric to a hyperbolic ball.

Constraints 1, 2, 3 are linear, hence the set $\mathcal{A}(T)$ of angle
structures of a triangulation is the relative interior of a polytope in
${\mathbb{R}}^{3n}$. It is of dimension $n+|\partial M|$ where $|\partial M|$
is the number of cusps of $M$; in the case of knot complements, there is a
single cusp. This polytope satisfies:

######  Theorem 3.1.

(Casson; see [8]) Let $T$ be an ideal triangulation of $M$, an orientable
3–manifold with toric cusps. If $\mathcal{A}(T)\neq\emptyset$, then $M$ admits
a complete hyperbolic metric.

Furthermore, [4] describes precisely a generating family for the tangent space
of $\mathcal{A}(T)$ that can be computed in polynomial time.

###  3.2 Maximizing the hyperbolic volume

The hyperbolic volume of an ideal hyperbolic tetrahedron with dihedral angles
$(\alpha,\beta,\gamma)\in(0,\pi)^{3}$ is given by the function $\mathrm{vol}$:

$\mathrm{vol}(\alpha,\beta,\gamma)=\text{L}(\alpha)+\text{L}(\beta)+\text{L}(\gamma)$

where L is the Lobachevsky function $\text{L}(x)=-\int_{0}^{x}\log|2\sin
t|\,\mathrm{d}t$. The volume functional can be extended to the whole polytope
$\mathcal{A}(T)$ by summing the volumes of the hyperbolic tetrahedra. The
following result is due independently to Casson and Rivin.

######  Theorem 3.2 (Casson, Rivin[13]).

Let $T$ be an ideal triangulation with $n$ tetrahedra of $M$, an orientable
3–manifold with boundary consisting of tori. Then a point
$p\in\mathcal{A}(T)\subset{\mathbb{R}}^{3n}$ corresponds to a complete
hyperbolic metric on the interior of $M$ if and only if $p$ is a critical
point of the function $\mathrm{vol}$.

Additionally, the volume functional is concave on $\mathcal{A}(T)$ and the
maximum can be computed via convex optimization methods.

We can finally define CHS that represent combinatorially complete hyperbolic
metrics:

######  Definition 3.2 (Complete Hyperbolic Structure).

A Complete Hyperbolic Structure (CHS) is a triangulation $T$ equipped with an
angle structure corresponding to the maximum of $\mathrm{vol}$ over
$\mathcal{A}(T)$.

######  Remark 3.1.

Angle structures are not CHS as they do not prevent shearing singularities,
see Figure 2 (right), and consequently may represent _non-complete_ hyperbolic
metrics.

##  4 Behavior of the optimization

Exploiting the formalism of the previous section, one can design an algorithm
[4] to find and compute a CHS by first finding a point in the polytope
$\mathcal{A}(T)$, then computing a basis of the tangent space of
$\mathcal{A}(T)$, and then maximizing the (concave) hyperbolic volume
functional on this subspace. If $\mathcal{A}(T)\neq\emptyset$ and the
procedure finds the point maximizing the volume in the inside of the polytope,
then the triangulation admits a complete hyperbolic structure represented by
this point. Finally, if $\mathcal{A}(T)\neq\emptyset$ but the maximum of the
volume functional is on the boundary, then one needs to re-triangulate the
manifold in order to search for a triangulation admitting a CHS.

In this section, we study the outcome of the hyperbolic volume maximization on
the space of angle structures. The result is given by the following lemma,
which also appeared independently in [11]:

######  Lemma 4.1.

Let $T$ be an ideal triangulation of $M$, a non-compact orientable 3–manifold
with toric cusps. Let $p$ be the point maximizing $\mathrm{vol}$ over
$\overline{\mathcal{A}(T)}$ (the topological closure of $\mathcal{A}(T)$),
then at $p$, if a tetrahedron has an angle equal to $0$, then all other angles
of the tetrahedron are in $\\{0,\pi\\}$.

In other words, either the maximization succeeds in $\mathcal{A}(T)$, or there
is at least one tetrahedron with angles $(0,0,\pi)$ in $p$. Such tetrahedron
is called _flat_. It is not possible to have a tetrahedron with angles
$(0,a,b)$, $(a,b)\in(0,\pi)$ in $p$.

###### Proof.

Let $\vec{p}\in\mathcal{A}(T)$ be an angle structure over $T$. Let
$\vec{w}\in\mathbb{R}^{3n}$ be a vector tangent to $\mathcal{A}(T)$, then:

$\frac{\partial\mathrm{vol}(\vec{p})}{\partial\vec{w}}=\sum_{i=1}^{3n}-w_{i}\log\sin(p_{i}).$

Note that because $\vec{w}$ is tangent to $\mathcal{A}(T)$, its restriction to
a single tetrahedron is three elements whose sum is equal to zero. Let us
assume the maximum of the volume is reached on some point $\hat{p}$. Let us
assume a tetrahedron $t$ has angles $(0,a,b)$ at $\hat{p}$, $0<a\leq b$, let
$\hat{w}$ be a vector tangent to $\mathcal{A}(T)$ pointing towards the
interior of the polytope with respect to $\hat{p}$ (and thus $w_{1}>0$) and
$\hat{w}_{t}=(w_{1},w_{2},w_{3})$ be its restriction to $t$. Then:

$\lim_{\epsilon\rightarrow
0^{+}}\frac{\partial\mathrm{vol}(t+\epsilon\hat{w}_{t})}{\partial\hat{w}_{t}}=-w_{2}\log\sin(a)-w_{3}\log\sin(b)+\lim_{\epsilon\rightarrow
0^{+}}-w_{1}\log\sin(\epsilon w_{1})=+\infty.$

Now, with the same notations, let us assume $t$ has angles $(0,0,\pi)$ at
$\hat{p}$, then using $\log(\sin(x))=\log(x)+O_{x\rightarrow 0}(x^{2})$:

$\displaystyle\lim_{\epsilon\rightarrow 0^{+}}$
$\displaystyle\frac{\partial\mathrm{vol}(t+\epsilon\hat{w}_{t})}{\partial\hat{w}_{t}}$
$\displaystyle=-w_{1}\log\sin(\epsilon w_{1})-w_{2}\log\sin(\epsilon
w_{2})+(w_{1}+w_{2})\log\sin(\pi-\epsilon(w_{1}+w_{2}))$
$\displaystyle=-w_{1}\log(\epsilon w_{1})-w_{2}\log(\epsilon
w_{2})+(w_{1}+w_{2})\log(\epsilon(w_{1}+w_{2}))+O_{\epsilon\rightarrow
0}(\epsilon^{2})$
$\displaystyle=-w_{1}\log(w_{1})-w_{2}\log(w_{2})+(w_{1}+w_{2})\log(w_{1}+w_{2})+O_{\epsilon\rightarrow
0}(\epsilon^{2})$ $\displaystyle>-\infty$

Since for all tetrahedra with non zero angles or angles of the form
$(0,0,\pi)$, the derivative of the volume is bounded, the existence of a
tetrahedra with exactly one angle set to zero contradicts the maximality of
the volume at that point. Indeed, when getting away from this point, there
exists a small neighborhood where the derivative of the volume can be
arbitrarily high.

∎

We introduce in the next section an algorithm that performs localized
combinatorial modifications on a triangulation equipped with an angle
structure, in order to get rid of flat tetrahedra.

##  5 Localized combinatorial modifications of triangulations

According to Lemma 4.1, the volume maximization either leads to a solution to
the gluing equations, or to flat tetrahedra. In this section, we discuss a
method to get rid of flat tetrahedra by combinatorial modifications of the
triangulation, while attempting to maintain the value of the volume
functional. This would lead to a triangulation admitting an angle structure
with larger volume than the previous one, allowing to resume the maximization.
In order to maintain the value of the volume functional, we introduce
geometric Pachner moves.

###  5.1 Geometric Pachner moves

We perform Pachner moves on the triangulation that preserves the partial
geometric data computed. More precisely, we define:

######  Definition 5.1 (Geometric Pachner move).

A _geometric Pachner move_ in a triangulation $T$ with angle structure is a
Pachner move in $T$ such that the resulting triangulation admits an angle
structure with identical dihedral angles for the tetrahedra not involved in
the move.

To check if a geometric Pachner move can be done, it must be valid
combinatorially, and it should be possible to assign dihedral angles to the
new tetrahedra without altering the rest of the angle structure. The
combinatorial conditions are simple [12]: for the 2-3 moves (resp. 3-2 move)
in Figure 4, the tetrahedra sharing the common triangle $BDC$ (resp. common
edge $AE$) must be distinct. The conditions on the angle structures are given
by the following lemma:

######  Lemma 5.1.

Given a triangulation with angle structure $T$, if the combinatorial
conditions given above are satisfied:

  * •

a 3-2 move is always geometric;

  * •

a 2-3 move is geometric if and only if the sum of the two dihedral angles
around each edge of the common face is smaller than $\pi$.

Furthermore the new angles can be computed using formulas.

Figure 4: Pachner moves 2-3 and 3-2 with vertices corresponding to the
notations of the lemmas of this section.

The idea to prove the lemma is to study the system of equations given by the
assumption that the angle structure is not altered outside the move, and we
break the lemma into three: Lemma 5.2 concerns the 3-2 move, Lemma 5.3
concerns the 2-3, and Lemma 5.5 concerns the values of the angles.

With the notations of Figure 4, let $t_{1}=ABCD$ and $t_{2}=BCDE$ be two
distinct tetrahedra glued along their common face, and $t_{3}=ABCE$,
$t_{4}=ACDE$ and $t_{5}=ADBE$ three distinct tetrahedra glued around their
common edge. Assuming there is an underlying angle structure, let us denote
$t(e)$ the value of the dihedral angle around the edge $e$ in the tetrahedron
$t$. Passing from $(t_{1},t_{2})$ to $(t_{3},t_{4},t_{5})$, i.e. doing a 3-2
move, or the opposite without altering the rest of the angle structure leads
to the following equations:

$t_{1}(BC)+t_{2}(BC)=t_{3}(BC),\ t_{1}(CD)+t_{2}(CD)=t_{4}(CD),\
t_{1}(DB)+t_{2}(DB)=t_{5}(DB)$ (1)

$t_{1}(AB)=t_{3}(AB)+t_{5}(AB),\ t_{1}(AC)=t_{3}(AC)+t_{4}(AC),\
t_{1}(AD)=t_{4}(AD)+t_{5}(AD)$ (2)

$t_{2}(EB)=t_{3}(EB)+t_{5}(EB),\ t_{2}(EC)=t_{3}(EC)+t_{4}(EC),\
t_{2}(ED)=t_{4}(ED)+t_{5}(ED)$ (3)

The constants, the angles that are fixed in the angle structure, and the
unknowns, the angles we want to compute, depends on which move is studied. By
equality of the dihedral angles, $t(UV)=t(XY)$ for $UVXY$ vertices of $t$.
This equality is enforced by the framework of the angle structures: for the
constraints system, $t(UV)$ and $t(XY)$ both represent the same variables.

####  5.1.1 The 3-2 move

######  Lemma 5.2.

Let $T$ be a triangulation with angle structure, let $t_{3}=ABCE$,
$t_{4}=ACDE$ and $t_{5}=ADBE$ be three distinct tetrahedra glued around an
edge. A 3-2 geometric Pachner move can always be performed at $AE$.

###### Proof.

The 3-2 move is combinatorially valid as the resulting tetrahedra and their
gluings are well defined. Let us have $t_{1}=ABCD$ and $t_{2}=BCDE$.

The system composed of Equations 1, 2 and 3 is linear with 6 unknowns: the
dihedral angles of $t_{1}$ and $t_{2}$. We need to check if this system admits
a solution where the angles are in $(0,\pi)$. Note that as $t_{3}(AE)$,
$t_{4}(AE)$, and $t_{5}(AE)$ are the angles around the edge $AE$, their sum is
$2\pi$.

$t_{1}(AB)=$ | $t_{3}(AB)+t_{5}(AB)$ | $=$ | $2\pi-(t_{3}(AE)+t_{3}(AC)+t_{5}(AE)+t_{5}(AD))$
---|---|---|---
$<$ | $2\pi-(t_{3}(AE)+t_{5}(AE))$ | $=$ | $t_{4}(AE)$
$<$ | $\pi$ |  |

$t_{1}(AB)=t_{3}(AB)+t_{5}(AB)>0$
---

Thus, if the system admits a solution, the dihedral angles of $t_{1}$ and
$t_{2}$ are in $(0,\pi)$ by symmetry.

Now let us check if all the equations can be satisfied. From Equation 2 and 3:

$t_{1}(AC)+t_{2}(EC)=t_{3}(AC)+t_{4}(AC)+t_{3}(EC)+t_{4}(EC)=2\pi-(t_{3}(AE)+t_{4}(AE))=t_{5}(AE)$

Thus the equations of Equation 1 are consequences of the other two. The system
is linear with 6 independent equations and 6 unknowns, it admits a solution.

Finally we have:

$t_{1}(AB)+t_{1}(AC)+t_{1}(AD)=$ | $t_{3}(AB)+t_{5}(AB)+t_{3}(AC)+t_{4}(AC)+t_{4}(AD)+t_{5}(AD)$
---|---
$=$ | $3\pi-(t_{3}(AE)+t_{4}(AE)+t_{5}(AE))$
$=$ | $\pi$

Thus the dihedral angles of $t_{1}$ sum to $\pi$, the same goes for $t_{2}$ by
symmetry. And the move is always valid as long as we started with an angle
structure. ∎

####  5.1.2 The 2-3 move

######  Lemma 5.3.

Let $T$ be a triangulation with angle structure, let $t_{1}=ABCD$ and
$t_{2}=BCDE$ be two distinct tetrahedra glued along a face. A 2-3 geometric
Pachner move can be performed at $BCD$ if and only if $\forall
e\in\\{BC,CD,DB\\},\ t_{1}(e)+t_{2}(e)<\pi$.

###### Proof.

The 2-3 move is combinatorially valid as the resulting tetrahedra and their
gluings are well defined. Let us have $t_{3}=ABCE$, $t_{4}=ACDE$ and
$t_{5}=ADBE$.

From Equations 1 we have that if $\exists e\in\\{BC,CD,DB\\},\
t_{1}(e)+t_{2}(e)\geq\pi$ then one of the resulting tetrahedra will have a
dihedral angle larger than $\pi$, preventing the validity of the move. Let us
assume it is not the case.

The system composed of Equations 1, 2 and 3 is linear with 9 unknowns: the
dihedral angles of $t_{3}$, $t_{4}$ and $t_{5}$. From Equations 1, the values
of $t_{3}(AE)$, $t_{4}(AE)$, and $t_{5}(AE)$ are fixed and we have:

$t_{3}(AE)+t_{4}(AE)+t_{5}(AE)\\\
=t_{1}(BC)+t_{2}(BC)+t_{1}(CD)+t_{2}(CD)+t_{1}(DB)+t_{2}(DB)=2\pi$

Thus the constraint around the edge $AE$ is satisfied.

Summing the equations of Equation 2 on the one hand, and the equations of
Equation 3 on the other hand both gives:

$\pi=t_{3}(AB)+t_{5}(AB)+t_{3}(AC)+t_{4}(AC)+t_{4}(AD)+t_{5}(AD)$

We also have that the sum of the dihedral angles of $t_{3}$ is $\pi$ (and by
symmetry it is true for $t_{4}$ and $t_{5}$):

$\displaystyle t_{3}(AB)+t_{3}(AC)+t_{3}(AE)$
$\displaystyle=t_{1}(BC)+t_{2}(BC)+t_{1}(AB)+t_{1}(AC)-t_{5}(AB)-t_{4}(AC)$
$\displaystyle=t_{1}(BC)+t_{1}(AB)+t_{1}(AC)+t_{2}(BC)-t_{2}(ED)$
$\displaystyle=\pi$

Now, let us note that fixing any unknown fixes the values of the others, thus
the system has exactly one degree of freedom: if
$(t_{3}(AB),t_{3}(AC),t_{4}(AC),t_{4}(AD),t_{5}(AD),t_{5}(AB))$ is a solution
to the system, then for any $\lambda\in{\mathbb{R}}$,
$(t_{3}(AB)+\lambda,t_{3}(AC)-\lambda,t_{4}(AC)+\lambda,t_{4}(AD)-\lambda,t_{5}(AD)+\lambda,t_{5}(AB)-\lambda)$
is also a solution. This is discussed in Remark 5.1. We denote by
$\Lambda_{+}$ the set of dihedral angle for which $\lambda$ is added, and by
$\Lambda_{-}$ its complement.

The last step of the proof is to show that there exists a choice of $\lambda$
such that all the dihedral angles are positive, given that the dihedral angles
of the tetrahedra sum to $\pi$ is enough to have them in $(0,\pi)$.

By Lemma 5.4, there is always a choice of dihedral angle such that fixing it
to 0 ensures that the elements of $\Lambda_{+}$ are positive and those of
$\Lambda_{-}$ are non negative (or $\Lambda_{-}$ and $\Lambda_{+}$
respectively, depending on whether the selected angle is in $\Lambda_{+}$ or
$\Lambda_{-}$). If the elements of $\Lambda_{+}$ are positive, then picking a
sufficiently small negative $\lambda$ to modify the dihedral angles will
maintain the elements of $\Lambda_{+}$ positive, and will turn positive those
of $\Lambda_{-}$, in the other case, in the same way a small positive
$\lambda$ will ensure the positivity of all the angles.

Thus there exists an angle structure corresponding to and coherent with the
result of the Pachner move.

∎

######  Lemma 5.4.

With the notations and the assumption of the proof of Lemma 5.3, there is
always a choice of dihedral angle $x$ in
$\\{t_{3}(AB),t_{3}(AC),t_{4}(AC),t_{4}(AD),t_{5}(AB),t_{5}(AD)\\}$ such that
fixing $x=0$ leads to the elements of $\Lambda_{+}$ being positive and those
of $\Lambda_{-}$ being non negative (or $\Lambda_{-}$ and $\Lambda_{+}$
respectively).

###### Proof.

Let us assume the Lemma is false.

If $t_{3}(AC)$ is set to $0$, then

$t_{3}(AB)=\pi-t_{3}(AE)>0,$ $t_{5}(EB)=t_{2}(EB)>0$ $t_{4}(ED)=t_{1}(AC)>0.$

Furthermore

$t_{4}(AD)=\pi-(t_{4}(AC)+t_{4}(AE))=t_{1}(AD)-t_{2}(EB)$

and

$t_{5}(AB)=t_{2}(ED)-t_{1}(AC).$

Note that the elements of $\Lambda_{+}$ are positive and the elements of
$\Lambda_{-}$ are non-negative if $t_{4}(AD)$ and $t_{5}(AB)$ are non-
negative. Let us assume $t_{1}(AD)-t_{2}(EB)<0$.

If $t_{4}(AD)$ is set to $0$, like previously the elements of $\Lambda_{+}$
are positive, and

$t_{3}(AC)=t_{2}(EB)-t_{1}(AD)>0$ $t_{5}(AB)=t_{1}(AB)-t_{2}(EC).$

If the lemma is false, then $t_{1}(AB)-t_{2}(EC)<0$.

If $t_{5}(AB)$ is set to $0$, like previously the elements of $\Lambda_{+}$
are positive, and

$t_{3}(AC)=t_{1}(AC)-t_{2}(ED)$ $t_{4}(AD)=t_{2}(EC)-t_{1}(AB)<0$

. If the lemma is false, then $t_{1}(AC)-t_{2}(ED)<0$.

The previous deductions give

$t_{1}(AB)+t_{1}(AC)+t_{1}(AD)<t_{2}(EB)+t_{2}(EC)+t_{2}(ED),$

which is absurd as these two sums are equal to $\pi$. The assumption that
$t_{1}(AD)-t_{2}(EB)<0$ must be false, and by symmetry $t_{2}(ED)-t_{1}(AC)<0$
will be false too.

This is absurd and the Lemma is true.

∎

######  Remark 5.1 (Non-conservation of the volume).

In Lemma 5.3, several angle structures are possible after the 2-3 move as it
can be seen in the poof of Lemma 5.4. Conversely, several angle structures
lead to the same result for the 3-2 move. This is because angle structures are
blind to shearing singularities. If such a singularity exists, a 3-2 move can
still be performed while maintaining the whole angle structure, but it will
decrease the hyperbolic volume of the structure. It is due to the fact that
the singularity, which maximized the volume, is fixed in the process.

####  5.1.3 Computing the angles

######  Lemma 5.5.

Let $T$ be a triangulation with angle structure, if a geometric Pachner move
can be performed, then the new angle structure can be computed in a constant
number of arithmetical operations.

###### Proof.

3-2 move.  From the proof of Lemma 5.2, the new angles can be obtained by
solving a system of equations with six equations.

2-3 move.  First, we remind that it is possible to convert angles to edge
parameters et vice versa [4]. We still consider the notations of Figure 4, but
we consider the situation in the upper halfspace with $B$ sent to $0$, $E$
sent to $1$, and $D$ sent to $\infty$ (Figure 5). The position of $C$ depends
on the edge parameter of $BD$ in $BCDE$, and the position of $A$ depends on
the edge parameter of $BD$ in $ABCD$. Then the edge parameter of $BD$ in
$ABDE$ is just the product of the two edge parameter. This construction is
valid by Lemma 5.3, as the sum of the arguments of both shape parameters is
smaller than $\pi$. We can use the same procedure to find an edge parameter
for $ABCE$ and $ACDE$.

Note that this construction leaves no shearing singularities around $AE$. ∎

Figure 5: Representation of the proof of Lemma 5.5, seen from the point $D$.
The 2-3 move will delete the line $BC$ in favor of $AE$.

###  5.2 Getting rid of flat tetrahedra

Let $T$ be a triangulation with an angle structure admitting a flat
tetrahedron $t$, and $(e,e^{\prime})$ the associated edges with dihedral angle
$\pi$. By Lemma 5.1, it is not possible to get rid of a tetrahedron with a
geometric 2-3 move. Indeed, if an edge of $t$ has a dihedral angle equal to
$\pi$, the second tetrahedron concerned by the 2-3 move must have a dihedral
angle equal to 0, and consequently a 2-3 move will produce a flat tetrahedron.

In order to get rid of $t$, our strategy is to turn either $e$ or $e^{\prime}$
into an edge of degree three, i.e. to put it at the center of three tetrahedra
on which a 3-2 move can be performed (central edge of Figure 4, right). Note
that edges of degree 2 prevent the existence of an angle structure as they
force the two tetrahedra sharing the edge to be flat. As a consequence, we
assume these configurations are removed, which can be done by SnapPy’s
simplification for instance, see Section 6. Consequently, all edges have
degree at least 3, and our strategy focuses on _reducing_ the incidence degree
of angle $\pi$ edges.

Figure 6: A sequence of moves getting rid of a flat tetrahedron. The bold
vertical edge is contained in 4 tetrahedra, three of which are represented
behind the edge. The fourth tetrahedron is implicit, situated at the front,
and flat. First, create the red edge of the first drawing with a 2-3 move.
Then delete the red edge of the second drawing with a 3-2 move.

To reduce the degree of an edge $e$, the strategy is to perform 2-3 moves on
the tetrahedra containing $e$, such that each move reduces the degree of the
edge by one. Either for topological (tetrahedron glued to itself) or
geometrical reasons, these moves will not always be possible, and the order in
which they are done matters. This is described in Figure 6.

######  Remark 5.2.

Doing a 2-3 move does not always reduce the degree of the edge, e.g., when a
tetrahedron is represented several times around an edge. However, doing a move
that does not decrease the degree of the edge may delete the multiple
occurrences of a tetrahedron around the edge and allow to continue with the
simplification.

######  Remark 5.3.

When a 2-3 move is performed around $e$, the value of the dihedral angle of
the new tetrahedron around $e$ is equal to the sum of the previous two
dihedral angles around $e$. Since the sum of the dihedral angles of a
tetrahedron is equal to $\pi$, this means successfully reducing the degree of
$e$ may unlock previously forbidden moves.

Recursive moves. If no geometric Pachner move is possible on an edge $e$, we
attempt to reduce the degree of a nearby edge $e_{f}$, i.e., an edge for which
there exists a tetrahedron containing $e$ and $e_{f}$; see Figure 7 where
$e_{f}$ is represented by the red dot, as it is seen from above. More
precisely, let $f$ be a triangle containing $e$ such that the associated 2-3
move is forbidden for geometric reasons, then by Lemma 5.1 there is an edge
$e_{f}$ of $f$ for which the associated dihedral angle is at least $\pi$.

In consequence, reducing the degree of $e_{f}$ to three without modifying the
tetrahedra containing $f$ and then performing a 3-2 move at $e_{f}$ will
reduce the degree of $e$ by one; see Figure 7.

Note that in case a geometric Pachner move is not possible on any pair of
tetrahedra containing $e_{f}$ either, we can call the procedure recursively in
a neighborhood of $e_{f}$, using the argument above to pass from $e$ to
$e_{f}$.

Figure 7: Sequences of moves to reduce the degree of an edge when no 2-3 move
is available, seen from above (triangles represent sections of tetrahedra,
edge section of triangles and dots section of edges). In all drawings, the
gray rectangle represents the flat tetrahedron, $e$ (light blue dot) is the
edge of the flat tetrahedron of which we want to reduce the degree. Assuming
no 2-3 move is available, we pick another edge $e_{f}$ (red dot), if a 3-2
move can be performed on $e_{f}$, we are in case (B) and the 3-2 move reduces
the degree of $e$ (case (C)). Otherwise, we are in case (A) we need to create
the red dashed tetrahedron by decreasing the degree of $e_{f}$, this can be
done by calling recursively our procedure on $e_{f}$.

######  Remark 5.4.

The choice of $e_{f}$ in $f$ is unique: it is not possible to have two
dihedral angles larger than $\pi$; selecting a move to reduce the degree of
$e$ boils down to selecting a triangle containing $e$.

Procedure summary. The procedure to get rid of a flat tetrahedron is a tree-
like _backtracking search_ , each branch corresponding to a choice of 2-3 or
3-2 move to perform. The moves used are geometric to preserve the volume
maximization advancement. A flat tetrahedron can only be deleted with a 3-2
move: the first step is to _reduce the degree_ of one of its edges to three.
To reduce the degree of an edge, the solution is to perform 2-3 moves on two
tetrahedra containing this edge. If a 2-3 move is not possible for geometric
reasons, we are in the situation of Figure 7, and the procedure to _reduce the
degree_ can be called _recursively_ on a neighboring edge. Performing a 3-2
move on this latter edge resulting in a decrease of the degree of the initial
one. The procedure ends when the initial edge has degree three and the flat
tetrahedron can be removed or when no move can be applied.

###  5.3 Implementation details

The implementation faces several practical challenges.

Breaking infinite loops. As such, the algorithm may loop infinitely on some
instance, as 2-3 and 3-2 moves may reverse themselves. To counteract this
phenomenon and avoid redundant modifications, a solution is to store at each
modification the _isomorphism signature_ of the triangulation [2],
characterizing uniquely the isomorphism type of the triangulation. This allows
us to recognize already processed triangulations and break branches of the
backtracking algorithms. Note however that this is costly compared to the
Pachner moves, and it makes the procedure no longer local.

Selecting the edge $e$. When attempting to remove a flat tetrahedron, one
needs to choose between the two $\pi$-angled edges $e$ and $e^{\prime}$ to
reduce. In our implementation, we consider both edges, and select the one that
has the larger smallest angle in its link. This performs better than a random
choice in our experiments.

Pruning the backtrack search. Some triangulations may have edges of large
degree (more than 10), which produces wide search trees. Additionally, the
recursive calls to the degree reduction procedure may induce trees of large
depth. Experimentally, the better strategy consists of exploring exhaustively
the first levels of the tree. We set the width of exploration to $8$ and
explore the first $2$ levels of the search tree.

Sequencing of Pachner moves. Different sequencings of Pachner moves lead to
different triangulations. In practice, we favor 2-3 moves over recursive ones,
which performs best. Additionally, among the 2-3 moves performed to get rid of
a flat tetrahedron containing edge $e$ of angle $\pi$, we prioritize moves
that eliminate tetrahedra for which $e$ has a small dihedral angle. Among the
recursive ones, we used the same method as the one choosing between $e$ and
$e^{\prime}$, indeed choosing a recursive move boils down to choosing an edge
to minimize its degree.

##  6 Experiments

In this section, we study the experimental performance of our approach to find
a triangulation with a CHS, and compare its behavior with the software SnapPy
[3]. SnapPy is the state-of-the-art software to study the geometric properties
of knots and 3-manifolds, and is widely used in the low-dimensional topology
community. Following this analysis, we propose and study a hybrid method with
practical interest.

SnapPy. SnapPy’s method is based on an implementation from Weeks, its uses a
random re-triangulation followed by a simplification. The first step is
constituted of $4n$ random 2-3 moves, where $n$ is the number of tetrahedra in
the input triangulation. The simplification performs non-deterministic
modifications to decrease the number of tetrahedra in the triangulation.
Notably it removes some configurations that prevent angle structures from
existing. The verification of the existence of a CHS is based on a Newton’s
method to solve the gluing equations [17].

Data set. We apply our algorithms to the census of prime hyperbolic knots with
up to 19 crossings. This census has been constituted by the efforts of many
researchers in the field, and recently completed with the exhaustive
enumeration of all knots with crossing number smaller than 20111The census is
available at https://regina-normal.github.io/data.html [1].

For each knot, given by a knot diagram, we compute a triangulation of the knot
complement using Regina [2], and simplify it with SnapPy. We then keep the
triangulations admitting an angle structure but not a CHS, and hence requiring
re-triangulation. All the others, without angle structure or admitting a CHS,
were discarded. They constitute the _Failure on first try_ data of Table 1.

The knots are grouped by crossing numbers (from 14 to 19), and on whether they
are alternating or not. In Section 6.2, our experiments are run on the first
2500 knots of each of the 12 groups having an angle structure but no CHS, and
in Section 6.3 we use the first 10.000 knots of the same data sets.

###  6.1 Our program and current state of the code

Our implementation is done in Python, to be able to interface with Regina and
SnapPy (and not SnaPpea). The idea behind using Python was also to enable
shorter development times. The basic tools to handle triangulations where
rewritten to have more control on the behavior of the simulation and we rely
on external libraries only for the isomorphism signatures, some linear
programming to find a point to begin the optimization, and to perform the
optimization (we use the SLSQP [7] optimization method from SciPy [15]).

However this choice of language leads to costly computation times, decision
was taken to not render public this version of the code and switch to C++,
using another library to manipulate the triangulations.

The version of the program presented in Section 6.3222Available at
https://github.com/orouille/SnapPy, which outperforms the state of the art in
terms of performance, is a simpler C version directly implemented in SnapPy’s
kernel. We are in discussion with the developers to integrate this code in the
library SnapPy.

###  6.2 Success rate and combinatorial performance

In terms of Pachner moves. Figure 8 (left) represents the rate of knots on
which our algorithm succeeds to find a triangulation with a CHS, for all the
12 groups of knots. It represents the success rate as a function of the number
of Pachner moves, the higher success rate groups are the alternating knots.
With Figure 8 (right), which represents this success rate for optimization
phases, we see that few Pachner moves are done on average between the phases,
hinting that the convex optimization will likely be the time bottleneck.

While there is a significant difference between the efficiency on alternating
and non-alternating knots, all the curves have the same behavior: the first
few Pachner moves are very effective. An important point is that, compared to
the number of moves done by SnapPy, at least four times the number of
tetrahedra and then roughly the same number for the simplification times the
number of re-triangulation, our method uses a much lower number of Pachner
moves on average (see Figure 8, left). This is of _combinatorial interest_ ,
as doing a small number of moves can mitigate the impact of the procedure on
the properties of the triangulation, such as keeping a small treewidth for
instance.

After the first few steps, the growth of the success rate slows down
drastically. An interpretation of this phenomenon is that the tetrahedra
created by the 2-3 moves tend to to be more flat than the original ones. This
leads to an increase in the number of required re-triangulations and issues
with floating point arithmetic.

Figure 8: On the sample of triangulations, success rate in finding a CHS for different crossing numbers and alternabilities. The limit rate of success for alternating knots is 0.98, the limit rate for non-alternating knots is 0.87. Left: success rate against the number of Pachner moves. Right: success rate against the number of optimization phases. | # of optim phases |  | # of moves
---|---|---|---
| 1 | 2 | 3 | 4 | 5 |  | 2 | 3 | 4 | 10 | 30
Efficiency alter | 0.73 | 0.89 | 0.95 | 0.97 | 0.97 |  | 0.36 | 0.63 | 0.72 | 0.92 | 0.97
Efficiency non-alter | 0.53 | 0.67 | 0.74 | 0.78 | 0.81 |  | 0.30 | 0.47 | 0.52 | 0.70 | 0.81
Figure 9: Average of success rates for the alternating and non-alternating
groups of knots of Figure 8.

In terms of _re-triangulation_. In order to compare the behaviors of our
method and of the state of the art: we want to know if the two methods
struggle on the same triangulations. Both methods using sequences of re-
triangulations, we use the length of these sequences as metrics.

We compare in Figure 10 the average number of randomized re-triangulations
required by SnapPy to find a triangulation with a CHS, and the (deterministic)
number of re-triangulations of our method. The majority of the CHS are found
within few re-triangulations on average: both methods manage to find 78% of
them within two steps.

However, it appears that for the manifolds requiring a large number of re-
triangulations with either method, the other one will perform well on it: the
performance of the methods are substantially orthogonal for the difficult
cases. This phenomenon highlights the fact that the methods are
_complementary_.

Figure 10: Number of re-triangulations required by our method over the number
of re-triangulations required by SnapPy on the 12 groups of knots, gathered in
bins (the color indicates the log of the number of points inside). Since our
method is deterministic, the points have discrete ordinates. The diagonal
$x=y$ appears in red.

Limitations. It is to be noted that our method suffers from several phenomena:
there are manifolds on which our method fails, the reasons to this being
convergence problems because of the floating point arithmetic, the creation of
flatter tetrahedra or the volume modification mentioned in Remark 5.1.
Furthermore some triangulations have several very high degree edges and
looking for the correct sequence of geometric Pachner moves can be costly. In
terms of time performance, the constraint convex optimization is an important
bottleneck, as well as the computation and book-keeping of a large number of
isomorphism signatures.

However, in light of the previous results, we design a hybrid algorithm,
mixing SnapPy’s pipeline together with the heuristics based on localized
geometric Pachner moves, in order to outperform both methods in terms of
experimental timings.

###  6.3 Hybrid algorithm and time performance

As is, SnapPy does not produce angle structures as defined in Section 3, as it
constructs _negative tetrahedra_ , i.e., tetrahedra whose shape has a negative
imaginary part. Our hybrid method consists of calling SnapPy to randomize and
simplify the triangulation, and, in case the triangulation admits only few
negative tetrahedra (less than four in our experiments), call the resolution
of flat tetrahedra with geometric Pachner moves on these negative tetrahedra,
as introduced in Section 5. Loop until a CHS is found.

Our implementation is done in C directly in SnapPy’s kernel and consists of a
simplified version of Section 5’s method. For each negative tetrahedron, we
try to reduce the degree of one of its edges down to 3 in turn using only
geometric 2-3 moves, we try all the edges, and we do not use recursive moves
or heuristics to perform the geometric moves.

The results of this method are summarized in Figure 11 (Figure 12 is given to
help read Figure 11), where the time to compute the CHS are compared for the
new hybrid method and the usual SnapPy pipeline, on the first 10.000
triangulations of the datasets of Section 6. The hybrid method using localized
geometric Pachner moves shows, at worst, similar time performance compared to
SnapPy, and performs significantly better overall. On average, the hybrid
method is 20% faster, however, it is to be noted that, on most cases, few re-
triangulations are required in which case the hybrid method and SnapPy show
naturally similar time performance. More interestingly, on cases requiring
many more re-triangulations in the SnapPy pipeline, the hybrid methods
performs much better than SnapPy, with running times up to 18 times faster.

As indicator of the global behavior of both algorithms, we indicate the linear
regression of the data points (black dashed line, slope $=0.4$) in Figure 11,
to illustrate that the dense region is not concentrated on the $x=y$ diagonal.
Furthermore, this line suggests that if we are not in the dense region where
almost all CHS are found instantly, then the gain is not 20%, but 60%.

Figure 11: Time comparison of the hybrid method and SnapPy in seconds for the
10.000 first triangulations of each crossing number of the data set of the
previous section, gathered in bins (the color indicates the log of the number
of points inside). Line $x=y$ is plain and red, a linear regression through
the data set is dashed and black. For readability, four outlayers favoring the
hybrid method are omitted.  Figure 12: Distribution, in logarithmic scale, of
the time required by SnapPy over the time required by the hybrid method.

##  7 Conclusion

The problematic of this work is the computation of complete hyperbolic
structures on triangulated 3-manifolds. This can be done on computers by
solving Thurston’s gluing equations. The state-of-the-art method, implemented
by SnapPy, consists in trying to simplify the triangulation before using
Newton’s method to solve the equations. If this does not succeed, then the
only strategy is to find another triangulation of the manifold and to start
over. This work well in most cases, but can require a large number of re-
triangulations in some cases. The key difficulty is to find a triangulation
where the equations admits a solution

We have introduced a new method based on results form Casson and Rivin stating
that maximizing the hyperbolic volume yields CHS. When it is not possible to
find a solution to the equations on a triangulation, we use the result of the
volume maximization to propose another one.

Our algorithm is complementary to SnapPy, and in particular performs much
better on almost all cases we have studied where SnapPy struggled to find a
solution to the CHS problem. While our current implementation of the whole
method is not competitive due to being implemented in Python, we introduced a
hybrid method that improves the state-of-the-art in most cases, specifically
in difficult ones, without sacrificing SnapPy’s efficiency when it performs
well.

Our method proposes an improvement in term of computation time over the state-
of-the-art, this is very important when running large scale experiments or
when the manifolds studied are large. But our method has another interest:
since it tries to avoid performing a large number of re-triangulations, we
tend to preserve the triangulation whenever possible.

The method and its implementation can be improved: when used alone, it can
simply fail, being able to continue with the information we gained would be a
great plus. Also, we plan on transferring the full method in C++, to check its
actual performances when Python is not used, this will rise questions about
the convex programming and arithmetic precision.

The method opens the door for new questions concerning its behavior. A group
of questions arises when looking at the volume maximization like a flow: What
is the impact of the volume maximization on the triangulation? Which
tetrahedra shrink first? Are there combinatorial configurations blocking our
method? These questions are close to the study of the _combinatorial Ricci
flow_ [18], which is dual to our approach, meets similar problems and has
known some development in the recent years. Closely linked to the flow, there
are questions on the path we take on path in the Pachner graph: these are
quite popular in computational topology and geometry, with the problem to know
what is happening at knocking nodes for instance. More generally, our graph is
less connected than the topological one, how connected is it?

There are other questions concerning our results: given the fact that the
dataset we used is actually sorted by the volume, does the hyperbolic volume
have an impact on the success rate of the methods? How does the approach
perform on very large knots that are not on the census? Are the triangulations
we output suitable for Weeks’ program that computes canonical triangulations
from complete hyperbolic structures [16]?

There are also questions about SnapPy’s method we can look at with ours:
SnapPy re-triangulation is a large number of 2-3 moves performed at random,
why does it miss so much on some specific examples? What makes triangulations
difficult to handle for SnapPy and easy for us, et vise versa?

In conclusion, while this work has a concrete and direct use in 3-manifolds
study, it raises many interesting questions in mathematics, combinatorics,
data science, and algorithmics.

## References

  * [1] Benjamin A. Burton.  The Next 350 Million Knots.  In Sergio Cabello and Danny Z. Chen, editors, 36th International Symposium on Computational Geometry (SoCG 2020), volume 164 of Leibniz International Proceedings in Informatics (LIPIcs), pages 25:1–25:17, Dagstuhl, Germany, 2020. Schloss Dagstuhl–Leibniz-Zentrum für Informatik.
  * [2] Benjamin A. Burton, Ryan Budney, William Pettersson, et al.  Regina: Software for low-dimensional topology.  http://regina-normal.github.io/, 1999–2021.
  * [3] Marc Culler, Nathan M. Dunfield, and Jeffrey R. Weeks.  SnapPy, a computer program for studying the geometry and topology of 3-manifolds.  https://snappy.math.uic.edu/, 1991–2021.
  * [4] David Futer and François Guéritaud.  From angled triangulations to hyperbolic structures, 2010.
  * [5] Joel Hass, Jeffrey C. Lagarias, and Nicholas Pippenger.  The computational complexity of knot and link problems.  J. ACM, 46(2):185–211, 1999.
  * [6] William Jaco and J. Hyam Rubinstein.  0-Efficient Triangulations of 3-Manifolds.  Journal of Differential Geometry, 65(1):61 – 168, 2003.
  * [7] D. Kraft.  A Software Package for Sequential Quadratic Programming.  Deutsche Forschungs- und Versuchsanstalt für Luft- und Raumfahrt Köln: Forschungsbericht. Wiss. Berichtswesen d. DFVLR, 1988.
  * [8] Marc Lackenby.  Word hyperbolic dehn surgery.  Inventiones mathematicae, pages 243–282, 2000.
  * [9] George D. Mostow.  Strong Rigidity of Locally Symmetric Spaces.  Princeton University Press, 1973.
  * [10] Hitoshi Murakami and Yoshiyuki Yokota.  Volume Conjecture for Knots.  Springer, 2018.
  * [11] Barbara Nimershiem.  Geometric triangulations of a family of hyperbolic 3-braids, 2021.
  * [12] Udo Pachner.  P.l. homeomorphic manifolds are equivalent by elementary shellings.  European Journal of Combinatorics, 12(2):129–145, 1991.
  * [13] Igor Rivin.  Euclidean structures on simplicial surfaces and hyperbolic volume.  Annals of Mathematics, 139:553–580, 1994.
  * [14] W. P. Thurston.  The geometry and topology of 3-manifolds, volume 1.  Princeton University Press, Princeton, N.J., 1980.  Electronic version 1.1 - March 2002.
  * [15] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [16] Jeffrey R. Weeks.  Convex hulls and isometries of cusped hyperbolic 3-manifolds.  Topology and its Applications, 52(2):127–149, 1993.
  * [17] Jeffrey R. Weeks.  Computation of hyperbolic structures in knot theory.  Handbook of Knot Theory, 10 2003.
  * [18] Xu Xu.  Combinatorial ricci flow on cusped 3-manifolds.  arXiv preprint arXiv:2009.05477, 2020.


* [17] Rupert L. Frank, María Mar González, Dario D. Monticelli and Jinggang Tan  “An extension problem for the CR fractional Laplacian”  In _Advances in Mathematics_ 270, 2015  DOI: 10.1016/j.aim.2014.09.026
  * [18] Gabriele Grillo, Matteo Muratori and Juan Luis Vázquez  “The porous medium equation on Riemannian manifolds with negative curvature. The large-time behaviour”  In _Advances in Mathematics_ 314, 2017  DOI: 10.1016/j.aim.2017.04.023
  * [19] Pak Tung Ho and Rong Tang  “Fractional Yamabe solitons and fractional Nirenberg problem”  In _Communications on Pure and Applied Analysis_ 20, 2021  DOI: 10.3934/cpaa.2021103
  * [20] Tianling Jin and Jingang Xiong  “A fractional Yamabe flow and some applications”  In _Journal für die reine und angewandte Mathematik (Crelles Journal)_ , 2014  DOI: 10.1515/crelle-2012-0110
  * [21] Arturo De Pablo, Fernando Quirós and Ana Rodríguez  “Nonlocal filtration equations with rough kernels”  In _Nonlinear Analysis, Theory, Methods and Applications_ 137, 2016  DOI: 10.1016/j.na.2016.01.026
  * [22] Arturo De Pablo, Fernando Quirós and Ana Rodríguez  “Regularity theory for singular nonlocal diffusion equations”  In _Calculus of Variations and Partial Differential Equations_ 57, 2018  DOI: 10.1007/s00526-018-1410-2
  * [23] Arturo De Pablo, Fernando Quirós, Ana Rodríguez and Juan Luis Vázquez  “A fractional porous medium equation”  In _Advances in Mathematics_ 226, 2011  DOI: 10.1016/j.aim.2010.07.017
  * [24] Arturo De Pablo, Fernando Quirós, Ana Rodríguez and Juan Luis Vázquez  “A General Fractional Porous Medium Equation”  In _Communications on Pure and Applied Mathematics_ 65, 2012  DOI: 10.1002/cpa.21408
  * [25] Luz Roncal and Sundaram Thangavelu  “Hardy’s inequality for fractional powers of the sublaplacian on the Heisenberg group”  In _Advances in Mathematics_ 302, 2016  DOI: 10.1016/j.aim.2016.07.010
  * [26] Pablo Raúl Stinga and José Luis Torrea  “Extension problem and Harnack’s inequality for some fractional operators”  In _Communications in Partial Differential Equations_ 35, 2010  DOI: 10.1080/03605301003735680
  * [27] Peter Tankov  “Financial Modelling with Jump Processes”  In _Financial Modelling with Jump Processes_ , 2003  DOI: 10.1201/9780203485217
  * [28] Juan Luis Vázquez  “The Porous Medium Equation : Mathematical Theory”  In _The Porous Medium Equation : Mathematical Theory_ , 2007  DOI: 10.1093/acprof:oso/9780198569039.001.0001

Department of Mathematics, Johns Hopkins University, 3400 N. North Charles,
Baltimore, MD 21218

E-mail address<EMAIL_ADDRESS>