# The convergence of discrete uniformizations for genus zero surfaces

Yanwen Luo  Department of Mathematics, Rutgers University, Piscataway, NJ,
08854<EMAIL_ADDRESS>,  Tianqi Wu  Center of Mathematical Sciences
and Applications, Harvard University, Cambridge, MA 02138
<EMAIL_ADDRESS>and  Xiaoping Zhu  Department of Mathematics,
Rutgers University, Piscataway, NJ, 08854<EMAIL_ADDRESS>
###### Abstract.

The notion of discrete conformality proposed by Luo [1] and Bobenko-Pinkall-
Springborn [2] on triangle meshes has rich mathematical theories and wide
applications. Gu et al. [3][4] proved that the discrete uniformizations
approximate the continuous uniformization for closed surfaces of genus $\geq
1$, given that the approximating triangle meshes are reasonably good. In this
paper, we generalize this result to the remaining case of genus zero surfaces,
by reducing it to planar cases via stereographic projections.

###### Contents

  1. 1 Introduction
    1. 1.1 Set up and the main theorem
    2. 1.2 An equivalent formulation
    3. 1.3 Organization of the paper
    4. 1.4 Acknowledgement
  2. 2 Discrete Calculus on Graphs
    1. 2.1 Isoperimetric conditions and a discrete elliptic estimate
    2. 2.2 The differential of the discrete curvature
  3. 3 Stereographic Projections of Triangle Meshes
  4. 4 Estimates on Triangulations
    1. 4.1 Elementary error estimates on triangles
    2. 4.2 Geodesic triangulations of surfaces
  5. 5 Proof of the Main Theorem 1.7
    1. 5.1 Step 1: The estimate of curvatures
    2. 5.2 Step 2: The construction of the flow
    3. 5.3 Step 3: The normalization and the inverse of the stereographic projection

##  1\. Introduction

This work is a continuation of [4], which studies the convergence of discrete
uniformizations to the continuous uniformization for closed surfaces of genus
greater than zero. The notion of discrete conformality and discrete
uniformization discussed here, which is called vertex scaling, were introduced
by Luo [1] and Bobenko-Pinkall-Springborn [2]. Gu-Luo-Wu [3] first proved the
convergence of discrete uniformizations for the unit disk and tori.

In this paper, we prove the convergence result for surfaces of genus zero.
Specifically, for a reasonable geodesic triangulation $T$ on a Riemannian
surface $(M,g)$, the discrete uniformization of its induced polyhedral metric
approximates the uniformization for $(M,g)$ by an error in the order of the
maximal edge length of edges in $T$. There has been extensive study about the
convergence of discrete conformality since Rodin-Sullivan’s seminal work [5]
on the convergence of circle packings. Other convergence results about the
vertex scaling can also be found in [6][7][8].

###  1.1. Set up and the main theorem

Suppose $M$ is a compact orientable surface, possibly with boundary, and $T$
is a triangulation of $M$, which is always assumed to be a simplicial complex
such that its one-skeleton is a $4$-vertex-connected graph. Denote
$V(T),E(T),F(T)$ as the sets of vertices, edges, and triangles of $T$
respectively. Further denote $\text{int}(T)$ as the set of interior vertices
of $T$, and $\text{bdy}(T)$ as the set of boundary vertices of $T$. If $M$ is
equipped with a smooth Riemannian metric $g$, $T$ is a geodesic triangulation
if any edge in $T$ is a shortest geodesic arc in $(M,g)$.

An _admissible edge length function_ $l\in\mathbb{R}_{>0}^{E(T)}$ of $T$
satisfies the triangle inequality $l_{ij}+l_{jk}>l_{ik}$ for any triangle
$\triangle ijk\in F(T)$. For a geodesic triangulation $T$ on a Riemannian
surface $(M,g)$, we can naturally define an edge length function $l$ using the
geodesic lengths of the edges. Given an admissible edge length function $l$,
we can construct a piecewise Euclidean triangle mesh $(T,l)_{E}$ by
isometrically gluing the Euclidean triangles with the edge lengths defined by
$l$ along pairs of edges. Similarly, a spherical triangle mesh $(T,l)_{S}$ can
be constructed by replacing Euclidean triangles with spherical triangles of
the same edge lengths, provided that $l_{ij}+l_{jk}+l_{ki}<\pi$ for any
triangle $\triangle ijk$ in $F(T)$. For given $(T,l)_{E}$ or $(T,l)_{S}$, we
often use $\theta^{i}_{jk}$ to denote the inner angle at $i$ in the triangle
$\triangle ijk$. We also define the discrete curvature $K_{i}$ at a vertex
$i\in V(T)$ as

$K_{i}=\left\\{\begin{array}[]{ll}2\pi-\sum_{jk\in E:ijk\in
F}\theta^{i}_{jk},&\mbox{ if $i\in\text{int}(T)$},\\\ \pi-\sum_{jk\in E:ijk\in
F}\theta^{i}_{jk},&\mbox{ if $i\in\text{bdy}(T)$}.\end{array}\right.$

A Euclidean (resp. spherical) triangle mesh is globally flat (resp. globally
spherical) if and only if $K_{i}=0$ for any $i\in\text{int}(T)$.

######  Definition 1.1.

Given a triangulation $T$, a _discrete conformal factor_ $u$ is a real-valued
function on $V(T)$. For the Euclidean case, $(T,l)_{E}$ and
$(T,l^{\prime})_{E}$ are called _discrete conformal_ if for some discrete
conformal factor $u$,

(1) $l^{\prime}_{ij}=e^{\frac{1}{2}(u_{i}+u_{j})}l_{ij}$

for any $ij\in E(T)$. For the spherical case, $(T,l)_{S}$ and
$(T,l^{\prime})_{S}$ are _discrete conformal_ if for some discrete conformal
factor $u$,

(2)
$\sin\frac{l^{\prime}_{ij}}{2}=e^{\frac{1}{2}(u_{i}+u_{j})}\sin\frac{l_{ij}}{2}$

for any $ij\in E(T)$.

We denote $l^{\prime}=u*l$ if equation (1) holds, and $l^{\prime}=u*_{s}l$ if
equation (2) holds. Given a triangle mesh $(T,l)_{E}$ (resp. $(T,l)_{S}$), let
$\theta^{i}_{jk}(u)$ and $K_{i}(u)$ denote the corresponding inner angle at
$i$ in triangle $\triangle ijk$ and the discrete curvature at $i$ respectively
in $(T,u*l)_{E}$ (resp. $(T,u*_{s}l)_{S}$). If $(T,l)_{S}$ is a topological
sphere, $u$ is called a _discrete uniformization factor_ for $(T,l)_{S}$, if
$(T,u*_{s}l)_{S}$ is isometric to the unit sphere, which is equivalent to that
the discrete curvature $K(u):=[K_{i}(u)]_{i\in V}$ is zero. A triangle mesh
$(T,l)_{E}$ (resp. $(T,l)_{S}$) is called _strictly Delaunay_ if for any edge
$ij$ in T, two adjacent triangles $\triangle ijk,\triangle ijk^{\prime}$
containing $ij$ satisfies

(3)
$\theta^{k}_{ij}+\theta^{k^{\prime}}_{ij}<\theta^{i}_{jk}+\theta^{i}_{jk^{\prime}}+\theta^{j}_{ik}+\theta^{j}_{ik^{\prime}}.$

We quantify the regularity of a triangle mesh as follows.

######  Definition 1.2.

A triangle mesh $(T,l)_{E}$ (resp. $(T,l)_{S}$) is called _$\epsilon$
-regular_ if

  1. (a)

any inner angle $\theta^{i}_{jk}\geq\epsilon$, and

  2. (b)

for any adjacent triangles $\triangle ijk$ and $\triangle ijk^{\prime}$,
$\theta^{k}_{ij}+\theta^{k^{\prime}}_{ij}\leq\pi-\epsilon$.

Condition (a) requires that any triangle is away from degenerating, and
condition (b) requires the triangle mesh to be uniformly strictly Delaunay.
Let $|x|=\max_{i\in I}|x_{i}|$ denote the maximal norm of a vector
$x\in\mathbb{R}^{I}$ in a finite dimensional vector space. Let
$\hat{\mathbb{C}}$ be the standard Riemann sphere, which can be identified
with the unit sphere $\mathbb{S}^{2}$ in $\mathbb{R}^{3}$ by the stereographic
projection. The main theorem of the paper is

######  Theorem 1.3.

Suppose $(M,g)$ is a closed smooth Riemannian surface of genus zero with three
marked points $X,Y,Z$, and $\bar{u}\in C^{\infty}(M)$ is the unique
uniformization conformal factor such that $(M,e^{2\bar{u}}g)$ is isometric to
the unit sphere $\mathbb{S}^{2}\cong\hat{\mathbb{C}}$ through map $\phi$,
where $\phi(Z)=0$, $\phi(Y)=1$, $\phi(X)=\infty$. Assume $T$ is a geodesic
triangulation of $(M,g)$, and $l\in\mathbb{R}_{>0}^{E(T)}$ denotes its edge
length in $(M,g)$. Then for any $\epsilon>0$, there exists
$\delta=\delta(M,g,X,Y,Z,\epsilon)>0$ such that if $(T,l)_{S}$ is
$\epsilon$-regular and $|l|\leq\delta$, then

  1. (a)

there exists a unique discrete conformal factor $u$ on $V(T)$, such that
$(T,u*_{s}l)_{S}$ is strictly Delaunay and isometric to the unit sphere
through a map $\psi$ such that $\psi(Z)=0$, $\psi(Y)=1$, and $\psi(X)=\infty$,
and

  2. (b)

$|u-\bar{u}|_{V(T)}|\leq C|l|$ for some constant $C=C(M,g,X,Y,Z,\epsilon)>0$.

######  Remark 1.4.

The uniqueness part of the theorem is already known as a consequence of
Springborn’s Theorem 10.5 in [9], which is equivalent to Rivin’s earlier
result on hyperbolic polyhedral realization in [10].

###  1.2. An equivalent formulation

Springborn [11][9] and Bobenko et al. [2] proposed an another notion of
discrete uniformization, which is for Euclidean triangle meshes that are
homeomorphic to a sphere. We adapt their definitions as follows.

Let $\mathcal{P}$ be the set of the compact convex polyhedral surfaces $P$,
satisfying that

(a) $P$ is the boundary of the convex hull of a finite subset of
$\mathbb{S}^{2}$, and

(b) $0$ is strictly inside $P$, and

(c) each face of $P$ is a triangle.

Given $P\in\mathcal{P}$, denote $V(P)$ as the set of its vertices, and $T_{P}$
as the natural triangulation of $P$ where each triangle is a face of $P$, and
$l_{P}\in\mathbb{R}^{E(T_{P})}$ as the edge length of $T_{P}$ on $P$. For a
Euclidean triangle mesh $(T,l)_{E}$, which is a topological sphere, we say
that $u$ is a _discrete uniformization factor_ of $(T,l)_{E}$ if $(T,u*l)_{E}$
is isometric to some $P\in\mathcal{P}$, through a map $\varphi$ such that
$\varphi(T)=T_{P}$.

We call a geodesic triangulation of the unit sphere _strictly Delaunay_ if the
circumference circle of each triangle contains no other vertex. It is well-
known that this spherical empty circle condition is equivalent to the
condition that if for any edge $ij$ in T, two adjacent triangles $\triangle
ijk,\triangle ijk^{\prime}$ containing $ij$ satisfies

(4)
$\theta^{k}_{ij}+\theta^{k^{\prime}}_{ij}<\theta^{i}_{jk}+\theta^{i}_{jk^{\prime}}+\theta^{j}_{ik}+\theta^{j}_{ik^{\prime}}.$

If we consider this geodesic triangulation as a sphercial triangle mesh
$(T,l)_{S}$ where $l(e)$ is the geodesic arc length of edge $e$, then the
above condition is just condition (3).

The central projection $p:(x,y,z)\mapsto(x,y,z)/\sqrt{x^{2}+y^{2}+z^{2}}$
naturally gives rise to a bijection between $\mathcal{P}$ and the set of
strictly Delaunay triangulations of the unit sphere. Here we always assume
that a triangulation of $\mathbb{S}^{2}$ is a geodesic triangulation and each
triangle is a proper subset of a hemisphere. For a vector
$x\in\mathbb{R}^{I}$, we denote $\sin x$ as the vector in $\mathbb{R}^{I}$
such that $(\sin x)_{i}=\sin(x_{i})$.

######  Proposition 1.5.

Let $P\mapsto p(T_{P})$ be a bijection from $\mathcal{P}$ to the set of
strictly Delaunay triangulations of $\mathbb{S}^{2}$. Further we have that
$l_{P}=2\sin\frac{l}{2}$ where $l$ denotes the geodesic edge lengths of
$p(T_{P})$ on $\mathbb{S}^{2}$.

###### Proof.

Given $P\in\mathcal{P}$, for any triangle $\triangle ijk$, other vertices are
on one side of the plane $\triangle ijk$ lies in by the convexity of $P$.
Therefore other vertices lies outside the circumference circle of spherical
triangle $p(\triangle ijk)$. So $p(T_{P})$ is a strictly Delaunay
triangulation of $\mathbb{S}^{2}$. See Figure 1 for illustrations.

If $T$ is a strictly Delaunay triangulation of $\mathbb{S}^{2}$, we construct
a polyhedral surface $P$ as the union of all flat triangles $\triangle ijk$
where $i,j,k$ are the three vertices of a triangle in $T$. Then $p|_{P}$ is a
homeomorphism from $P$ to $\mathbb{S}^{2}$. Since $T$ satisfies the empty
circle property, we conclude that the dihedral angle on any edge $ij\in E(T)$
is less than $\pi$ by the similar argument in the above paragraph. So $P$ is a
convex polyhedral surface (See Lemma 6.1 in [12] for example).

Since all the vertices of $P$ is on the unit sphere, therefore $P$ satisfies
the condition $(a)$ of set $\mathcal{P}$. It is not hard to see that $P$
satisfies condition $(b)$ and $(c)$. The condition $l_{P}=2\sin\frac{l}{2}$
follows easily from the construction of $P$. ∎

Figure 1. Equivalence between local Delaunay condition and local convexity.

As a consequence, we obtain an equivalence between the two notions of discrete
uniformizations.

######  Corollary 1.6.

Assume $T$ is topologically a sphere, then $u$ is a discrete uniformization
factor of $(T,2\sin\frac{l}{2})_{E}$ if and only if $u$ is a discrete
uniformization factor of $(T,l)_{S}$ and $(T,u*_{s}l)_{S}$ is strictly
Delaunay.

By such an equivalence, we can reformulate our main theorem as follows.

######  Theorem 1.7.

Suppose $(M,g)$ is a closed smooth Riemannian surface of genus zero with three
marked points $X,Y,Z$, and $\bar{u}\in C^{\infty}(M)$ is the unique
uniformization conformal factor such that $(M,e^{2\bar{u}}g)$ is isometric to
the unit sphere $\mathbb{S}^{2}\cong\hat{\mathbb{C}}$ through map $\phi$, and
$\phi(Z)=0$, $\phi(Y)=1$, $\phi(X)=\infty$. Assume $T$ is a geodesic
triangulation of $(M,g)$, and $l\in\mathbb{R}_{>0}^{E(T)}$ denotes its edge
length in $(M,g)$. Then for any $\epsilon>0$, there exists
$\delta=\delta(M,g,X,Y,Z,\epsilon)>0$ such that if $(T,l)_{S}$ is
$\epsilon$-regular and $|l|\leq\delta$, then

  1. (a)

there exists a unique discrete conformal factor $u$ on $V(T)$, such that
$(T,u*(2\sin\frac{l}{2}))_{E}$ is isometric to some $P\in\mathcal{P}$ through
a map $\psi$ such that $\psi(Z)=0$, $\psi(Y)=1$, and $\psi(X)=\infty$, and

  2. (b)

$|u-\bar{u}|_{V(T)}|\leq C|l|$ for some constant $C=C(M,g,X,Y,Z,\epsilon)>0$.

We will prove this new version of our main theorem. By the stereographic
projection, we can consider triangulations of a flat polygon, instead of the
polyhedrons inscribed in the unit sphere. To obtain a satisfactory flat
triangle mesh, we adapt the idea in [4] and construct a discrete curvature
flow on the triangle mesh. The estimate in part (b) essentially follows from a
discrete elliptic estimate on the flow.

###  1.3. Organization of the paper

In Section 2, we will introduce the discrete calculus on graphs and the key
discrete elliptic estimate. In Section 3 we will discuss the stereographic
projection and the one-to-one correspondence between the convex polyhedral
surfaces inscribed in the unit sphere and the the Delaunay triangulations of
convex polygons. In Section 4, we will introduce some useful and elementary
estimates about triangulations. The proof of the main theorem 1.7 will be
given in Section 5.

###  1.4. Acknowledgement

The work is supported in part by NSF 1719582, NSF 1760471, NSF 1760527, NSF
DMS 1737876, and NSF 1811878.

##  2\. Discrete Calculus on Graphs

Assume $G=(V,E)$ is an undirected connected simple graph, on which we will
frequently consider vectors in $\mathbb{R}^{V}$, $\mathbb{R}^{E}$ and
$\mathbb{R}^{E}_{A}$. Here $\mathbb{R}^{E}$ and $\mathbb{R}^{E}_{A}$ are both
vector spaces of dimension $|E|$ such that

  1. (a)

a vector $x$ in $\mathbb{R}^{E}$ is represented symmetrically, i.e.,
$x_{ij}=x_{ji}$, and

  2. (b)

a vector $x$ in $\mathbb{R}^{E}_{A}$ is represented anti-symmetrically, i.e.,
$x_{ij}=-x_{ji}$.

A vector $x$ in $\mathbb{R}^{E}_{A}$ is also called a _flow_ on $G$. An _edge
weight_ $\eta$ on $G$ is a vector in $\mathbb{R}^{E}$. Given an edge weight
$\eta$ and a vector $x\in\mathbb{R}^{V}$, its _gradient_ $\nabla
x=\nabla_{\eta}x$ is a flow in $\mathbb{R}^{E}_{A}$ defined as

$(\nabla x)_{ij}=\eta_{ij}(x_{j}-x_{i}).$

Given a flow $x\in\mathbb{R}^{E}_{A}$, its _divergence_ $div(x)$ is a vector
in $\mathbb{R}^{V}$ defined as

$div(x)_{i}=\sum_{j\sim i}x_{ij}.$

Given an edge weight $\eta$, the associated _Laplacian_
$\Delta=\Delta_{\eta}:\mathbb{R}^{V}\to\mathbb{R}^{V}$ is defined as $\Delta
x=\Delta_{\eta}x=div(\nabla_{\eta}x)$, i.e.,

$(\Delta x)_{i}=\sum_{j\sim i}(\nabla x)_{ij}=\sum_{j\sim
i}\eta_{ij}(x_{j}-x_{i}).$

The Laplacian is a linear transformation on $\mathbb{R}^{V}$, and could be
identified as a symmetric $|V|\times|V|$ matrix.

###  2.1. Isoperimetric conditions and a discrete elliptic estimate

Analogous to the isoperimetric condition on a Riemannian surface, we introduce
the notion of $C$-isoperimetric conditions for a graph $G=(V,E)$ with a
positive vector $l\in\mathbb{R}^{E}_{>0}$. Given $V_{0}\subset V$, we denote

$\partial V_{0}=\\{ij\in E:i\in V_{0},\text{ and }j\notin V_{0}\\}$

(See Figure 2 below), and define the l-perimeter of $V_{0}$ and the l-area of
$V_{0}$ as

$|\partial V_{0}|_{l}=\sum_{ij\in\partial
V_{0}}l_{ij},\quad\text{and}\quad|V_{0}|_{l}=\sum_{ij\in E:i,j\in
V_{0}}l_{ij}^{2}.$

Given a constant $C>0$, such a pair $(G,l)$ is called $C$-_isoperimetric_ if
for any $V_{0}\subset V$,

$\min\\{|V_{0}|_{l},|V|_{l}-|V_{0}|_{l}\\}\leq C|\partial V_{0}|_{l}^{2}.$

Figure 2. The $C$-isoperimetric condition on graphs.

We will see in Section 4 that a uniform $C$-_isoperimetric_ condition is
satisfied by $\epsilon-$ regular triangle meshes approximating a closed smooth
surface. The following discrete elliptic estimate is a key tool to prove our
theorem of convergence, and is reformulated from the second conclusion of
Lemma 2.3 in [4].

######  Lemma 2.1.

Given a constant $C>0$ and a $C$-isoperimetric pair $(G,l)$, consider the
equation

(5) $(D-\Delta_{\eta})u=div(x)+y$

on $G$ where

  1. (i)

$\eta\in\mathbb{R}^{E}$ is an edge weight such that for any $ij\in E$

$\eta_{ij}\geq\frac{1}{C},$

and

  2. (ii)

$x\in\mathbb{R}^{E}_{A}$ is a flow such that for any $ij\in E$

$|x_{ij}|\leq Cl_{ij}^{2},$

and

  3. (iii)

$D\in\mathbb{R}^{V\times V}$ is a nonzero nonnegative diagonal matrix, and

  4. (iv)

$y\in\mathbb{R}^{V}$ satisfies that for any $i\in V$

$|y_{i}|\leq C\cdot D_{ii}|l|\cdot|V|_{l}^{\frac{1}{2}}.$

Then the solution $u\in\mathbb{R}^{V}$ of equation (5) satisfies that

$|u|\leq C^{\prime}|l|\cdot|V|_{l}^{\frac{1}{2}}.$

for some constant $C^{\prime}=C^{\prime}(C)>0$.

###  2.2. The differential of the discrete curvature

There is an explicit formula for the infinitesimal changes of the discrete
curvature, as the mesh deforms in its discrete conformal class.

######  Proposition 2.2 (Proposition 4.1.6 in [2]).

Given a triangle mesh $(T,l)_{E}$, if $u*l$ is an admissible edge length,
denote $\theta^{i}_{jk}(u)$ as the inner angle at $i$ in triangle $\triangle
ijk$ in $(T,u*l)_{E}$, and $K(u)\in\mathbb{R}^{V}$ as the discrete curvature
in $(T,u*l)_{E}$. We also define the cotangent weight $\eta\in\mathbb{R}^{E}$
as

$\eta_{ij}(u)=\frac{1}{2}\cot\theta^{k}_{ij}(u)+\frac{1}{2}\cot\theta^{k^{\prime}}_{ij}(u)$

if $\triangle ijk,\triangle ijk^{\prime}\in F(T)$ sharing edge $ij$, and

$\eta_{ij}(u)=\frac{1}{2}\cot\theta_{ij}^{k}$

if $\triangle ijk\in F(T)$ and $ij$ is a boundary edge of $T$. Then we have

$\frac{\partial K}{\partial u}(u)=-\Delta_{\eta(u)}.$

##  3\. Stereographic Projections of Triangle Meshes

We will use the stereographic projection to connect convex polyhedral surfaces
inscribed in the unit spheres with triangulations of planar convex polygons.
Denote $N$ as the north pole $(0,0,1)$ of the unit sphere $\mathbb{S}^{2}$.
The stereographic projection $p_{N}$ is a map from
$\mathbb{R}^{3}\backslash\\{z=1\\}$ to the $xy$-plane, which is identified as
$\mathbb{R}^{2}$ or $\mathbb{C}$. The map $p_{N}$ is defined as

$p_{N}(x,y,z)=\frac{x}{1-z}+i\frac{y}{1-z}.$

It is well known that the restriction of $p_{N}$ on
$\mathbb{S}^{2}\backslash\\{N\\}$ is a conformal diffeomorphism to
$\mathbb{R}^{2}$, and maps any circle to a circle or a straight line. Given a
convex polyhedral surface $P$ in $\mathcal{P}$ such that $N$ is a vertex of
$P$, denote $\mathring{T}_{P}$ as the subtriangulation of $T_{P}$ with the
open 1-star neighborhood of $N$ in $T_{P}$ being removed. In this case,
$\mathring{P}$ denotes the carrier of $\mathring{T}_{P}$ and is a topological
closed disk. We use $|\cdot|_{2}$ to denote the standard $l^{2}$-norm.

######  Lemma 3.1.

(a) Assume $P\in\mathcal{P}$ and $P$ contains $N$ as a vertex, then $p_{N}$ is
injective on $\mathring{P}$, and $Q=p_{N}(\mathring{P})$ is a convex polygon,
and $T_{Q}=p_{N}(\mathring{T}_{P})$ is a geodesic triangulation of $Q$.
Further, if we naturally identify $\mathring{T}_{P}$ and $T_{Q}$, and denote
$l_{P}$ (resp. $l_{Q}$) as its edge length on $P$ (resp. $Q$), then

(6) $l_{Q}=w*l_{P}\text{\quad}{where}\quad
w_{i}=\log\frac{2}{|i-N|_{2}^{2}}=\log\frac{|p_{N}(i)|_{2}^{2}+1}{2},\quad\forall
i\in V(\mathring{T}_{P}).$

(b) Assume $Q$ is a convex polygon in $\mathbb{R}^{2}$, and $T_{Q}$ is a
strictly Delaunay triangulation of $Q$ such that $0$ is an interior vertex and
$K_{i}>0$ for any boundary vertex $i$ in $V(T_{Q})$. Then there exists a
convex polyhedral surface $P\in\mathcal{P}$ such that $N\in P$ and
$p_{N}(\mathring{P})=Q$ and $p_{N}(\mathring{T}_{P})=T_{Q}$.

###### Proof.

(a) Let us first prove the injectivity by contradiction. Suppose $x,y$ are two
different points on $\mathring{P}$ such that $p_{N}(x)=p_{N}(y)$. Then $N,x,y$
are co-linear and pairwise different. Without loss of generality, assume $y$
is between $x$ and $N$. Then it is not difficult to show that the line segment
$\overline{xN}\subset P$. Suppose $\triangle ijk$ is a face of $P$ containing
the line segment $\overline{xN}$. Then $N$ has to be one of the vertex of
$\triangle ijk$, and $x$, $y$ are contained in the edge in $\triangle ijk$
opposite to $N$. But this implies that $N,x,y$ are not co-linear, which leads
to a contradiction.

So $Q$ is a polygon, and $T_{Q}$ is a geodesic triangulation of $Q$. Any inner
angle of the polygon $Q$ is less than $\pi$, since the dihedral angle on any
edge $Ni\in E(T_{P})$ is less than $\pi$. Equation (6) can be proved by a
standard computation.

(b) We will first construct a polyhedron $P$ and then show that it is
satisfactory. The set of vertices of $P$ is given by
$V_{P}=(p_{N}|_{\mathbb{S}^{2}})^{-1}(V(T_{Q}))\cup\\{N\\}$, and the set of
faces of $P$ is given by a set of flat triangles in $\mathbb{R}^{3}$ with
vertices in $V_{P}$. The $\triangle ijk$ is a triangle in $P$ if and only if
$p_{N}(\triangle ijk)$ is a triangle in $T_{Q}$, or $\\{i,j,k\\}=\\{N,x,y\\}$
where $p_{N}(xy)$ is a boundary edge of $T_{Q}$.

It is ordinary to verify that the union $P$ of such triangles is a topological
sphere, and these flat triangles naturally give a triangulation $T_{P}$ of
$P$. It remains to show that $P\in\mathcal{P}$, or indeed that any dihedral
angle in $T_{P}$ is less than $\pi$.

Assume $ij$ is an edge in $T_{P}$. If $i=N$, the dihedral angle at $ij$ is
less than $\pi$ because the discrete curvature at $p_{N}(j)$ in $T_{Q}$ is
greater than $0$. If $p_{N}(ij)$ is a boundary edge in $T_{Q}$, assume
$\triangle ijk\in F(T_{P})$ and $k\neq N$, and then the dihedral angle at $ij$
is less than $\pi$ because $p_{N}(k)$ and $0$ are in the same half plane
divided by $p_{N}(ij)$. Now we can assume that $p_{N}(ij)$ is an inner edge in
$T_{Q}$, and $p_{N}(\triangle ijk),p_{N}(\triangle ijk^{\prime})$ are two
triangles in $T_{Q}$. Since $T_{Q}$ is strictly Delaunay, $p_{N}(k^{\prime})$
is strictly outside of the circumcircle of $p_{N}(\triangle ijk)$. Since
$p_{N}|_{\mathbb{S}^{2}}$ preserves circles, $k^{\prime}$ is strictly outside
of the spherical circumcircle of $\\{i,j,k\\}$ on $\mathbb{S}^{2}$. So the
dihedral angle at $ij$ is less than $\pi$. ∎

In the following lemma we prove that the stereogarphic projection preserves
the $\epsilon$-regularity.

######  Lemma 3.2.

Assume $P\in\mathcal{P}$, and $N\in P$, and $T=p(T_{P})$ is a geodeisc
triangulation of $\mathbb{S}^{2}$, and $Q=p_{N}(\mathring{P})$, and
$T_{Q}=p_{N}(\mathring{T}_{P})$, and $l$ (resp. $l_{Q}$) denotes the edge
length of $T$ (resp. $T_{Q}$) on $\mathbb{S}^{2}$ (resp. $Q$). Then for any
$\epsilon>0$, there exists constants
$\epsilon^{\prime}=\epsilon^{\prime}(\epsilon)>0$ and
$\delta=\delta(\epsilon)>0$ such that if $(T,l)_{S}$ is $\epsilon$-regular and
$|l|<\delta$, then $(T_{Q},l_{Q})_{E}$ is $\epsilon^{\prime}$-regular, and
$K_{i}\geq\epsilon^{\prime}$ for any boundary vertex $i$ in $T_{Q}$.

###### Proof.

Let $\theta_{ij}^{k}$ denote the inner angles in $(T,l)_{S}$, and
$\phi_{ij}^{k}$ denote the inner angles in $(T_{Q},l_{Q})_{E}$. We need to
prove following three statements: (a) $\phi_{ij}^{k}$ are bounded below by
$\epsilon^{\prime}>0$, and (b) $T_{Q}$ is strictly Delaunay with angle sums
$\phi_{ij}^{k}+\phi_{ij}^{k^{\prime}}$ bounded above by
$\pi-\epsilon^{\prime}$, and (c) $K_{i}\geq\epsilon^{\prime}$ for any boundary
vertex $i$ in $T_{Q}$.

Consider a pair of triangles $\triangle ijk$ and $\triangle ijk^{\prime}$ in
$(T,l)_{S}$, and by assumption
$\theta_{ij}^{k}+\theta_{ij}^{k^{\prime}}\leq\pi-\epsilon$. Let $\Theta_{ij}$
be the intersecting angle of the circumcircles of two triangles $\triangle
ijk$ and $\triangle ijk^{\prime}$ on $\mathbb{S}^{2}$. It is elementary to
show that

$\Theta_{ij}=\theta_{jk}^{i}+\theta_{ik}^{j}+\theta_{jk^{\prime}}^{i}+\theta^{j}_{ik^{\prime}}-(\theta_{ij}^{k}+\theta_{ij}^{k^{\prime}})>2\pi-2(\theta_{ij}^{k}+\theta_{ij}^{k^{\prime}})\geq
2\epsilon.$

The stereographic projection preserves angles and circles, so the intersecting
angle of the circumcircles of $p_{N}(\triangle ijk)$ and $p_{N}(\triangle
ijk^{\prime})$ in $T_{Q}$ is also $\Theta_{ij}$, if $N$ is not contained in
$\triangle ijk\cup\triangle ijk^{\prime}$. Then it is also ordinary to show
that this intersecting angle is

$\Theta_{ij}=\phi_{jk}^{i}+\phi_{ik}^{j}+\phi_{jk^{\prime}}^{i}+\phi^{j}_{ik^{\prime}}-(\phi_{ij}^{k}+\phi_{ij}^{k^{\prime}})=2\pi-2(\phi_{ij}^{k}+\phi_{ij}^{k^{\prime}}).$

Therefore part (b) is true by

$\phi_{ij}^{k}+\phi_{ij}^{k^{\prime}}=\pi-\frac{\Theta_{ij}}{2}\leq\pi-\epsilon.$

If $i=N$, then the circumcircles of $\triangle ijk$ and $\triangle
ijk^{\prime}$ are mapped to straight lines $p_{N}(jk)$ and
$p_{N}(jk^{\prime})$. So the angle between $p_{N}(jk)$ and
$p_{N}(jk^{\prime})$, or the inner angle of the polygon $Q$ at $p_{N}(j)$, is
equal to $\pi-\Theta_{ij}\leq\pi-2\epsilon$. So part (c) is true.

Figure 3. Projection of angles under stereographic projection

Now we prove part (a). Assume $\triangle ijk$ is a triangle in $T$ not
containing $N$, and $C$ is the circumcircle of $\triangle ijk$ on
$\mathbb{S}^{2}$, and $C^{\prime}$ is the circle on $\mathbb{S}^{2}$
containing $\\{j,k,N\\}$. Then $p_{N}(C)$ is the circumcircle of
$p_{N}(\triangle ijk)$ and $p_{N}(C^{\prime})$ is the straight line
$p_{N}(jk)$, and the intersecting angle of them is equal to the intersecting
angle of $C$ and $C^{\prime}$. It is elementary to show that $\phi^{i}_{jk}$
is equal to an intersecting angle of $p_{N}(C)$ and $p_{N}(C^{\prime})$, i.e.,
an intersecting angle of $C$ and $C^{\prime}$. See Figure 3 for illustration.
We only need to show that the intersecting angle of $C$ and $C^{\prime}$ are
at least $\epsilon^{\prime}$ for some constant
$\epsilon^{\prime}(\epsilon)>0$, when $|l|<\delta$ for some constant
$\delta(\epsilon)>0$.

Denote $R$ as the spherical radius of the cirlce $C$. Since $(T,l)_{S}$ is
$\epsilon$-regular and $|l|<\delta$, the degree (valence) of any vertex in $T$
is at most $\lfloor 2\pi/\epsilon\rfloor$, and $l_{ij},l_{ik},l_{jk}$ are at
least $r_{1}R$ for some constant $r_{1}=r_{1}(\epsilon)$. Further it is not
difficult to show that there exists a constant $r_{2}(\epsilon)>0$ such that
the 1-star neighborhood of $i$ in $T$ contains the open spherical disk $U_{i}$
centered at $i$ with radius $r_{2}R$. So $N\notin U_{i}$. We define $U_{j}$
and $U_{k}$ similarly.

Assume $C_{k}$ is the circumcircle of the triangle in $T$ that is adjacent to
$\triangle ijk$ along the edge $ij$. Then the intersecting angle between $C$
and $C_{k}$ is $\Theta_{ij}\geq 2\epsilon$. Assume $C_{k}^{\prime}$ is the
circle on $\mathbb{S}^{2}$ such that $i,j\in C_{k}^{\prime}$ and the
intersecting angle between $C$ and $C_{k}^{\prime}$ is equal to $2\epsilon$.
Denote $D_{k}$ (resp. $D_{k}^{\prime}$, $D$) as the open spherical disk
bounded by $C_{k}$ (resp. $C_{k}^{\prime}$, $C$). Then $D_{k}^{\prime}\subset
D\cup D_{k}$ and has a diameter less than $r_{3}R$ for some constant
$r_{3}(\epsilon)>0$. Then $N\notin D_{k}$ and $N\notin D$ by the convexity of
$P$, and so $N\notin D_{k}^{\prime}$. Define $D_{i}^{\prime}$ and
$D_{j}^{\prime}$ similarly.

Without loss of generality, assume $j^{\prime}$ is the opposite point of $j$
on $\mathbb{S}^{2}$, and $j^{\prime}\neq N$. Denote $X$ as the tangent plane
of $\mathbb{S}^{2}$ at $j$ in $\mathbb{R}^{3}$, and $p^{\prime}$ as the
projection map centered at $j^{\prime}$ and mapping
$\mathbb{S}^{2}\backslash\\{j^{\prime}\\}$ to $X$. Since $p^{\prime}$
preserves the angles and disks, $p^{\prime}(D)$, $p^{\prime}(D_{i}^{\prime})$,
$p^{\prime}(D_{j}^{\prime})$, $p^{\prime}(D_{k}^{\prime})$,
$p^{\prime}(U_{i})$, $p^{\prime}(U_{j})$, $p^{\prime}(U_{k})$ are all disks on
$X$, and the intersecting angle between $p^{\prime}(D)$ and
$p^{\prime}(D_{i})$ (resp. $p(D_{j})$, $p^{\prime}(D_{k})$) is $2\epsilon$,
and we only need to show that the intersecting angle between $p^{\prime}(C)$
and $p^{\prime}(C^{\prime})$ is at least $\epsilon^{\prime}$ for some constant
$\epsilon^{\prime}(\epsilon)>0$. The projection $p^{\prime}$ is very close to
an isometry near the point $j$, so if $\delta=\delta(\epsilon)>0$ is
sufficiently small,

$\frac{1}{2}d(x,y)\leq d(p^{\prime}(x),p^{\prime}(y))\leq 2d(x,y)$

for any $x,y\in D\cup D_{i}^{\prime}\cup D_{j}^{\prime}\cup D_{k}^{\prime}\cup
U_{i}\cup U_{j}\cup U_{k}$. Assume $R^{\prime}$ is the radius of
$p^{\prime}(D)$ and it is not difficult to show that the radii of
$p^{\prime}(U_{i})$ and $p^{\prime}(U_{j})$ and $p^{\prime}(U_{k})$ are at
least $r_{2}R^{\prime}/4$, and $d(p^{\prime}(a),p^{\prime}(b))\geq
r_{1}R^{\prime}/4$ for any two different vertices $a,b$ in $\\{i,j,k\\}$.

Figure 4. Seven circles.

By a scaling, it suffices to prove the following claim. Assume $\triangle xyz$
is a triangle inscribed in the unit circle $\mathbb{S}$ in the plane, and all
the edge lengths are at least $r_{1}/4$, and $C_{x}$ is the circle such that
$y,z\in C_{x}$ and $x$ is not inside $C_{x}$ and the intersecting angle
between $C$ and $C_{x}$ is $2\epsilon$, and $C_{y}$ and $C_{z}$ are defined
similarly, and $U_{x}$ is the open disk centered at $x$ with radius $r_{2}/4$,
and $U_{y}$ and $U_{z}$ are defined similarly, and $N^{\prime}$ is a point
that is not strictly inside of the unit circle or $C_{x}$ or $C_{y}$ or
$C_{z}$ or $U_{x}$ or $U_{y}$ or $U_{z}$, and $C_{N^{\prime}}$ is the circle
(or the straight line) passing through $y,z,N^{\prime}$, then the intersecting
angle $\theta$ between $C$ and $C_{N^{\prime}}$ is at least
$\epsilon^{\prime}$ for some constant $\epsilon^{\prime}(\epsilon)>0$. See
Figure 4 for illustration.

If the above claim is not true, for some $\epsilon>0$, one can pick a sequence
$(x_{n},y_{n},z_{n},N^{\prime}_{n})\in(\mathbb{R}^{2})^{4}$ such that the
resulted intersecting angle $\theta_{n}$ goes to $0$. By picking a
subsequence, we may assume $x_{n}\rightarrow x\in\mathbb{S}$ and
$y_{n}\rightarrow y\in\mathbb{S}$ and $z_{n}\rightarrow z\in\mathbb{S}$ and
$N^{\prime}_{n}\rightarrow N^{\prime}\in\mathbb{R}^{2}\cup\\{\infty\\}$, and
then by the continuity $(x,y,z,N^{\prime})$ satisfies the conditions in the
claim, and the resulted intersecting angle is $\theta=0$. This means that
$N^{\prime}$ is on the unit circle. However this is impossible because by the
continuity

$N^{\prime}\notin D_{x}\cup D_{y}\cup D_{z}\cup U_{x}\cup U_{y}\cup
U_{z}\supset\mathbb{S}$

where $D_{x}$ (resp. $D_{y}$, $D_{z}$) is the open disk bounded by $C_{x}$
(resp. $C_{y}$, $C_{z}$). ∎

##  4\. Estimates on Triangulations

In this section we will introduce some elementary estimates on triangles, and
properties of triangulations on a Rimannian surface.

###  4.1. Elementary error estimates on triangles

Denote $|\triangle ABC|$ (resp. $|\triangle A^{\prime}B^{\prime}C^{\prime}|$)
as the area of $\triangle ABC$ (resp. $\triangle
A^{\prime}B^{\prime}C^{\prime}$). We have the following estimate on the change
of the angles and area of a Euclidean triangle.

######  Lemma 4.1 (Lemma 3.5 in [4]).

Given a Euclidean triangle $\triangle ABC$ with edge lengths $a$, $b$, and
$c$, if all the angles in $\triangle ABC$ are at least $\epsilon>0$, and
$\delta<\epsilon^{2}/48$, and

$|a^{\prime}-a|\leq\delta a,\quad|b^{\prime}-b|\leq\delta
a,\quad|c^{\prime}-c|\leq\delta c,$

then $a^{\prime},b^{\prime},c^{\prime}$ form a Euclidean triangle with
opposite inner angles $A^{\prime},B^{\prime},C^{\prime}$ respectively such
that

$|A^{\prime}-A|\leq\frac{24}{\epsilon}\delta,$

and

$\bigg{|}|\triangle A^{\prime}B^{\prime}C^{\prime}|-|\triangle
ABC|\bigg{|}\leq\frac{576}{\epsilon^{2}}\delta\cdot|\triangle ABC|.$

By the following lemma we show that the linear map between two Euclidean
triangles is close to isometry if their corresponding edge lengths are close.

######  Lemma 4.2.

Assume $\triangle ABC$ ($\triangle A^{\prime}B^{\prime}C^{\prime}$) is a
Euclidean triangle with edge lengths $a,b,c$ (resp.
$a^{\prime},b^{\prime},c^{\prime}$), and all the angles in $\triangle ABC$ are
at least $\epsilon>0$, and $\delta<\epsilon^{2}/576$, and

$|a^{\prime}-a|\leq\delta a,\quad|b^{\prime}-b|\leq\delta
a,\quad|c^{\prime}-c|\leq\delta c,$

and $\lambda_{1},\lambda_{2}$ are the two singular values of the unique linear
map sending $\triangle ABC$ to $\triangle A^{\prime}B^{\prime}C^{\prime}$
preserving the correspondence of the vertices. Then

$1-\frac{10^{4}}{\epsilon^{4}}\delta\leq\lambda_{i}\leq
1+\frac{10^{4}}{\epsilon^{4}}\delta,\quad i=1,2.$

###### Proof.

By Lemma 4.1, $|A-A^{\prime}|,|B-B^{\prime}|,|C-C^{\prime}|$ are all less or
equal to $24\delta/\epsilon<\epsilon/2$, and thus
$A^{\prime},B^{\prime},C^{\prime}$ are all at least $\epsilon/2$. Then again
by Lemma 4.1 it is easy to show that

$a^{\prime 2}\leq\frac{2|\triangle
A^{\prime}B^{\prime}C^{\prime}|}{\sin^{2}(\epsilon/2)}\leq\frac{64|\triangle
ABC|}{\epsilon^{2}}.$

It is well known that

$\lambda_{1}\lambda_{2}=\frac{|\triangle
A^{\prime}B^{\prime}C^{\prime}|}{|\triangle ABC|}$

and thus by Lemma 4.1

(7) $|\lambda_{1}\lambda_{2}-1|<\frac{576}{\epsilon^{2}}\delta.$

In [13] we can find the formula

$\lambda_{1}^{2}+\lambda_{2}^{2}=\frac{a^{\prime 2}\cot A+b^{\prime 2}\cot
B+c^{\prime 2}\cot C}{2|\triangle ABC|}.$

Applying this formula to the special case $\triangle
A^{\prime}B^{\prime}C^{\prime}=\triangle ABC$, we get

$2=\frac{a^{\prime 2}\cot A^{\prime}+b^{\prime 2}\cot B^{\prime}+c^{\prime
2}\cot C^{\prime}}{2|\triangle A^{\prime}B^{\prime}C^{\prime}|},$

, which implies

(8) $\lambda_{1}^{2}+\lambda_{2}^{2}-2\lambda_{1}\lambda_{2}=\frac{a^{\prime
2}(\cot A-\cot A^{\prime})+b^{\prime 2}(\cot B-\cot B^{\prime})+c^{\prime
2}(\cot C-\cot C^{\prime})}{2|\triangle ABC|}.$

Denote $f(x)=\cot x$, then $f^{\prime}(x)=-1/\sin^{2}x$ and
$f^{\prime\prime}(x)=2\cos x/\sin^{3}x$. By Taylor’s expansion, there exists
$\xi_{A}$ between $A$ and $A^{\prime}$ such that

$a^{\prime 2}(\cot A-\cot A^{\prime})=a^{\prime
2}[f^{\prime}(A^{\prime})(A-A^{\prime})+\frac{1}{2}f^{\prime\prime}(\xi_{A})(A-A^{\prime})^{2}]$
$=-\frac{a^{\prime 2}}{\sin^{2}A^{\prime}}(A-A^{\prime})+\frac{a^{\prime
2}}{2}f^{\prime\prime}(\xi_{A})(A-A^{\prime})^{2}=-(2R)^{2}(A-A^{\prime})+\frac{a^{\prime
2}}{2}f^{\prime\prime}(\xi_{A})(A-A^{\prime})^{2}$

where $R$ is the radius of the cicumcircle of $\triangle
A^{\prime}B^{\prime}C^{\prime}$, and

$\left|\frac{a^{\prime
2}}{2}f^{\prime\prime}(\xi_{A})(A-A^{\prime})^{2}\right|\leq\frac{64|\triangle
ABC|}{\epsilon^{2}}\cdot\frac{2}{\sin^{3}(\epsilon/2)}\cdot\left(\frac{24\delta}{\epsilon}\right)^{2}\leq|\triangle
ABC|\cdot\frac{10^{6}\cdot\delta^{2}}{\epsilon^{7}}.$

Combining the similar computation for $B$ and $C$, we get that the right hand
side of equation (8) is less or equal to $3\times
10^{6}\delta^{2}/\epsilon^{7}$, and thus
$|\lambda_{1}-\lambda_{2}|\leq\sqrt{3\times 10^{6}\delta^{2}/\epsilon^{7}}\leq
10^{4}\delta/\epsilon^{4}$. Then by equation (7) and the fact that
$10^{4}\delta/\epsilon^{4}\geq 576\delta/\epsilon^{2}$, it is easy to prove
that

$1-\frac{10^{4}}{\epsilon^{4}}\delta\leq\lambda_{i}\leq
1+\frac{10^{4}}{\epsilon^{4}}\delta.$

∎

We also need the following lemma to compare the angle differences between
spherical and Euclidean triangles of the same edge lengths.

######  Lemma 4.3 (Special case of Lemma 5.3 in [4]).

Assume $\triangle ABC$ (resp. $\triangle A^{\prime}B^{\prime}C^{\prime}$) is a
Euclidean (resp. spherical) triangle with the edge lengths $a,b,c$, and
$diam(\triangle A^{\prime}B^{\prime}C^{\prime})<{\pi}/3$. Then

$|A^{\prime}-A|\leq 2(a+b+c)^{2}.$

###  4.2. Geodesic triangulations of surfaces

The following cubic estimate quantifies the error between the discrete
conformal change and the smooth conformal change.

######  Lemma 4.4 (Lemma 4.3 in [4]).

Suppose $(M,g)$ is a closed Riemannian surface, and $u\in C^{\infty}(M)$ is a
conformal factor. Then there exists $C=C(M,g,u)>0$ such that for any $x,y\in
M$,

$|d_{e^{2u}g}(x,y)-e^{\frac{1}{2}(u(x)+u(y))}d_{g}(x,y)|\leq Cd_{g}(x,y)^{3}.$

The following lemma shows that a smooth conformal change will preserve the
existence of the geodesic triangulation, and its $\epsilon$-regularity.

######  Lemma 4.5 (Part (a) of Lemma 4.4 in [4]).

Suppose $(M,g)$ is a closed Riemannian surface, and $T$ is a geodesic
triangulation of $(M,g)$. Let $l\in\mathbb{R}^{E(T)}$ denote the geodesic
lengths of the edges of $T$, and assume $(T,l)_{E}$ is $\epsilon$-regular.
Given a conformal factor $u\in C^{\infty}(M)$, there exits a constant
$\delta=\delta(M,g,u,\epsilon)$ such that if $|l|\leq\delta$ then there exists
a geodesic triangulation $T^{\prime}$ in $(M,e^{2u}g)$ such that $T^{\prime}$
is homotopic to $T$ relative to $V(T)=V(T^{\prime})$. Further
$(T^{\prime},\bar{l})_{E}$ is $(\epsilon/2)$-regular, where
$\bar{l}\in\mathbb{R}^{E(T^{\prime})}$ denotes the geodesic lengths of the
edges of $T^{\prime}$ in $(M,e^{2u}g)$.

The following lemma shows that an $\epsilon$-regular geodesic triangulation of
a closed surface satisfies the $C$-isoperimetric condition if the edge lengths
are sufficiently small.

######  Lemma 4.6 (Part (b) of Lemma 4.4 in [4]).

Suppose $(M,g)$ is a closed Riemannian surface, and $T$ is a geodesic
triangulation of $(M,g)$ with geodesic length $l$ such that $(T,l)_{E}$ is
$\epsilon$-regular. Then there exists a constant $\delta=\delta(M,g,\epsilon)$
such that if $|l|<\delta$, $(T,l)$ is $C$-isoperimetric for some constant
$C=C(M,g,\epsilon)>0$.

However, what we really need is the following modified version of Lemma 4.6.

######  Lemma 4.7.

Suppose $(M,g)$ is a closed Riemannian surface, and $T$ is a geodesic
triangulation of $(M,g)$ with geodesic length $l$ such that $(T,l)_{E}$ is
$\epsilon$-regular. Assume $v\in V$ and $star(v)\subset V$ contains $v$ and
its neighbors in $T$, and let $\hat{V}=V-star(v)$ and $G=(\hat{V},\hat{E})$ be
the subgraph of $(V(T),E(T))$ generated by $\hat{V}$. Then there exists a
constant $\delta=\delta(M,g,\epsilon)$ such that if $|l|<\delta$,
$(\hat{T},l|_{\hat{E}})$ is $C$-isoperimetric for some constant
$C=C(M,g,\epsilon)>0$.

###### Proof.

By Lemma 4.6, we can find constants $\delta(M,g,\epsilon)>0$ and
$C(M,g,\epsilon)>0$ such that $(T,l)$ is $C$-isoperimetric if $|l|<\delta$.
Now assume $B=\\{i\in\hat{V}:\exists j\in V-\hat{V}\text{ s.t. }ij\in E\\}$ is
the set of boundary vertices of $G$ in $T$, and $V_{0}\subset\hat{V}$, and
$\hat{\partial}V_{0}$ (resp. $\partial V$) is the boundary of $V_{0}$ in $G$
(resp. $T$). We consider the following three cases:

Case 1: $V_{0}\cap B=\emptyset$. Then $|\hat{\partial}V_{0}|_{l}=|\partial
V_{0}|_{l}$ and $|V|_{l}-|V_{0}|_{l}\geq|\hat{V}|_{l}-|V_{0}|_{l}$. Since
$(T,l)$ is $C$-isoperimetric, we have

$C|\hat{\partial}V_{0}|_{l}^{2}\geq\min\\{|V_{0}|_{l},|\hat{V}|_{l}-|V_{0}|_{l}\\}.$

Case 2: $B\subset V_{0}$. In this case,
$\hat{\partial}V_{0}=\partial(V_{0}\cup star(v))$. Since $(T,l)$ is
$C$-isoperimetric, we have

$C|\hat{\partial}V_{0}|_{l}^{2}=C|\partial(V_{0}\cup
star(v))|^{2}_{l}\geq\min\\{|V_{0}\cup star(v)|_{l},|V|_{l}-|V_{0}\cup
star(v)|_{l}\\}$

Clearly, $|V_{0}\cup star(v)|_{l}\geq|V_{0}|_{l}$, and $|V|_{l}-|V_{0}\cup
star(v)|_{l}=|\hat{V}|_{l}-|V_{0}|_{l}$. Then

$C|\hat{\partial}V_{0}|_{l}^{2}\geq\min\\{|V_{0}|_{l},|\hat{V}|_{l}-|V_{0}|_{l}\\}.$

Case 3: $V_{0}\cap B\neq\emptyset$ and $B\not\subset V_{0}$. It is not
difficult to show that $B$ is connected in $V$ since the $1$-skeleton of $T$
is $4$-vertex-connected, so there is an edge $ij\in\hat{\partial}V_{0}$ such
that $i\in B\cap V_{0}$ and $j\in B-V_{0}$. By the $\epsilon$-regularity, the
degree of each vertex in $T$ is bounded by $\lfloor 2\pi/\epsilon\rfloor$ if
$\delta(M,g,\epsilon)$ is sufficiently small, and the ratio of the two edge
lengths in a triangle of $T$ is at least $\sin\epsilon$. So there is a
constant $C_{1}(M,g,\epsilon)>0$ such that

$C_{1}l_{ij}\geq\sum_{xy\in E(T)-\hat{E}}l_{xy}\geq|\partial
V_{0}|_{l}-|\hat{\partial}V_{0}|_{l}.$

Then

$C(1+C_{1})^{2}|\hat{\partial}V_{0}|_{l}^{2}\geq C|\partial
V_{0}|_{l}^{2}\geq\min\\{|V_{0}|_{l},|\hat{V}|_{l}-|V_{0}|_{l}\\}.$

∎

##  5\. Proof of the Main Theorem 1.7

Assume $\epsilon>0$ is a fixed constant and $(T,l)_{S}$ is $\epsilon$-regular
and $|l|<\delta$ where

$\delta=\delta(M,g,X,Y,Z,\epsilon)>0$

is a sufficiently small constant to be determined. By Lemma 4.3 we may assume
that $(T,l)_{E}$ is $(\epsilon/2)$-regular. By Lemma 4.5, we may assume that
there exists a geodesic triangulation $T^{\prime}$ of $(M,e^{2\bar{u}}g)$ such
that $T^{\prime}$ is homotopic to $T$ relative to $V(T)=V(T^{\prime})$. Denote
$\bar{l}\in\mathbb{R}^{E(T)}\cong\mathbb{R}^{E(T^{\prime})}$ as the geodesic
edge length of $T^{\prime}$ in $(M,e^{2\bar{u}}g)$, and then $(T,\bar{l})_{S}$
is isometric to the unit sphere $(M,e^{2\bar{u}}g)$ and has zero discrete
curvatures. Again by Lemma 4.5 we may assume that $(T,\bar{l})_{E}$ is
$(\epsilon/4)$-regular. Then by Lemma 4.3 we may assume $(T,\bar{l})_{S}$ is
$(\epsilon/5)$-regular, and thus is strictly Delaunay, and then by Proposition
1.5 $p(T_{P})=\phi(T^{\prime})$ where $P\in\mathcal{P}$ is the boundary of the
convex hull of $\phi(V(T))$. By Lemma 3.1, $Q=p_{N}(\mathring{P})$ is a convex
polygon, and $T_{Q}=p_{N}(\mathring{T}_{P})$ is a geodesic triangulation of
$Q$. Denote $l_{Q}\in\mathbb{R}^{E(T_{Q})}$ as the edge lengths in $Q$, and
then by Lemma 3.2 there exists a constant
$\epsilon^{\prime}(M,g,X,Y,Z,\epsilon)>0$ such that $(T_{Q},l_{Q})_{E}$ is
$\epsilon^{\prime}$-regular and $K_{i}\geq\epsilon^{\prime}$ for any boundary
vertex in $(T_{Q},l_{Q})_{E}$. The combinatorial structures of $T$,
$T^{\prime}$, $T_{P}$ and $p(T_{P})=\phi(T^{\prime})$ are naturally
identified. We also identify the combinatorial structures of
$\mathring{T}_{P}$ and $T_{Q}$ and just denote it as $\mathring{T}$. Denote
$l_{P}=2\sin(\bar{l}/2)\in\mathbb{R}^{E(T)}$ as the edge length of $T_{P}$ on
$P$, and then by equation (6) on $\mathring{T}$ we have

$l_{Q}=w*l_{P}\quad\text{where}\quad
w_{i}=\log\frac{2}{|\phi(i)-N|_{2}^{2}}=\log\frac{|p_{N}(\phi(i))|_{2}^{2}+1}{2},\quad\forall
i\in V(\mathring{T}).$

Denote $l_{P}^{\prime}=\bar{u}*2\sin(l/2)\in\mathbb{R}^{E(T)}$ and
$l_{Q}^{\prime}=w*l_{P}^{\prime}=(\bar{u}+w)*2\sin(l/2)\in\mathbb{R}^{E(\mathring{T})}$,
and $K(u)\in\mathbb{R}^{V(\mathring{T})}$ as the discrete curvature in
$(\mathring{T},u*2\sin(l/2))_{E}$.

In the following proof, for simplicity we will use the notation $a=O(b)$ to
represent that if $\delta(M,g,X,Y,Z,\epsilon)$ is sufficiently small, then
$|a|\leq Cb$ for some constant $C(M,g,X,Y,Z,\epsilon)>0$. We summarize the
remaining part of the proof in three steps:

  1. (a)

Estimate the curvature $K(\bar{u}+w)$ of $(\mathring{T},l_{Q}^{\prime})_{E}$
for interior vertices.

  2. (b)

Construct a flow $u(t):[0,1]\to\mathbb{R}^{V(\mathring{T})}$, starting from
$u(0)=\bar{u}+w$, such that $u(t)$ linearly eliminates the discrete curvature
$K(u)$ for interior vertices, i.e.,

$K_{i}(u(t))=(1-t)K_{i}(\bar{u}+w)$

for any interior vertex $i$ of $\mathring{T}$. Furthermore, we will also show
that $|u^{\prime}(t)|=O(|l|)$, and $(\mathring{T},u(1)*2\sin(l/2))_{E}$ is
isometric to a convex polygon in the plane.

  3. (c)

After a proper normalization, which is a small perturbation, we use the
inverse of the stereographic projection to construct the desired polyhedral
surface $P\in\mathcal{P}$.

###  5.1. Step 1: The estimate of curvatures

By Lemma 4.4,

$|\bar{l}_{ij}-(\bar{u}*l)_{ij}|=O(l_{ij}^{3}).$

Notice the fact that $|x-2\sin(x/2)|\leq 10x^{3}$ if $|x|<0.01$, so

$|(l_{P})_{ij}-(l_{P}^{\prime})_{ij}|=O(l_{ij}^{3}),$

and

(9)
$\Bigg{|}\frac{(l_{Q}^{\prime})_{ij}-(l_{Q})_{ij}}{(l_{Q})_{ij}}\Bigg{|}=\Bigg{|}\frac{(l_{P}^{\prime})_{ij}-(l_{P})_{ij}}{(l_{P})_{ij}}\Bigg{|}=O(l_{ij}^{2}).$

Given a triangle $\triangle ijk\in F(\mathring{T})$, denote
$\theta^{i}_{jk}(u)$ (resp. $\bar{\theta}^{i}_{jk}$) as the inner angle at $i$
in $\triangle ijk$ in $(\mathring{T},u*2\sin(l/2))_{E}$ (resp.
$(\mathring{T},{l_{Q}})_{E}$), and $K_{i}(u)$ as the discrete curvature at $i$
in $\triangle ijk$ in $(\mathring{T},u*2\sin(l/2))_{E}$.

Since $(\mathring{T},l_{Q})_{E}$ is $\epsilon^{\prime}$-regular, by equation
(9) and Lemma 4.1,

$\alpha^{i}_{jk}:=\bar{\theta}^{i}_{jk}-\theta^{i}_{jk}(\bar{u}+w)=O(l_{ij}^{2}).$

So for sufficiently small $\delta(M,g,\epsilon)$, we have

$|\alpha^{i}_{jk}|\leq\frac{\epsilon^{\prime}}{4}.$

Then $(\mathring{T},{l_{Q}^{\prime}})_{E}$ is $(\epsilon^{\prime}/2)$-regular.
Since $(\mathring{T},{l_{Q}})_{E}$ is globally flat, for any
$i\in\text{int}(\mathring{T})$,

$\sum_{ijk\in{F}}\bar{\theta}^{i}_{jk}=2\pi.$

So

$K_{i}(\bar{u}+w)=2\pi-\sum_{ijk\in F}\theta^{i}_{jk}(\bar{u}+w)=\sum_{ijk\in
F}(\bar{\theta}^{i}_{jk}-\theta^{i}_{jk}(\bar{u}+w))=\sum_{ijk\in
F}\alpha^{i}_{jk}\quad\text{if $i\in\text{int}(\mathring{T})$}.$

Set $x\in\mathbb{R}^{E(\mathring{T})}_{A}$ be such that if $ij$ is an interior
edge of $\mathring{T}$

$x_{ij}=\frac{\alpha^{i}_{jk}-\alpha^{j}_{ik}}{3}+\frac{\alpha^{i}_{jk^{\prime}}-\alpha^{j}_{ik^{\prime}}}{3},$

where $\triangle ijk$ and $\triangle ijk^{\prime}$ are adjacent triangles in
$\mathring{T}$, and if $ij$ is a boundary edge of $\mathring{T}$

$x_{ij}=\frac{\alpha^{i}_{jk}-\alpha^{j}_{ik}}{3},$

where $\triangle ijk$ is a triangle in $\mathring{T}$. Then

(10) $x_{ij}=O(l_{ij}^{2}),$

and it is straightforward to verify that for any
$i\in\text{int}(\mathring{T})$

$div(x)_{i}=\sum_{j:j\sim i}x_{ij}=\sum_{ijk\in
F}\alpha_{jk}^{i}=K_{i}(\bar{u}+w)$

using $\alpha^{i}_{jk}+\alpha^{j}_{ik}+\alpha^{k}_{ij}=0$.

###  5.2. Step 2: The construction of the flow

Consider the sets defined by

$\tilde{\Omega}=\\{u\in\mathbb{R}^{V(\mathring{T})}:\text{$(\mathring{T},u*2\sin\frac{l}{2})_{E}$
satisfies the triangle inequality and is strictly Delaunay}\\},$

and

$\Omega=\\{u\in\tilde{\Omega}:\text{ $(\mathring{T},u*2\sin\frac{l}{2})_{E}$
is $\frac{\epsilon^{\prime}}{4}$-regular},|u-(\bar{u}+w)|\leq 1\\}.$

Notice that $\tilde{\Omega}$ is an open domain in
$\mathbb{R}^{V(\mathring{T})}$ and $\Omega$ is a compact subset of
$\tilde{\Omega}$. By the construction, $(\bar{u}+w)$ is in the interior of
$\Omega$, since $(\mathring{T},{l_{Q}^{\prime}})_{E}$ is
$(\epsilon^{\prime}/2)$-regular. Given $u\in\tilde{\Omega}$ and an interior
edge $ij$ in $\mathring{T}$, denote

$\eta_{ij}(u)=\frac{1}{2}(\cot\theta^{k}_{ij}(u)+\cot\theta^{k^{\prime}}_{ij}(u))$

where $\triangle ijk$ and $\triangle ijk^{\prime}$ are adjacent triangles in
$\mathring{T}$. Then for $u\in\Omega$,

(11)
$2\eta_{ij}(u)=\cot\theta^{k}_{ij}(u)+\cot\theta^{k^{\prime}}_{ij}(u)=\frac{\sin(\theta^{k}_{ij}(u)+\theta^{k^{\prime}}_{ij}(u))}{\sin\theta^{k}_{ij}(u)\sin\theta^{k^{\prime}}_{ij}(u)}\geq\sin(\theta^{k}_{ij}(u)+\theta^{k^{\prime}}_{ij}(u))\geq\sin\frac{\epsilon^{\prime}}{4}$

for any interior edge $ij$ in $\mathring{T}$.

Consider the following system of differential equations on $\tilde{\Omega}$,

(12) $\displaystyle\frac{\partial K_{i}}{\partial u}\frac{du}{dt}$
$\displaystyle=-K_{i}(\bar{u}+w)=-div(x)_{i},\quad\quad$ $\displaystyle\text{
$i\in\text{int}(\mathring{T})$},$ $\displaystyle\frac{du_{i}}{dt}$
$\displaystyle=(\log
2-2\log(2\sin\frac{l_{iX}}{2})-\bar{u}_{X})-(\bar{u}_{i}+w_{i}),\quad\quad$
$\displaystyle\quad{i\in\text{bdy}(\mathring{T})},$ $\displaystyle u(0)$
$\displaystyle=\bar{u}+w,$

where $\bar{u}_{X}$ is the value of $\bar{u}$ at the marked point $X$ sent to
the north pole, and $l_{iX}$ is the length of the edge $iX$ given by $l$. We
want to show that the solution $u(t)$ exists on $[0,1]$, then it is easy to
see that $K_{i}(u(1))=0$ for an interior vertex $i$ of $\mathring{T}$, and

$u_{i}(1)=\log 2-2\log(2\sin\frac{l_{iX}}{2})-\bar{u}_{X}$

for a boundary vertex $i$ of $\mathring{T}$.

For a boundary vertex $i$ of $\mathring{T}$, $u_{i}(t)$ can be easily solved
as

$u_{i}(t)=t(\log
2-2\log(2\sin\frac{l_{iX}}{2})-\bar{u}_{X})+(1-t)(\bar{u}_{i}+w_{i}),$

and

$\displaystyle\begin{split}\frac{du_{i}}{dt}=&(\log
2-2\log(2\sin\frac{l_{iX}}{2})-\bar{u}_{X})-(\bar{u}_{i}+w_{i})\\\ =&\log
2-2\log(2\sin\frac{l_{iX}}{2})-\bar{u}_{X}-\bar{u}_{i}-\log\frac{2}{(l_{P})_{iX}^{2}}\\\
=&-{2\log(2\sin\frac{l_{iX}}{2})}-\bar{u}_{X}-\bar{u}_{i}+2\log{(l_{P})_{iX}}\\\
=&-2\log(l_{P}^{\prime})_{iX}+2\log(l_{P})_{iX}\\\
=&O(l_{iX}^{2}).\end{split}$

Now let us focus on solving $u_{i}(t)$ for all the interior vertices of
$\mathring{T}$. Let $\hat{V}$ be the set of interior vertices of
$\mathring{T}$, and $G=(\hat{V},\hat{E})$ be the subgraph of
$(V(\mathring{T}),E(\mathring{T}))$ generated by $\hat{V}$. It is easy to show
that $G$ is nonempty and connected. Let $\hat{u}\in\mathbb{R}^{\hat{V}}$ and
$\hat{x}\in\mathbb{R}^{\hat{E}}_{A}$ and $\hat{\eta}\in\mathbb{R}^{\hat{E}}$
be the restrictions of $u$ and $x$ and $\eta$ respectively on
$G=(\hat{V},\hat{E})$, and $\hat{\Delta}=\Delta_{\hat{\eta}}$ be the
associated discrete Laplacian on $G$. Then by Proposition 2.2 it is
straightforward to verify that Equation (12) can be rewritten as

(13) $(D-\hat{\Delta})\frac{d\hat{u}}{dt}=-div(\hat{x})+y,$

where

(a) $D\in\mathbb{R}^{\hat{V}\times\hat{V}}$ is a nonzero diagonal matrix and

$D_{ii}=\sum_{j\sim i:j\notin\hat{V}}\eta_{ij}\geq 0,$

and

(b)

$y_{i}=\sum_{j\sim i:j\notin\hat{V}}\eta_{ij}\frac{du_{j}}{dt}-\sum_{j\sim
i:j\notin\hat{V}}x_{ij}.$

For $u\in\tilde{\Omega}$, it is easy to show that $(D-\hat{\Delta})$ is
positive definite, by the fact that $G$ is connected, $\eta_{ij}>0$ for any
$ij\in\hat{E}$, and $D$ is nonzero and non-negative. So equation (13) locally
has a unique solution in $\tilde{\Omega}$.

Assume the maximum existing open interval for the solution
$\hat{u}(t)\in\Omega$ is $[0,T_{0})$ where $0<T_{0}\leq+\infty$. For
$t\in[0,T_{0})$, we have

$\frac{d\hat{u}}{dt}(t)=O(|l|\cdot|\hat{V}|_{l}^{1/2})$

by Lemma 2.1, and Lemma 4.7, and equation (11) and (10), and the fact that

$y_{i}=\sum_{j\sim
i:j\notin\hat{V}}(\eta_{ij}\frac{du_{j}}{dt}-x_{ij})=O\left(\sum_{j\sim
i:j\notin\hat{V}}(\eta_{ij}l_{jX}^{2}+l_{ij}^{2})\right)=O(D_{ii}|l|^{2})=O(D_{ii}|l|\cdot|\hat{V}|_{l}^{1/2}).$

Further

$\displaystyle\begin{split}&|\hat{V}|_{l}\leq|V|_{l}=\sum_{ij\in
E}l_{ij}^{2}=O(\sum_{ij\in E}\bar{l}_{ij}^{2})=O(\sum_{ijk\in
F}(\bar{l}_{ij}^{2}+\bar{l}_{jk}^{2}+\bar{l}_{ik}^{2}))\\\ =&O(\sum_{ijk\in
F}Area(\triangle
ijk,\bar{l})_{S})=O(Area((T,\bar{l})_{S}))=O(Area(\mathbb{S}^{2}))=O(1),\end{split}$

and thus $(du/dt)(t)=O(|l|)$ for $t\in[0,T_{0})$.

If $T_{0}\leq 1$, combining Lemma 4.1, we have

$|u(T_{0})-(\bar{u}+w)|=O(|l|),\quad\text{and
}\quad|\theta_{jk}^{i}(u(T_{0}))-\theta_{jk}^{i}(\bar{u}+w)|=O(|l|).$

This implies that $u(T_{0})\in\text{int}(\Omega)$ if $\delta$ is sufficiently
small, which contradicts to the maximality of $T_{0}$. Thus, $T_{0}>1$ and
$u(1)$ is well-defined. Further we have that

(a) $K_{i}(u(1))=0$ for any interior vertex $i$ of $\mathring{T}$, and

(b) $u_{i}(1)=\log 2-2\log(2\sin\frac{l_{iX}}{2})-\bar{u}_{X}$ for any
boundary vertex $i$ of $\mathring{T}$, and

(c) $u(1)-(\bar{u}+w)=O(|l|)$, and

(d) $(\mathring{T},u(1)*2\sin\frac{l}{2})_{E}$ is strictly Delaunay, and

(e) $K_{i}(u(1))>0$ for any boundary vertex $i$ in $\mathring{T}$.

###  5.3. Step 3: The normalization and the inverse of the stereographic
projection

We know that $(\mathring{T},u(1)*2\sin\frac{l}{2})_{E}$ is isometric to a
closed convex polygon in $\mathbb{C}$. Let $f$ be the piecewise linear map
from $(\mathring{T},u(1)*2\sin\frac{l}{2})_{E}$ to $(\mathring{T},l_{Q})_{E}$
that preserves the triangulation and is linear on each triangle. From equation
(9) and the fact that $u(1)-(\bar{u}+w)=O(|l|)$, we can deduce that

$\Big{|}\frac{(u(1)*2\sin\frac{l}{2})_{ij}-({l_{Q}})_{ij}}{({l_{Q}})_{ij}}\Big{|}=\Big{|}\frac{(u(1)*2\sin\frac{l}{2})_{ij}-((\bar{u}+w)*2\sin\frac{l}{2})_{ij}}{({l_{Q}})_{ij}}\Big{|}+O(|l|^{2})=O(|l|).$

Then by Lemma 4.2, $\|Df\|_{2}$ and $\|Df^{-1}\|_{2}$ are both
$(1+C|l|)$-Lipschitz for some constant $C(M,g,X,Y,Z,\epsilon)>0$. So the
distance $d_{YZ}$ between $Y$ and $Z$ in $(\mathring{T},u(1)*2\sin(l/2))_{E}$
lies in $[1-C|l|,1+C|l|]$. So we can scale
$(\mathring{T},u(1)*2\sin(l/2))_{E}$ by letting $\tilde{u}=u(1)-\log d_{YZ}$,
and then $(\mathring{T},\tilde{u}*2\sin(l/2))_{E}$ is still isometric to a
convex polygon and the distance between $Y$ and $Z$ is $1$, and

$\Big{|}\frac{\tilde{u}*2\sin\frac{l}{2}-{(l_{Q})}_{ij}}{{(l_{Q})}_{ij}}\Big{|}=\Big{|}\frac{\tilde{u}*2\sin\frac{l}{2}-u(1)*2\sin\frac{l}{2}}{{(l_{Q})}_{ij}}\Big{|}+O(|l|)=O(|l|).$

Let $g$ be the isometry from $(\mathring{T},\tilde{u}*2\sin\frac{l}{2})_{E}$
to a closed convex polygon $Q_{1}$ in $\mathbb{C}$ such that $g(Z)=0$ and
$g(Y)=1$. Then for any $i\in\mathring{V}$, the above bi-Lipschitz property of
$f$ implies that

$\left|\log\frac{|g(i)|_{2}}{|p_{N}(\phi(i))|_{2}}\right|=O(|l|).$

Now we are ready to project the points in the plane back to the sphere. Let

${V_{1}}=(p_{N}|_{\mathbb{S}^{2}})^{-1}(g({V(\mathring{T})}))\cup\\{N\\}$

and ${P_{1}}$ be the convex hull of ${V_{1}}$. Then by part (b) of Lemma 3.1,
$P_{1}\in\mathcal{P}$ and $p_{N}(\mathring{P}_{1})=Q_{1}$ and
$p_{N}(\mathring{T}_{P_{1}})=g(\mathring{T})$. Naturally identify the
combinatorial structures of $T$ and $T_{P_{1}}$, and denote
$l_{P_{1}}\in\mathbb{R}^{E(T)}$ as the edge length on $P_{1}$. We will verify
that

$l_{P_{1}}=u*2\sin\frac{l}{2}$

where $u_{X}=\bar{u}_{X}+\log d_{YZ}$ and

$u_{i}=\tilde{u}_{i}-w_{i}^{\prime},\quad\text{where}\quad
w_{i}^{\prime}=\log\frac{|g(i)|^{2}+1}{2}$

if $i\in V(\mathring{T})$. If $ij\in E(\mathring{T})$,

$(l_{P_{1}})_{ij}=(u*2\sin\frac{l}{2})_{ij}$

by implementing Lemma 3.1 on $l_{P_{1}}$ and $l_{Q_{1}}$. For edge $iX\in
E(T)$, we have that

$\displaystyle\begin{split}&\log(u*2\sin\frac{l}{2})_{iX}\\\
=&\log(2\sin\frac{l_{iX}}{2})+{\frac{1}{2}(\bar{u}_{X}+\log d_{YZ}+\log
2-2\log(2\sin\frac{l_{iX}}{2})-\bar{u}_{X}-\log d_{YZ}}-w_{i}^{\prime})\\\
=&{\frac{1}{2}(\log
2-w_{i}^{\prime}})=\frac{1}{2}\log\frac{4}{|g(i)|^{2}+1}=\frac{1}{2}\log(l_{P_{1}})_{iX}^{2}=\log(l_{P_{1}})_{iX}.\end{split}$

So $u$ is our desired discrete conformal factor. As we mentioned in Remark
1.4, such $u$ is known to be unique. It remains to show
$u_{i}-\bar{u}_{i}=O(|l|)$ for any $i\in V$. Notice that

$|w_{i}^{\prime}-w_{i}|=\Big{|}\log\frac{|g(i)|_{2}^{2}+1}{2}-\log\frac{|p_{N}(\phi(i))|_{2}^{2}+1}{2}\Big{|}=\Big{|}\log\frac{|g(i)|_{2}^{2}+1}{|p_{N}(\phi(i)|_{2}^{2}+1}\Big{|}=O(|l|).$

So restricted on $V(\mathring{T})$, we have that

$u-\bar{u}=(\tilde{u}-w^{\prime})-\bar{u}=u(1)-\log
d_{YZ}-w^{\prime}-\bar{u}=(u(1)-\bar{u}-w)+(w-w^{\prime})-\log d_{YZ}=O(|l|).$

On vertex $X$ we have that $u_{X}-\bar{u}_{X}=\log d_{YZ}=O(|l|)$.

## References

  * [1] Feng Luo.  Combinatorial yamabe flow on surfaces.  Communications in Contemporary Mathematics, 6(05):765–780, 2004\.
  * [2] Alexander I Bobenko, Ulrich Pinkall, and Boris A Springborn.  Discrete conformal maps and ideal hyperbolic polyhedra.  Geometry & Topology, 19(4):2155–2215, 2015.
  * [3] David Gu, Feng Luo, and Tianqi Wu.  Convergence of discrete conformal geometry and computation of uniformization maps.  Asian Journal of Mathematics, 23(1):21–34, 2019.
  * [4] Tianqi Wu and Xiaoping Zhu.  The convergence of discrete uniformizations for closed surfaces, 2020\.
  * [5] Burt Rodin and Dennis Sullivan.  The convergence of circle packings to the riemann mapping.  Journal of Differential Geometry, 26(2):349–360, 1987.
  * [6] Feng Luo, Jian Sun, and Tianqi Wu.  Discrete conformal geometry of polyhedral surfaces and its convergence.  arXiv preprint arXiv:2009.12706, 2020.
  * [7] Ulrike Bücking.  $c^{\infty}$-convergence of conformal mappings for conformally equivalent triangular lattices.  Results in Mathematics, 73(2):1–21, 2018.
  * [8] Ulrike Bücking.  Approximation of conformal mappings using conformally equivalent triangular lattices.  In Advances in Discrete Differential Geometry, pages 133–149. Springer, Berlin, Heidelberg, 2016.
  * [9] Boris Springborn.  Ideal hyperbolic polyhedra and discrete uniformization.  Discrete & Computational Geometry, pages 1–46, 2019.
  * [10] Igor Rivin.  Intrinsic geometry of convex ideal polyhedra in hyperbolic 3-space.  Lecture Notes in Pure and Applied Mathematics, pages 275–291, 1994\.
  * [11] Boris Springborn, Peter Schröder, and Ulrich Pinkall.  Conformal equivalence of triangle meshes.  ACM Transactions on Graphics, 27(3):Art–No, 2008.
  * [12] Satyan L Devadoss and Joseph O’Rourke.  Discrete and computational geometry.  Princeton University Press, 2011.
  * [13] Ulrich Pinkall and Konrad Polthier.  Computing discrete minimal surfaces and their conjugates.  Experimental mathematics, 2(1):15–36, 1993.


P. Gijsbers and J. Vanschoren

11institutetext: Eindhoven University of Technology, Eindhoven, Netherlands
11email<EMAIL_ADDRESS>
# GAMA: a General Automated Machine learning Assistant

Pieter Gijsbers[✉]  Joaquin Vanschoren

###### Abstract

The General Automated Machine learning Assistant (GAMA) is a modular AutoML
system developed to empower users to track and control how AutoML algorithms
search for optimal machine learning pipelines, and facilitate AutoML research
itself. In contrast to current, often black-box systems, GAMA allows users to
plug in different AutoML and post-processing techniques, logs and visualizes
the search process, and supports easy benchmarking. It currently features
three AutoML search algorithms, two model post-processing steps, and is
designed to allow for more components to be added.

##  1 Introduction

Automated Machine Learning (AutoML) aims to automate the process of building
machine learning models, for instance, by automating the selection and tuning
of preprocessing and learning algorithms in machine learning pipelines. In
recent years, many AutoML systems have been developed, such as Auto-WEKA [10],
auto-sklearn [4], TPOT [9] and ML-Plan [8]. They vary in the types of
pipelines they build (e.g. fixed or variable length), how they optimize them
(e.g. using evolutionary or Bayesian optimization), and whether or how they
employ meta-learning (e.g. warm-starting) or post-processing (e.g.
ensembling).

We demonstrate111A video demonstration can be found at
https://youtu.be/angsGMvEd1w a new open-source AutoML system, GAMA222Code and
documentation can be found at https://github.com/PGijsbers/gama/ [6], which
distinguishes itself by it modularity (allowing users to compose AutoML
systems from sub-components), extensibility (allowing new components to be
added), transparency (tracking and visualizing the search process to better
understand what the AutoML system is doing), and support for research, such as
integration with the AutoML benchmark [5]. The main difference to our earlier
publication ([6]) is the redesign to allow for a modular AutoML pipeline and
the addition of a graphical user interface.

As such, it caters to a wide range of users, from people without a deep
machine learning background who want an easy-to-use AutoML tool, to those who
want better control and understanding of the AutoML process, and especially
researchers who want to perform systematic AutoML research.

Currently, three different search algorithms and two post-processing
techniques are available, but we welcome and plan to include more techniques
in the future. For novice users, GAMA offers a default configuration shown to
perform well in our benchmarks.

##  2 System Overview

### Modular AutoML Pipeline

Rather than prescribing a specific combination of AutoML techniques, GAMA
allows users to combine different search and post-processing algorithms into a
flexible AutoML ‘pipeline’ that can be tuned to the problem at hand.

There are three optimization algorithms currently implemented in GAMA to
search for optimal machine learning pipelines: random search [1], an
asynchronous successive halving algorithm (ASHA) [7] which uses low-fidelity
estimates to filter out bad pipelines early, and an asynchronous multi-
objective evolutionary algorithm.

After the pipeline search has completed, a post-processing technique will be
executed to construct the final model. It is currently possible to either
train the single best pipeline or create an ensemble out of pipelines
evaluated during search, as described in [2]. In subsequent work, we plan to
expand the number of search and post-processing techniques available out-of-
the-box.

Listing 1 shows how to configure GAMA with non-default search and
postprocessing methods and use it as a drop-in replacement for scikit-learn
estimators.333An always up-to-date version of this listing can be found at
https://pgijsbers.github.io/gama/master/citing.html

New AutoML algorithms or variations to existing ones can be included and
tested with relative ease. For instance, each of the search algorithms
described above has been implemented and integrated in GAMA with less than 170
lines of code, and they can all make use of shared functions for logging,
parallel pipeline evaluation and adhering to runtime constraints. It also
allows users to research other questions, such as how to choose the search
algorithm for AutoML.

Listing 1: Configuring an AutoML pipeline with GAMA

⬇

from gama import GamaClassifier

from gama.search_methods import AsynchronousSuccessiveHalving

from gama.postprocessing import EnsemblePostProcessing

automl = GamaClassifier(

search=AsynchronousSuccessiveHalving(),

post_processing=EnsemblePostProcessing()

)

automl.fit(X, y)

automl.predict(X_test)

automl.fit(X_test, y_test)

### Interface

GAMA comes with a graphical web interface which allows novice users to start
and configure GAMA.Moreover, it visualizes the AutoML process to enable
researchers to easily monitor and analyse the behavior of specific AutoML
configurations.

GAMA logs the creation and evaluation of each pipeline, including meta-data
such as creation time and evaluation duration. For pipelines created through
evolution, it also records the parent pipelines and how they differ. One can
also compare multiple logs at once, creating figures such as Figure 1 that
shows the convergence rate of five different GAMA runs over time on the
airline dataset444https://www.openml.org/d/1169.

Figure 1: Visualization of logs

### Benchmarking

GAMA in integrated with the open-source AutoML Benchmark introduced in [5].
Figure 2 shows the results of running GAMA with its default settings some of
the biggest and most challenging datasets for which each other framework had
results in the original work.555Although we could not run these experiments on
the same (AWS) hardware, we took care to use the same computational
constraints. The full and latest results will be made available in the GAMA
documentation.

Figure 2: Performance benchmark results

##  3 Related work

GAMA compares most closely to auto-sklearn and TPOT as they also optimize
scikit-learn pipelines. Auto-sklearn and GAMA both implement the same
ensembling technique [2]. GAMA and TPOT both feature evolutionary search with
NSGA2 selection [3], although GAMA’s implementation uses asynchronous
evolution, which is often faster. While TPOT and auto-sklearn have a fixed
AutoML pipeline, they do allow modifications to their search space. To the
best of our knowledge, GAMA is the only AutoML framework that offers a modular
and extensible composition of AutoML systems, and extensive support for AutoML
research.

##  4 Conclusion

In this proposal we presented GAMA, an open-source AutoML tool which
facilitates AutoML research and skillful use through its modular design and
built-in logging and visualization. Novice users can make use of the graphical
interface to start GAMA, or simply use the default configuration which is
shown to generate models of similar performance to other AutoML frameworks.
Researchers can leverage GAMA’s modularity to integrate and test new AutoML
search procedures in combination with other readily available building blocks,
and then log, visualize, and analyze their behavior, or run extensive
benchmarks. In future work, we aim to integrate additional search techniques
as well as extend the AutoML pipeline with additional steps, such as warm-
starting the pipeline search with meta-data.

## Acknowledgements

This software was developed with support from the Data Driven Discovery of
Models (D3M) program run by DARPA and the Air Force Research Laboratory.

## References

  * [1] Bergstra, J., Bengio, Y.: Random search for hyper-parameter optimization. Journal of Machine Learning Research 13(Feb), 281–305 (2012)
  * [2] Caruana, R., Niculescu-Mizil, A., Crew, G., Ksikes, A.: Ensemble selection from libraries of models. In: Proceedings of the twenty-first international conference on Machine learning. p. 18 (2004)
  * [3] Deb, K., Pratap, A., Agarwal, S., Meyarivan, T.: A fast and elitist multiobjective genetic algorithm: Nsga-ii. IEEE transactions on evolutionary computation 6(2), 182–197 (2002)
  * [4] Feurer, M., Klein, A., Eggensperger, K., Springenberg, J., Blum, M., Hutter, F.: Efficient and robust automated machine learning. In: Advances in neural information processing systems. pp. 2962–2970 (2015)
  * [5] Gijsbers, P., LeDell, E., Thomas, J., Poirier, S., Bischl, B., Vanschoren, J.: An open source automl benchmark. arXiv preprint arXiv:1907.00909 (2019)
  * [6] Gijsbers, P., Vanschoren, J.: Gama: genetic automated machine learning assistant. Journal of Open Source Software 4(33), 1132 (2019)
  * [7] Li, L., Jamieson, K., Rostamizadeh, A., Gonina, E., Hardt, M., Recht, B., Talwalkar, A.: Massively parallel hyperparameter tuning. arXiv preprint arXiv:1810.05934 (2018)
  * [8] Mohr, F., Wever, M., Hüllermeier, E.: Ml-plan: Automated machine learning via hierarchical planning. Machine Learning 107(8-10), 1495–1515 (2018)
  * [9] Olson, R.S., Urbanowicz, R.J., Andrews, P.C., Lavender, N.A., Kidd, L.C., Moore, J.H.: Applications of Evolutionary Computation: 19th European Conference, EvoApplications 2016, Porto, Portugal, March 30 – April 1, 2016, Proceedings, Part I, chap. Automating Biomedical Data Science Through Tree-Based Pipeline Optimization, pp. 123–137. Springer International Publishing (2016)
  * [10] Thornton, C., Hutter, F., Hoos, H.H., Leyton-Brown, K.: Auto-WEKA: Combined selection and hyperparameter optimization of classification algorithms. In: Proc. of KDD-2013. pp. 847–855 (2013)


# An Empirical Study on Internet Traffic Prediction Using Statistical Rolling
Model

Sajal Saha, Anwar Haque, and *Greg Sidebottom  Department of Computer Science
University of Western Ontario, London, ON, Canada
*Juniper Networks, Kanata, ON, Canada
Email:{ssaha59<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Real-world IP network traffic is susceptible to external and internal factors
such as new internet service integration, traffic migration, internet
application, etc. Due to these factors, the actual internet traffic is non-
linear and challenging to analyze using a statistical model for future
prediction. In this paper, we investigated and evaluated the performance of
different statistical prediction models for real IP network traffic; and
showed a significant improvement in prediction using the rolling prediction
technique. Initially, a set of best hyper-parameters for the corresponding
prediction model is identified by analyzing the traffic characteristics and
implementing a grid search algorithm based on the minimum Akaike Information
Criterion (AIC). Then, we performed a comparative performance analysis among
AutoRegressive Integrated Moving Average (ARIMA), Seasonal ARIMA (SARIMA),
SARIMA with eXogenous factors (SARIMAX), and Holt-Winter for single-step
prediction. The seasonality of our traffic has been explicitly modeled using
SARIMA, which reduces the rolling prediction Mean Average Percentage Error
(MAPE) by more than 4% compared to ARIMA (incapable of handling the
seasonality). We further improved traffic prediction using SARIMAX to learn
different exogenous factors extracted from the original traffic, which yielded
the best rolling prediction results with a MAPE of 6.83%. Finally, we applied
the exponential smoothing technique to handle the variability in traffic
following the Holt-Winter model, which exhibited a better prediction than
ARIMA (around 1.5% less MAPE). The rolling prediction technique reduced
prediction error using real Internet Service Provider (ISP) traffic data by
more than 50% compared to the standard prediction method.

###### Index Terms:

IP traffic prediction, ISP, internet traffic, rolling prediction, statistical
model

##  I Introduction

Internet traffic engineering deals with the technology, principle, technique,
and tool that assist network administrators in evaluating and optimizing the
operational IP network performance. Also the traffic prediction and
forecasting are some of the most crucial parts affecting network performance
[1]. It helps enhance the network Quality of Service (QoS) and Quality of
Experience (QoE). Also, traffic forecasting assists ISP providers in their
business decisions such as new product development and service decommissions,
advertising, pricing, traffic migration, etc., based on traffic forecasting
results. In addition, accurate forecasting helps service providers in capacity
planning and investment optimization. Therefore, selecting an appropriate
methodology for network traffic prediction is critical for ISP business.

Currently, ISP providers highly depend on experienced network administrators,
and they follow an instinctive approach to forecast future traffic using the
market analysis data such as a possible number of customers and their usage
behavior [2]. The factors they considered for their prediction can be divided
into two main categories: internal and external factors. Internal factors are
related to ISP companies, such as introducing new services, traffic migration,
speed up-gradation, etc. In contrast, external factors come from outside, such
as new internet applications, regional economic factors, seasonal effects,
etc. Therefore, the intuitive method can only predict a rough estimation of
future traffic, which is inadequate to make business decisions.

On the other hand, Operational Research, Statistics, and Computer Science
contributions led to reliable prediction methods that replaced intuition-based
ones. In particular, Time Series Forecasting (TSF), also termed as univariate
forecasting, discusses scientific ways to predict chronologically ordered
data, called time-series data [3]. The ultimate objective of TSF is to model a
complex forecasting system, predicting future behavior based on the historical
observation. The TSF models can be categorized into three main categories
based on their learning techniques: statistical model, machine learning, and
deep learning model. The prediction model from these categories requires
historical data to learn the general trend in time series and make inferences
about the future. The learning models also demand different settings of hyper-
parameters based on the dataset size and complexity. The prediction model’s
learning capability and accuracy directly rely on these parameter
configurations.

Traditional statistical forecasting models such as Auto-Regressive (AR), Auto-
Regressive Moving Average (ARMA), Auto-Regressive Integrated Moving Average
(ARIMA), AutoRegressive Integrated Moving Average (ARIMA), Seasonal
Autoregressive Integrated Moving Average (SARIMA), Seasonal Auto-Regressive
Integrated Moving Average with eXogenous factors (SARIMAX), Holt-Winter, etc.,
were studied extensively at different time-series domains. These classical
models have been used in traffic load forecasting [4], cloud traffic
prediction [5], electricity load forecasting [6], and so on. However, the
statistical models are best at capturing the Linear and Short Range
Dependencies (SRD) in time-series data but exhibit the least performance in
handling Long Range Dependencies (LRD), which results in poor time-series
prediction and forecasting [7]. In addition, they are also incapable of
learning the non-linear attribute of the time-series data. As a result, many
variations of the classical forecasting model such as Fractionally Integrated
Autoregressive Moving Average (FARIMA) [8], or hybrid models such as ARIMA-
GARCH [9] have been proposed to improve the statistical forecasting model’s
performance. This research explored several state-of-the-art statistical
prediction methods such as ARIMA, SARIMA, SARIMAX, and Holt-Winter to predict
internet traffic volume. Also, we implemented a technique called rolling
prediction to improve the performance of the standard prediction. In rolling
prediction, the model used the validation data of the recent prediction for
re-training the corresponding model. The comparative analysis among different
prediction models summarizes their overall performance in traffic prediction.
This research shows a significant improvement in traffic prediction using the
rolling technique for our statistical prediction models. The main
contributions of this work are as follows:

  * •

A comparative performance analysis between classical forecasting model in
traffic prediction.

  * •

Extracting new features from time series data for better prediction.

  * •

Achieved significant performance improvement for the traditional model by
applying the rolling prediction technique.

This paper is organized as follows. Section II describes the literature review
of current traffic prediction using statistical models. Section III presents
the proposed methodology introducing the rolling prediction model to the
traditional model. Section IV presents the performance results of the
different prediction methods and draws a comparative picture between standard
prediction and rolling prediction. Finally, section V concludes our paper and
sheds light on future research directions.

##  II Literature Review

Khashei and Bijari [10] proposed an ensemble forecasting model to improve
accuracy. Their model is comprised of a statistical model called Auto-
Regressive Integrated Moving Average (ARIMA) and Artificial Neural Network
(ANN) model. They identified the limitation of the ANN model in handling
linear data, which motivated them to apply a hybrid model based on Multi-Layer
Perceptron (MLP) to process the non-linear part of the time-series data. ARIMA
model handles the linear component in the time series data. A hybrid
forecasting model [10] combining the ARIMA and ANN has been developed to
improve overall forecasting accuracy. Their model has been tested on three
well-known data sets and indicated an improved performance for all of them.
The first stage of their methodology was to generate the required data using
the ARIMA model, and then this data is fed into the ANN model to predict the
future. U. Kumar and V. Jain [11] investigated the performance of the Auto-
Regressive Moving Average (ARMA) and ARIMA model to predict the future in
advance. They fine-tune the model parameters $p$, $q$, and $d$ by
experimenting with different information criteria such as AIC (Akaike
Information Criterion), HIC (Hannon–Quinn Information Criterion), BIC
(Bayesian Information criterion), and FPE (Final Prediction Error). They also
consider AutoCorrelation function (ACF) and Partial AutoCorrelation Function
(PACF) plot information to identify the best performing model. Different model
performance evaluation metric such as MAPE (Mean Absolute Percentage Error),
MAE (Mean Absolute Error) and RMSE (Root Mean Squared Error) has been
considered in their work. M. Dastorani et al. [12] performed a comparative
analysis among different statistical models such AR (Auto-Regressive), MA
(Moving Average), ARMA, ARIMA, and SARIMA. They decomposed the original time
series to process the random component using AR, MA, and ARMA models. A trial
and error approach has been adapted in their research to find out the best-
performing model, and they identified the stochastic model most appropriate
for their problem. H. Liu et al. [13] compared the performance of two
different hybrid models such as ARIMA-ANN and ARIMA-Kalman. These models have
been applied to process the non-stationary data, and both showed better
performance. The ARIMA part of the hybrid model is used to identify the
architecture of the ANN in the case of the ARIMA-ANN model. In contrast, it is
used to initialize the Kalman measurement for the ARIMA-Kalman model. Huda M.
A. El Hag and Sami M. Sharif [14] identified the weakness of the ARIMA model
in long-term prediction and suggested an adjustment in the ARIMA model to
solve the problem. Their proposed Adjusted ARIMA (AARIMA) model can handle an
extra parameter called self-similarity compared to the traditional ARIMA
model. They used four different Hurst estimators to calculate the self-
similarity. Their proposed model shows better hourly internet traffic
prediction in comparison with the ARIMA model.

The modification of conventional statistical models or their combination has
been proposed earlier to improve the prediction in internet traffic based on
the publicly available dataset. In this research, we focused on improving the
performance of internet traffic prediction by adapting a new training strategy
instead of modifying or combining the state-of-the-art statistical models. Our
experiment results show a substantial performance improvement in traffic
prediction after applying the new training strategy

##  III Methodology

In this section, we first introduce the real IP traffic dataset used in our
experiment in subsection III-A. The data requires some preprocessing steps
explained in subsection III-B to clean and make it compatible with the
prediction model. Then, we describe the feature extraction techniques to
extricate the exogenous attributes from our original dataset in subsection
III-C. Next, the rolling prediction method that significantly improves our
model accuracy is introduced in the subsection III-D. After that, we explain
the mathematical background of the prediction model and the performance metric
to evaluate them in subsection III-E and III-F respectively. Finally, we
summarize the configuration of our experimental environment in subsection
III-G.

###  III-A Dataset

Real internet traffic telemetry on several high speed interface have been used
for this experiment. The data are collected every five minutes for a recent
thirty day time period. The data contains average generated traffic per five
minutes in bit per second (bps). There are 8563 data samples in our dataset
consisting of 29 days complete (288 data instances per day) and last day
incomplete data (211 data instance for 30th day).

###  III-B Data Preprocessing

The dataset consists of time series data for the entire thirty days collected
at every 5 minutes interval. The original dataset was collected in the JSON
format, which is incompatible for processing using the prediction model. So,
it was converted into CSV format first before starting the model prediction.
Only the timestamp (GMT) and traffic data (bps) are taken from the JSON file,
and all other information is discarded. The last day data in the dataset is
incomplete by the previous eight hours (approximately) data and is removed
from the original time series data for the experiment. Ultimately, a total of
29 days of data were considered for developing our prediction model. In
addition, there are some missing values in the dataset filled up using the
mean value. Finally, the time series value unit has been changed from bps to
Gbps (gigabit per second) as the original value is large for feeding into our
statistical model.

###  III-C Feature Extraction

Feature engineering is an essential part of any prediction model, whether a
statistical or machine learning model. Some features have been derived from
the time series data for better prediction. We divided a day into a total
seven different portion as mid-night (12 am-3 pm), late-night (3 am-6 am),
early morning (6 am - 9 am), morning (9 am – 12 am), afternoon (12 pm -15 pm),
late afternoon (15 pm - 18 pm), evening (18 pm - 21 pm), night (21 pm - 24
pm). We introduce another feature based on weekdays and weekends since there
is a chance of high traffic on weekdays. These features are particularly
provided in the SARIMAX model as it is capable of handling exogenous
attributes.

###  III-D Rolling Prediction

Our dataset is divided into train and test sets for training and evaluating
the prediction model, respectively. The proposed model used a prediction
function, which takes the testing data indices as a parameter for immediate
prediction. We looped through our test set entries to make inferences for all
test instances. After evaluating each test instance, the model has been re-
trained using the true observation from the test set in rolling prediction. In
contrast, the standard prediction technique estimates all predictions in the
test set after training the model once. This rolling prediction technique
exhibits improved performance for all prediction models.

###  III-E Time Series Forecasting Models

####  III-E1 ARIMA

ARIMA model was proposed by Box and Jenkins and is also known as the Box-
Jenkins methodology. This model predicts the future value based on the past
values of the time series, that is, its own lagged values and the lagged
forecast white noises. The time series need to be stationary before applying
the ARIMA model as it performs well when there is no correlation and
dependency among the predictors [15]. Total of three parameters, such as order
of AR term $(p)$, order of MA term $(q)$, and the number of differencing to
make the time series stationary $(d)$ are required to design the ARIMA
$(p,q,d)$ model. We can express the ARIMA model mathematically as follow [16]:

$\Phi(L)^{p}\ \Delta^{d}y_{t}=\upphi(L)^{q}\ \Delta^{d}\epsilon_{t}$ (1)
$\Delta^{d}y_{t}=y_{t}^{d-1}-y_{t-1}^{d-1}$ (2)

  * •

Here, $y_{t}$ is the time series.

  * •

$p$, $d$, and $q$ are referred to as the order of AR, I and MA components of
the ARIMA model.

  * •

$\Delta^{d}$ is an operator to make the $y_{t}$ stationary.

  * •

$\Phi(L)^{p}$ id the lag polynomials of order $p$, and $L$ is defined as the
lag operator.

  * •

$\epsilon_{t}$ is white noise.

####  III-E2 SARIMA

SARIMA is a generalized version of the ARIMA model which can handle the
seasonality in the time series data. The SARIMA model requires additional four
parameters such as seasonal autoregressive order $(P)$, seasonal moving
average order $(Q)$, the seasonal difference $(D)$, and the length of the
seasonality period $(S)$ to process the seasonal component in the series. As a
result, six parameters are necessary to define the SARIMA model where $p$,
$d$, and $q$ are the same as the ARIMA model. We can express SARIMA using the
following mathematical equation [17].

$\Phi(L)^{p}\Phi(L^{S})^{P}\
\Delta^{d}y_{t}\Delta_{S}^{D}y_{t}=\upphi(L)^{q}\upphi(L^{S})^{Q}\
\Delta^{d}\epsilon_{t}\Delta_{S}^{D}\epsilon_{t}$ (3)
$\Delta^{d}y_{t}=y_{t}^{d-1}-y_{t-1}^{d-1}$ (4)

  * •

Here, $y_{t}$ is a time series with seasonality $S$.

  * •

$P,D,\text{and }Q$ represent the similar meaning of $p,q,\text{and }d$ in the
ARIMA model but rather applicable for seasonal lags.

####  III-E3 SARIMAX

SARIMAX model capable of processing the exogenous features of the time series.
The exogenous attribute calculated at time $t$ impacts the not auto-regressive
time series value at time $t$. We can change the SARIMA equation 3 above to
make it an equivalent SARIMAX equation as follows [18].

$\displaystyle{\Phi(L)^{p}\Phi(L^{S})^{P}\
\Delta^{d}y_{t}\Delta_{S}^{D}y_{t}=\upphi(L)^{q}\upphi(L^{S})^{Q}\
\Delta^{d}\epsilon_{t}\Delta_{S}^{D}\epsilon_{t}+}$ (5) $\displaystyle{\ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\sum_{i=1}^{n}\beta_{i}x_{t}^{i}}$

  * •

Here, $x_{t}^{i}$ is the exogenous attribute at time $t$ and $n$ is the total
number of exogenous features.

  * •

$\beta_{i}$ is the coefficient for the variable $x_{i}$

####  III-E4 Holt-Winter

The Holt-Winters method is also known as Triple Exponential Smoothing, one of
the popular algorithms designed for time-series forecasting. The first studies
of Exponential Smoothing are back to the Simeon Poisson; after him in 1956,
Robert Brown has introduced its forecasting application. This forecasting
model is defined in three equations one for level $(l_{t})$, one for trend
$(b_{t})$, and one for seasonality $(s_{t})$ and also forecasting equation
with smoothing parameters $\alpha$, $\beta$, and $\gamma$, respectively. Holt-
Winters method has two variations, additive and multiplicative, that differ
like the seasonal data values. We used an additive version of the Holt-Winter
model for our experiment, and the equation can be defined as follow [19]:

$l_{t}=\alpha(y_{t}-s_{t-m})+(1-\alpha)(l_{t-1}+b_{t-1})$ (6)
$b_{t}=\beta(l_{t}-l_{t-1})+(1-\beta)b_{t-1}$ (7)
$s_{t}=\gamma(y_{t}-l_{t-1}-b_{t-1})+(1-\gamma)s_{t-m}$ (8) $\hat{y}_{t+h\mid
t}=l_{t}+hb_{t}+s_{t+h-m(k+1)}$ (9)

  * •

Here, $m$ represents seasonal period.

  * •

$\beta_{i}$ is the coefficient for the variable $x_{i}$

  * •

$h$ is forecast horizon

  * •

$k$ is the integer part of $(h-1)/m$.

###  III-F Evaluation Metrics

We used Mean Absolute Percentage Error (MAPE), to estimate the performance of
our traffic forecasting models. The performance metric identify the deviation
of the predicted result from the original data. For example, MAPE error
represents the average percentage of fluctuation between the actual value and
predicted value. Therefore, we can define our performance metric
mathematically as follow:

$MAPE=\dfrac{1}{n}\sum_{i=1}^{n}\bigg{|}\dfrac{p_{i}-o_{i}}{o_{i}}\bigg{|}\times
100\%$ (10)

Here, $p_{i}$ and $o_{i}$ are predicted and original value respectively, and
$n$ is the total number of test instance

###  III-G Software and Hardware Preliminaries

We have used Python and the statistical model library statsmodels [20] to
conduct the experiments. Our computer has the configuration of Intel (R)
i3-8130U<EMAIL_ADDRESS>8GB memory, and a 64-bit Windows operating system.

TABLE I: Top five best performing model parameter SL |  (p, d, q) |  AIC
---|---|---
1 | (13, 1, 16) | 3782.588307
2 | (21, 1, 2) | 3782.751381
3 | (14,1,17) | 3783.706900
4 | (13, 1, 18) | 3785.061129
5 | (21, 1, 3) | 3785.332434
6 | (22, 1, 3) | 3786.325772
7 | (16, 1, 17) | 3786.547150
8 | (17,1,17) | 3786.761965
9 | (22, 1, 4) | 3786.964070
10 | (20, 1, 2) | 3787.326614
Figure 1: Decomposed time series data

(a)

(b)

Figure 2: Partial and Auto Partial Correlation Analysis TABLE II: Performance summary for all model Model | MAPE (Standard Prediction) | MAPE (Rolling Prediction)
---|---|---
ARIMA | 20.09% | 7.71%
SARIMA | 21.03% | 7.34%
SARIMAX | 21.58% | 6.83%
Holt-Winter | 18.29% | 7.59%
TABLE III: Comparison between actual traffic and predicted traffic Prediction Model |  Prediction |  Rolling Prediction
---|---|---
ARIMA |  [-.5] |  [-.5]
SARIMA |  [-.5] |  [-.5]
SARIMAX |  [-.5] |  [-.5]
Holt-Winter |  [-.5] |  [-.5]

##  IV Results and Discussion

We considered 21 days of data for training our model and the last eight days
for testing. Before applying the prediction model, several time-series
characteristics, such as stationarity, seasonality, tread, etc., need to be
identified. For example, the ARIMA model performs better in stationary time-
series data. We can check the stationarity of the time series using the
Augmented Dicky Fuller (ADF) test. According to the test results, the ADF
statistics (-1.791022) is greater than the critical values (- 3.440), and the
$p$-value (0.38) is more significant than 0.05, which satisfies its non-
stationarity. To make the time series stationary, we took the first order log
difference of our data and again performed the ADF test. Since the $p$-value
(0.00) is less than 0.05 and ADF statistics (-19.49) is less than the critical
value (-3.440), so we concluded that the time series is stationary after
taking the first-order difference. This experiment helped us to set the
difference order $d$ to 1 in defining our ARIMA model. Next, we decomposed our
time series to identify the seasonality and trend in Fig. 1. There is a clear
seasonality in our time series, which repeats every 24 hours that is there is
a daily seasonality in our dataset.

After that, we plotted the ACF and PACF in Fig. 2 to figure out the hyper-
parameters $p$ and $q$ for the ARIMA model. It is difficult to define AR($p$)
and MA($q$) order from these plots as there is a sinusoidal pattern and no
clear indication of the lag, which is significant for our time series data.
That is why a grid search technique has been applied to identify the best
hyper-parameter combination of $p$, $q$, and $d$ based on the minimum Akaike
Information Criterion (AIC) value. We tested different $p$ and $q$
combinations from 0 to 24 with difference order 1 to find the best parameter
set for our ARIMA model. The top ten combinations based on minimum AIC value
are shown in Table I. Our grid search result indicates the order (13, 1, 16)
is the best combination with minimum AIC of 3782.588307 for the ARIMA model.
Similarly, the best model parameter we identified for SARIMA model is ($p$,
$d$, $q$) = (13, 1, 16) and ($P$, $D$, $Q$, $S$) = (1, 0, 1, 24). The best
performing parameter for SARIMAX model is ($p$, $d$, $q$) = (13, 1, 16) and
($P$, $D$, $Q$, $S$) = (1, 0, 1, 24). We also provided the exogenous features
extracted from our original traffic data according to the feature extraction
process discussed in subsection III-C into the SARIMAX model. Finally, we
fine-tuned three hyper-parameters parameters for the Holt-Winter model: trend
type, seasonality type, and seasonality period. The additive version of the
Holt-Winter model is used in our experiment since the seasonal variation is
relatively constant over time.

We experimented with two different types of prediction methodology. The
standard prediction method is mainly divided into two stages: train and test.
But rolling prediction model used the validation data of the latest prediction
to re-train the model after each inference. The result presented in Table II
shows improved performance for all models in rolling prediction if we compare
it with standard prediction. The model performance evaluation metrics MAPE in
standard prediction decreased to more than 50% in rolling prediction for every
model. Firstly, we applied ARIMA model with the best parameter combination
($p=13$, $d=1$, and $q=16$), and it resulted in average percentage prediction
error of 7.71 Gbps. Since ARIMA cannot handle the seasonality in the time-
series data, we implemented another model, SARIMA, accepting seasonality
information as a hyper-parameter. SARIMA provides better results with a 7.34%
average deviation from the original traffic in the rolling prediction,
reducing the prediction error by more than 4% compared to ARIMA. Next, we
extracted some features from our original dataset according to the feature
extraction process described in section III-C. The original traffic and their
exogenous attributes are then trained using SARIMAX, which can handle extra
features along with traffic seasonality. The SARIMAX yields the best
prediction result in our experiment with the lowest MAPE of 6.83%, decreasing
the error by more than 11% and 6% compared to ARIMA and SARIMA, respectively.
Finally, we implement Holt-Winter prediction model to handle the variability
in our traffic data by using the exponential smoothing technique. The Holt-
Winter shows 7.59% average fluctuation between predicted and actual traffic
using rolling prediction, and it is 1.5% less error than ARIMA. Our best-
performing model SARIMAX reduces the prediction error by more than 10% in
comparison to the Holt-winter. In Table III, we depicted the actual and
expected traffic for the last eight days using both standard prediction and
rolling prediction. The comparison between actual and predicted traffic shows
a significantly better fitting using the rolling prediction technique.

##  V Conclusion

In this research, we experimented with several internet traffic prediction
models. We tried to improve the performance of state-of-the-art prediction
models to learn the general trend in real IP traffic. Also, a comparative
performance analysis among several conventional statistical models has been
conducted in traffic prediction. Our experimental results show a significant
improvement in traffic prediction when we feed our model with the validation
data after each prediction. We considered four different prediction models,
e.g., ARIMA, SARIMA, SARIMAX, and Holt-Winter, showing a better accuracy using
the rolling prediction technique than the standard prediction. As part of our
future work, we would like to model the residual variation using
Autoregressive Conditional Heteroskedasticity (ARCH) and Generalized
AutoRegressive Conditional Heteroskedasticity (GARCH). Also, we plan to extend
our work from single-step traffic prediction to multi-step prediction.

## References

  * [1] T. Otoshi, Y. Ohsita, M. Murata, Y. Takahashi, K. Ishibashi, and K. Shiomoto, “Traffic prediction for dynamic traffic engineering,” _Computer Networks_ , vol. 85, pp. 36–50, 2015.
  * [2] K. Papagiannaki, N. Taft, Z.-L. Zhang, and C. Diot, “Long-term forecasting of internet backbone traffic,” _IEEE transactions on neural networks_ , vol. 16, no. 5, pp. 1110–1124, 2005.
  * [3] Z. Shen, Y. Zhang, J. Lu, J. Xu, and G. Xiao, “A novel time series forecasting model with deep learning,” _Neurocomputing_ , vol. 396, pp. 302–313, 2020\.
  * [4] S. Nyaramneni, M. A. Saifulla, and S. M. Shareef, “Arima for traffic load prediction in software defined networks,” in _Evolutionary Computing and Mobile Sustainable Networks_. Springer, 2021, pp. 815–824.
  * [5] H. Mehdi, Z. Pooranian, and P. G. Vinueza Naranjo, “Cloud traffic prediction based on fuzzy arima model with low dependence on historical data,” _Transactions on Emerging Telecommunications Technologies_ , vol. 33, no. 3, p. e3731, 2022.
  * [6] B. Nepal, M. Yamaha, A. Yokoe, and T. Yamaji, “Electricity load forecasting using clustering and arima model for energy management in buildings,” _Japan Architectural Review_ , vol. 3, no. 1, pp. 62–76, 2020.
  * [7] C. Katris and S. Daskalaki, “Comparing forecasting approaches for internet traffic,” _Expert Systems with Applications_ , vol. 42, no. 21, pp. 8172–8183, 2015.
  * [8] F. Wu, C. Cattani, W. Song, and E. Zio, “Fractional arima with an improved cuckoo search optimization for the efficient short-term power load forecasting,” _Alexandria Engineering Journal_ , vol. 59, no. 5, pp. 3111–3118, 2020.
  * [9] Q. T. Tran, Z. Ma, H. Li, L. Hao, and Q. K. Trinh, “A multiplicative seasonal arima/garch model in evn traffic prediction,” _International Journal of Communications, Network and System Sciences_ , vol. 8, no. 4, p. 43, 2015.
  * [10] M. Khashei and M. Bijari, “A novel hybridization of artificial neural networks and arima models for time series forecasting,” _Applied Soft Computing_ , vol. 11, no. 2, pp. 2664–2675, 2011.
  * [11] U. Kumar and V. Jain, “Arima forecasting of ambient air pollutants (o 3, no, no 2 and co),” _Stochastic Environmental Research and Risk Assessment_ , vol. 24, no. 5, pp. 751–760, 2010.
  * [12] M. Dastorani, M. Mirzavand, M. T. Dastorani, and S. J. Sadatinejad, “Comparative study among different time series models applied to monthly rainfall forecasting in semi-arid climate condition,” _Natural Hazards_ , vol. 81, no. 3, pp. 1811–1827, 2016.
  * [13] H. Liu, H.-q. Tian, and Y.-f. Li, “Comparison of two new arima-ann and arima-kalman hybrid methods for wind speed prediction,” _Applied Energy_ , vol. 98, pp. 415–424, 2012.
  * [14] H. M. El Hag and S. M. Sharif, “An adjusted arima model for internet traffic,” in _AFRICON 2007_. IEEE, 2007, pp. 1–6.
  * [15] R. Madan and P. S. Mangipudi, “Predicting computer network traffic: a time series forecasting approach using dwt, arima and rnn,” in _2018 Eleventh International Conference on Contemporary Computing (IC3)_. IEEE, 2018, pp. 1–5.
  * [16] G. E. Box, G. M. Jenkins, G. C. Reinsel, and G. M. Ljung, _Time series analysis: forecasting and control_. John Wiley & Sons, 2015.
  * [17] T. Otoshi, Y. Ohsita, M. Murata, Y. Takahashi, K. Ishibashi, and K. Shiomoto, “Traffic prediction for dynamic traffic engineering considering traffic variation,” in _2013 IEEE Global Communications Conference (GLOBECOM)_. IEEE, 2013, pp. 1570–1576.
  * [18] F. Sheng and L. Jia, “Short-term load forecasting based on sarimax-lstm,” in _2020 5th International Conference on Power and Renewable Energy (ICPRE)_ , 2020, pp. 90–94.
  * [19] N. Talkhi, N. A. Fatemi, Z. Ataei, and M. J. Nooghabi, “Modeling and forecasting number of confirmed and death caused covid-19 in iran: A comparison of time series forecasting methods,” _Biomedical Signal Processing and Control_ , vol. 66, p. 102494, 2021.
  * [20] S. Seabold and J. Perktold, “statsmodels: Econometric and statistical modeling with python,” in _9th Python in Science Conference_ , 2010.