# Strengthening a theorem of Meyniel

Quentin Deschamps Carl Feghali František Kardoš
Clément Legrand-Duchesne Théo Pierron∗ Univ. Lyon, Université Lyon 1, LIRIS
UMR CNRS 5205, F-69621, Lyon, France.Univ. Lyon, EnsL, UCBL, CNRS, LIP,
F-69342, Lyon Cedex 07, France. Faculty of Mathematics, Physics, and
Informatics, Comenius University, Bratislava, Slovakia.LaBRI, CNRS, Université
de Bordeaux, Bordeaux, France.

###### Abstract

For an integer $k\geq 1$ and a graph $G$, let $\mathcal{K}_{k}(G)$ be the
graph that has vertex set all proper $k$-colorings of $G$, and an edge between
two vertices $\alpha$ and $\beta$ whenever the coloring $\beta$ can be
obtained from $\alpha$ by a single Kempe change. A theorem of Meyniel from
1978 states that $\mathcal{K}_{5}(G)$ is connected with diameter
$O(5^{|V(G)|})$ for every planar graph $G$. We significantly strengthen this
result, by showing that there is a positive constant $c$ such that
$\mathcal{K}_{5}(G)$ has diameter $O(|V(G)|^{c})$ for every planar graph $G$.

##  1 Introduction

Let $k$ be a positive integer, and let $G$ be a graph. A _Kempe chain_ in
colors $\\{a,b\\}$ is a maximal connected subgraph $B$ of $G$ such that every
vertex of $B$ has color $a$ or $b$. By swapping the colors $a$ and $b$ on $B$,
a new coloring is obtained. This operation is called a _K-change_. Let
$\mathcal{K}_{k}(G)$ be the graph that has vertex set all proper $k$-colorings
of $G$, and an edge between two vertices $\alpha$ and $\beta$ whenever the
coloring $\beta$ can be obtained from $\alpha$ by a single K-change.

A graph $G$ is $d$-degenerate if every subgraph of $G$ contains a vertex of
degree at most $d$. Las Vergnas and Meyniel [1] proved the following result.

######  Theorem 1.1.

If $G$ is $d$-degenerate and $k>d$ is an integer, then $\mathcal{K}_{k}(G)$ is
connected.

Meyniel [2] strengthened this result for planar graphs by proving

######  Theorem 1.2.

If $G$ is a planar graph, then $\mathcal{K}_{5}(G)$ is connected.

Here, the number $5$ of colors cannot be replaced by $4$ [3].

The proof of Theorem 1.1 implies that $\mathcal{K}_{k}(G)$ has diameter
$O(d^{|V(G)|})$. Similarly, the proof of Theorem 1.2 implies that
$\mathcal{K}_{5}(G)$ has diameter $O(5^{|V(G)|})$. Bonamy, Bousquet, Feghali
and Johnson [4] conjectured that the former can be significantly improved.

######  Conjecture 1.3.

If $G=(V,E)$ is $d$-degenerate and $k>d$ is an integer, then
$\mathcal{K}_{k}(G)$ has diameter $O(|V|^{2})$.

For $k>d+1$ in Conjecture 1.3, a result of Bousquet and Heinrich [5] regarding
the reconfiguration graph for colorings of graphs with bounded degeneracy
implies the following

######  Theorem 1.4.

If $G=(V,E)$ is $d$-degenerate and $k>d+1$ is an integer, then
$\mathcal{K}_{k}(G)$ has diameter $O(|V|^{d+1})$.

The case $k=d+1$ in Conjecture 1.3 seems much more challenging. It was only
until very recently that Bonamy, Delecroix and Legrand-Duschene [6] addressed
this case for some proper subclasses of degenerate graphs such as graphs with
bounded treewidth, graphs with bounded maximum average degree and
$(\Delta-1)$-degenerate graphs, where $\Delta$ denotes the maximum degree.

The object of this paper is to break the $k>d+1$ (in fact, $k=d+1$) barrier
for the class of planar graphs, by proving the following strengthening of
Theorem 1.2.

######  Theorem 1.5.

If $G$ is a planar graph, then $\mathcal{K}_{5}(G)$ has diameter at most a
polynomial in the number of vertices of $G$.

The proof of Theorem 1.5 is based on a proof method introduced in [7] and some
ideas from [8], but of course has many of its own features. We also note that
the proof can be adapted for a larger number of colors.

The paper is organized as follows. In Section 2, we analyze a proof from [9]
to get an essential estimate (Proposition 2.2) that we use in Section 3 to
prove Theorem 1.5.

##  2 The case of 3-colorable planar graphs

In [9], Mohar has proved the following theorem.

######  Theorem 2.1.

Let $G$ be a $3$-colorable planar graph. Then $\mathcal{K}_{4}(G)$ is
connected.

In this section we show, by making a few simple observations, that the proof
of Theorem 2.1 in fact establishes the following stronger fact that is crucial
to our main theorem.

######  Proposition 2.2.

Let $G$ be a $3$-colorable planar graph on $n$ vertices. Then for every
$4$-colorings $\alpha$ and $\beta$ of $G$, there exists a sequence of
K-changes from $\alpha$ to $\beta$ that changes the color of each vertex at
most $O(n^{2})$ times.

In particular, note that Proposition 2.2 implies that $\mathcal{K}_{4}(G)$ has
diameter $O(n^{3})$ when $G$ is a 3-colorable planar graph on $n$ vertices.
The rest of this section is devoted to the proof of Proposition 2.2. We
basically follow the same steps as the proof of Theorem 2.1, except that we
add some complexity estimates to:

  * •

a result of Fisk [10] that handles the case of 3-colorable planar
triangulations, and

  * •

a reduction to the (preceding) case as done by Mohar [9].

### Analysing the result of Fisk.

We consider the theorem of Fisk [10] stated below.

######  Theorem 2.3.

Let $G$ be a $3$-colorable triangulation of the plane. Then
$\mathcal{K}_{4}(G)$ is connected.

We show that the proof of Theorem 2.3, word for word, gives us the following
estimate.

######  Lemma 2.4.

Let $G$ be a $3$-colorable triangulation of the plane on $n$ vertices. Then
for every $4$-colorings $\alpha$ and $\beta$ of $G$, there exists a sequence
of K-changes from $\alpha$ to $\beta$ that changes the color of each vertex at
most $O(n^{2})$ times.

Let $f$ be a $4$-coloring of a triangulation $G$ of the plane, and let $e=xy$
be an edge of $G$. We denote by $f(e)=\\{f(x),f(y)\\}$ the color of $e$ under
$f$. If $xyz$ and $xyw$ are the two triangles containing an edge $xy$, we say
that $xy$ is _singular under $f$_ if $f(w)=f(z)$.

To prove Lemma 2.4, we require the following key structural lemma extracted
from the proof of Theorem 2.3. Note that, for a triangulation $G$ of the
plane, if $G$ has a $3$-coloring, then this coloring is unique up to
permutations of colors, and all edges are singular.

######  Lemma 2.5.

Let $G$ be a $3$-colorable triangulation of the plane, and let $f$ be a
$4$-coloring of $G$. Then every monochromatic set of non-singular edges of $G$
contains a cycle that bounds some region of the plane.

Using this lemma, we can conclude the proof of Lemma 2.4.

###### Proof of Lemma 2.4.

We shall show, by exhibiting at most $n\cdot|E(G)|$ K-changes, that any
$4$-coloring of $G$ is K-equivalent to the (unique) $3$-coloring of $G$, which
will prove the lemma. To do so, given a $4$-coloring $f$ of $G$, we show how
to obtain, via at most $n$ K-changes, a new coloring $g$ that is
$K$-equivalent to $f$ and with fewer non-singular edges. As no edge in the
$3$-coloring of $G$ is non-singular, by iterating this argument at most
$|E(G)|$ times the result will follow.

Let $e$ be a non-singular edge of $G$. By Lemma 2.5, there is a cycle in $G$
whose edges have the same color as $e$ and bounding some region $D$ of the
plane. By interchanging the two colors in $\\{1,2,3,4\\}\setminus f(e)$ in the
interior of $D$, we obtain a new coloring $g$ with fewer non-singular edges
than $f$ (singular edges in the interior of $D$ stay singular, while the edges
of the cycle on the boundary of $D$ change from non-singular to singular). ∎

### Reduction to the triangulation case.

We first restate Proposition 4.3 from Mohar [9] except that we add some
observations about the number of vertices of the resulting triangulation and
the number of K-changes involved – these directly follow from Mohar’s proof.

######  Proposition 2.6.

Let $G$ be a planar graph with a facial cycle $C$ and two 4-colorings $c_{1}$,
$c_{2}$. Then there exists a graph $H$ formed from $G$ by adding a near-
triangulation of size $O(|C|)$ inside $C$ and two 4-colorings
$c^{\prime}_{1},c^{\prime}_{2}$ of $H$ such that $c^{\prime}_{1}|_{V(G)}$ and
$c^{\prime}_{2}|_{V(G)}$ are obtained from $c_{1},c_{2}$ using at most $O(1)$
K-changes. Moreover, if the restriction of $c_{1}$ to $C$ is a $3$-coloring,
then $c^{\prime}_{1}$ is a 3-coloring of $H$ that coincides with $c_{1}$ on
$V(G)$.

We may now prove Proposition 2.2.

###### Proof of Proposition 2.2.

The proof follows the same steps as Theorem 4.4 in [9]. We apply Proposition
2.6 to each face of $G$ (instead of Proposition 4.3 from [9]). We thus made
$O(n)$ K-changes, and the resulting triangulation $T$ has $O(n)$ vertices. We
then apply Lemma 2.4 (instead of Theorem 4.1 from [9]) to obtain a sequence of
K-changes between the two colorings of $T$ that changes the color of each
vertex at most $O(n^{2})$ times. ∎

##  3 Main Theorem

In this section we prove Theorem 1.5. Thanks to the celebrated Four Colour
Theorem, it suffices to prove the following result.

######  Theorem 3.1.

Let $G$ be a plane graph with $n$ vertices. For every $5$-coloring $\alpha$ of
$G$ and every $4$-coloring $\beta$ of $G$, there is a sequence of K-changes
from $\alpha$ to $\beta$ where each vertex is recolored polynomially many
times.

In the remainder of this section, we prove Theorem 3.1. Let us briefly sketch
the details of the approach. The proof proceeds by induction on the number of
vertices. Our aim is to describe a sequence of K-changes from $\alpha$ to
$\beta$ such that each vertex is recolored at most $f(n)$ times, where $f$
will satisfy a recurrence relation given at the end of the section. To
establish this, we roughly adopt the following strategy:

  1. 1.

We find a ‘large’ independent set $I$ that is monochromatic in both $\alpha$
and $\beta$ and that contains vertices of degree at most 6 in $G$ (that $I$ is
‘large’ will enable us to show that $f$ is a polynomial function).

  2. 2.

We introduce an operation at a vertex that we call _collapsing_ , which when
applied to each vertex of $I$ gives a new graph $H$ where the degree of each
vertex of $I$ is at most $4$ in $H$ and such that $F=H-I$ is planar. We use
these to show that any sequence of K-changes in $F$ extends to a sequence of
K-changes in $H$ and, in turn, in $G$.

  3. 3.

We apply induction to find a sequence of K-changes in $F$ from any
$5$-coloring of $F$ to a $4$-coloring of $F$ avoiding the color $\alpha(I)$.
Applying Step 2, this sequence extends to a sequence in $G$ ending at a
$5$-coloring, where color $5$ may appear only on $I$.

  4. 4.

By definition, $I\subset B$ for some color class $B$ of $\beta$; so we can
recolor each vertex of $B$ to color $5$. Finally, noting that $G-B$ is a
3-colorable planar graph, we then apply Proposition 2.2 to recolor the
remaining vertices in $G-B$ to their color in $\beta$.

In the rest of this section, we give the details and conclude with a small
analysis of the maximum number of times a vertex changes its color.

### Step 1: Constructing $I$.

We prove that the required independent set $I$ exists.

######  Lemma 3.2.

There exists an independent set $I$ of $G$ such that:

  * •

all the vertices of $I$ have degree at most 6,

  * •

$I$ is contained in a color class of $\alpha$ and of $\beta$,

  * •

$|I|\geqslant\frac{n}{140}$.

###### Proof.

Let $S$ be the set of vertices of degree at most $6$ in $G$. Then $|S|>n/7$
since otherwise

$\sum_{v\in V(G)}d(v)\geq\sum_{v\in V(G)-S}d(v)\geq 7(n-\frac{n}{7})=6n,$

which contradicts Euler’s formula.

For $i\in\\{1,\dots,5\\}$ and $j\in\\{1,\dots,4\\}$ define the set

$\displaystyle S_{i,j}=S\cap\alpha^{-1}(i)\cap\gamma^{-1}(j).$

Note that each $S_{i,j}$ satisfies all the criteria from the lemma, except
maybe the last. However, by the pigeonhole principle, there exists $i$ and $j$
such that $S_{i,j}$ contains at least $|S|/(5\times 4)\geq n/140$ vertices,
which concludes the proof. ∎

From now on, we fix a set $I$ satisfying the hypotheses of Lemma 3.2.

### Step 2: Constructing $H$ and extending recoloring sequences.

In order to construct $H$, we want to identify vertices in $N(I)$ that are
colored alike so that vertices of $I$ end up with degree 4 and the resulting
graph with $I$ excluded is planar. We show that we can modify the coloring
$\alpha$ so that such identification becomes possible.

Let $N$ be a plane graph. For a $5$-coloring $\varphi$ of $N$ and a vertex $v$
of $K$ with $d(v)=6$, we say that $v$ is _$\varphi$ -good_ if, in $\varphi$,
the vertex $v$ has three neighbors colored alike or two pairs $(a,b)$ and
$(c,d)$ of neighbors colored alike that are _non-overlapping_ , i.e. such that
$N-v+ab+cd$ is planar. A sequence of K-changes is said to _avoid_ color $a$ if
no vertex involved in some K-change in the sequence changes its color to $a$.

######  Lemma 3.3.

Let $N$ be a plane graph, $\alpha$ a 5-coloring of $K$, and $v\in V(N)$ such
that $d(v)=6$. There exists a sequence of at most three K-changes avoiding
$\alpha(v)$ that transforms $\alpha$ to a $5$-coloring $\beta$ of $N$ such
that $v$ is $\beta$-good.

###### Proof.

We can assume that $v$ is not $\alpha$-good. Then the neighbors of $v$ can be
colored in four possible ways (up to permutation of colors). We present a
visual proof in Figure 1, where the six neighbors of $v$ are represented by
circles from left to right in the cyclic ordering around $v$, the numbers
represent their color, a bold circle represents an attempt to perform a Kempe
change, a curved edge between two vertices $u$ and $w$ represents a Kempe
chain containing both $u$ and $w$, and dashed arrows represent the actual
Kempe changes. ∎

Case 1121324121323121324121314Case 2123124123121Solved by case
1423121123421123124123134Solved by case 1Case 4121234121231421231Solved by
case 1124231121234121224121324Solved by case 1Case 3123123Solved by case
2123124123123123113123213 Figure 1: The proof of Lemma 3.3

We may now successively process each vertex of $I$: for each $v\in I$, we
apply Lemma 3.3 to make $v$ $\alpha$-good, then we identify vertices in $N(v)$
so that $v$ has degree at most $4$ in the resulting graph $N^{\prime}$, and
moreover, $N^{\prime}-v$ is planar. Note that since we avoid the color of $v$
in a K-change by Lemma 3.3, the vertices of $I$ are never recolored at any
stage of the construction. We now formalize this intuition.

Let $N$ be a plane graph. Let $\varphi$ be a $5$-coloring of $N$. For a vertex
$v$ of $N$ of degree at most $6$ such that $v$ is $\varphi$-good if $d(v)=6$,
we say that $(N^{\prime},\varphi^{\prime})$ is the result of _collapsing_
$(N,v,\varphi)$ if

  * •

in the case $1\leq d(v)\leq 4$, $N^{\prime}=N$ and $\varphi^{\prime}=\varphi$;

  * •

in the case $d(v)=5$, $N^{\prime}$ is the graph obtained from $N$ by
identifying two neighbors $u$ and $w$ of $v$ with $\varphi(u)=\varphi(w)$ into
a new vertex $z$ and $\varphi^{\prime}$ is the coloring obtained from
$\varphi$ by setting $\varphi^{\prime}(z)=\varphi(u)=\varphi(w))$;

  * •

in the case $d(v)=6$, $N^{\prime}$ is the graph obtained from $N$ by either

    * –

identifying three neighbors $u$, $w$, $z$ of $v$ with
$\varphi(u)=\varphi(w)=\varphi(z)$ into a new vertex $x$. In this case,
$\varphi^{\prime}$ is obtained from $\varphi$ by setting
$\varphi^{\prime}(x)=\varphi(u)$, or

    * –

for two non-overlapping pairs $(u,w)$ and $(x,z)$ of neighbors of $v$ with
$\varphi(u)=\varphi(w)$ and $\varphi(x)=\varphi(z)$, identifying $u$, $w$ into
a new vertex $s$ and $x$, $z$ into a new vertex $t$, and defining
$\varphi^{\prime}$ by setting $\varphi^{\prime}(s)=\varphi(u)$ and
$\varphi^{\prime}(t)=\varphi(x)$.

Since $v$ is $\varphi$-good, it should be immediate that $N^{\prime}$ and
$\varphi^{\prime}$ are well-defined, and that $N^{\prime}-\\{v\\}$ is planar.
We now focus on showing that we can extend a given recoloring sequence of the
collapsed graph to a recoloring sequence of the original graph.

######  Lemma 3.4.

Let $N$ be a plane graph with a $5$-coloring $\varphi$. Let $v$ be a vertex of
$N$ of degree at most $6$ such that $v$ is $\varphi$-good if $d(v)=6$. Let
$(N^{\prime},\varphi^{\prime})$ be the result of collapsing $(N,v,\varphi)$.
Then every sequence of K-changes in $N^{\prime}-\\{v\\}$ starting from
$\varphi^{\prime}\restriction(V(N^{\prime})\setminus\\{v\\})$ extends to a
sequence of K-changes in $N$ starting from $\varphi$. Moreover, each vertex of
$N-v$ changes its color as many times as in $N^{\prime}-v$, and $v$ changes
its color at most once every time one of its neighbors in $N^{\prime}$ changes
its color.

###### Proof.

Each time a neighbor $w$ of $v$ is recolored in $N^{\prime}-v$, we may use the
same K-change in $N^{\prime}$ unless it involves the color of $v$ and there is
another neighbor $u$ of $v$ of the same color as $w$. In this case, since at
most 3 colors appear on $N_{N^{\prime}}(v)$ (as $v$ has degree 4 in
$N^{\prime}$) we can precede this K-change by first recoloring $v$ to a color
not appearing in its neighborhood. This shows that any sequence of K-changes
in $N^{\prime}-\\{v\\}$ extends to a sequence in $N^{\prime}$. To extend the
sequence to $N$, observe that we can simulate in $N$ a K-change in
$N^{\prime}$ at a vertex $w$ by performing a K-change at each vertex that was
identified to form $w$. Clearly, each vertex of $N-v$ changes its color as
many times as in $N^{\prime}-v$, and $v$ changes its color at most once every
time one of its neighbors in $N^{\prime}$ changes its color, which concludes.
∎

We now successively apply Lemma 3.4 to each vertex of $I$ in $G$, so that
every vertex of $I$ becomes $\alpha$-good.

######  Lemma 3.5.

Let $N$ be a plane graph with a $5$-coloring $\varphi$. If $I$ is an
$\varphi$-monochromatic independent set of vertices of degree at most $6$ in
$N$, then there is a $5$-coloring $\psi$ of $N$ for which $I$ is $\psi$-good
and a sequence of K-changes from $\varphi$ to $\psi$ that changes the color of
each vertex at most $3|I|$ times.

###### Proof.

We shall prove by induction on $|V(N)|$ the stronger claim that there is such
a sequence that avoids the color of $I$.

By Lemma 3.3, we can assume, up to at most three K-changes, that $I$ contains
a vertex $v$ of degree at most $5$ or a vertex of degree $6$ that is
$\varphi$-good. So we can let $(N^{\prime},\varphi^{\prime})$ be the result of
collapsing $(N,v,\varphi)$. Since $N^{\prime\prime}=N^{\prime}-\\{v\\}$ is
planar, we can apply our induction hypothesis to $N^{\prime\prime}$ (with
$I-\\{v\\}$ instead of $I$ and $\varphi^{\prime}\restriction N^{\prime\prime}$
instead of $\varphi$) to find a sequence of K-changes in $N^{\prime\prime}$
that avoids the color of $v$ and that transforms $\varphi^{\prime}\restriction
N^{\prime\prime}$ to some $5$-coloring $\varphi^{\prime\prime}$ of
$N^{\prime\prime}$ so that the following holds:

  * •

$I-\\{v\\}$ is $\varphi^{\prime\prime}$-good, and

  * •

each vertex changes its color at most $3(|I|-1)$ times.

By Lemma 3.4, this sequence extends to a sequence in $G$. Moreover, $v$ does
not change its color (as the sequences avoids $\alpha(v)$), and every other
vertex changes its color at most $3(|I|-1)+3=3|I|$ times. This completes the
proof. ∎

### Step 3: Induction.

By Lemma 3.5, we can assume that $I$ is $\alpha$-good. Let $n=|V(G)|$ and
$f(n)$ be the maximum number of times a vertex in $G$ is involved in a
K-change. Write $I=\\{v_{1},\dots,v_{m}\\}$, set $G_{1}=G$, $\psi_{1}=\alpha$
and, for $i=2,\dots,m+1$, let $H_{i-1}=G_{i}-\\{v_{1},\dots,v_{i-1}\\}$ where
$(G_{i},\psi_{i})$ is the result of collapsing $(G_{i-1},v_{i-1},\psi_{i-1})$.
Then the final graph $H_{m}$ is planar. Hence, by the induction hypothesis
combined with Lemma 3.2 (with $H_{m}$ in place of $G$), there is a sequence of
K-changes from $\psi_{m+1}\restriction V(H_{m})$ to some $4$-coloring
$\gamma^{\prime}$ of $H_{m}$ on colors $\\{1,\dots,5\\}\setminus\alpha(I)$
where each vertex changes its color at most $f(n-|I|)$ times. By successively
applying Lemma 3.4 to $H_{m}$ etc. up until $H_{1}$ this same sequence extends
to a sequence in $G$ from $\alpha$ to some $5$-coloring $\psi$ of $G$ where
$\psi\restriction(G-I)$ uses only colors $\\{1,\dots,4\\}$. Now, recalling
Step 4 verbatim, by definition, $I\subset B$ for some color class $B$ of
$\beta$; so we can recolor each vertex of $B$ to color $5$. Finally, noting
that $G-B$ is a 3-colorable planar graph, we then apply Proposition 2.2 to
recolor the remaining vertices in $G-B$ to their color in $\beta$.

### Complexity analysis.

By Lemma 3.5, the color of each vertex is changed at most $3|I|$ times in
order to reach a 5-coloring $\psi$ in which $I$ is $\psi$-good. During the
induction step, by Lemma 3.4, the color of each vertex in $I$ is changed at
most $4f(n-|I|)$ times, while the color of the other vertices changes at most
$f(n-|I|)$ times. By Proposition 2.2, the final step requires at most
$O(n^{2})$ changes of color per vertex. We deduce that $f(n)$ satisfies the
recurrence

$f(n)\leqslant 3|I|+4f(n-|I|)+O(n^{2})\leq
4f\left(n-\frac{n}{140}\right)+O(n^{2}).$

The master theorem then yields that each vertex changes its color
$O(n^{\log_{\frac{140}{139}}(4)})=O(n^{194})$ times, hence the sequence has
length at most $O(n^{195})$, which concludes the proof.

##  4 Acknowledgements

This work was supported by Agence Nationale de la Recherche (France) under
research grants ANR DIGRAPHS ANR-19-CE48-0013-01 and ANR prpdfoject GrR
(ANR-18-CE40-0032) and by APVV-19-0308 and VEGA 1/0743/21.

## References

  * [1] Michel Las Vergnas and Henri Meyniel  “Kempe classes and the Hadwiger Conjecture”  In _Journal of Combinatorial Theory, Series B_ 31.1, 1981, pp. 95–104
  * [2] Henry Meyniel  “Les 5-colorations d’un graphe planaire forment une classe de commutation unique”  In _Journal of Combinatorial Theory, Series B_ 24.3, 1978, pp. 251–257
  * [3] Bojan Mohar  “Akempic triangulations with 4 odd vertices”  In _Discrete mathematics_ 54.1 Elsevier, 1985, pp. 23–29
  * [4] Marthe Bonamy, Nicolas Bousquet, Carl Feghali and Matthew Johnson  “On a conjecture of Mohar concerning Kempe equivalence of regular graphs”  In _Journal of Combinatorial Theory, Series B_ 135 Elsevier, 2019, pp. 179–199
  * [5] Nicolas Bousquet and Marc Heinrich  “A polynomial version of Cereceda’s conjecture”  In _arXiv preprint arXiv:1903.05619_ , 2019
  * [6] Marthe Bonamy, Vincent Delecroix and Clément Legrand-Duchesne  “Kempe changes in degenerate graphs”  In _arXiv preprint arXiv:2112.02313_ , 2021
  * [7] Carl Feghali  “Paths between colourings of sparse graphs”  In _European Journal of Combinatorics_ 75 Elsevier, 2019, pp. 169–171
  * [8] Eduard Eiben and Carl Feghali  “Toward Cereceda’s conjecture for planar graphs”  In _Journal of Graph Theory_ 94.2 Wiley Online Library, 2020, pp. 267–277
  * [9] Bojan Mohar  “Kempe Equivalence of Colorings”  In _Graph Theory in Paris: Proceedings of a Conference in Memory of Claude Berge_ Basel: Birkhäuser Basel, 2007, pp. 287–297
  * [10] Steve Fisk  “Geometric coloring theory”  In _Advances in Mathematics_ 24.3, 1977, pp. 298–340


# Self-Distilled Quantization: Achieving High Compression Rates in
Transformer-Based Language Models

James O’ Neill  Sourav Dutta
Huawei Ireland Research Center
Georges Court, Townsend St, Dublin 2, Ireland
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

We investigate the effects of post-training quantization and quantization-
aware training on the generalization of Transformer language models. We
present a new method called self-distilled quantization (SDQ) that minimizes
accumulative quantization errors and outperforms baselines. We apply SDQ to
multilingual models XLM-R${}_{\text{Base}}$ and InfoXLM${}_{\text{Base}}$ and
demonstrate that both models can be reduced from 32-bit floating point weights
to 8-bit integer weights while maintaining a high level of performance on the
XGLUE benchmark. Our results also highlight the challenges of quantizing
multilingual models, which must generalize to languages they were not fine-
tuned on.

##  1 Introduction

A main aim of neural network quantization is to reduce the size and
computational demands of a model while maintaining its performance. There are
two main approaches: quantization-aware training (QAT) Banner et al. (2018);
Chin et al. (2020); Faghri et al. (2020); Kim et al. (2020); Wang et al.
(2018) and post-training quantization (PTQ) Neill (2020); Bondarenko et al.
(2021); Kim et al. (2021); Dettmers et al. (2022). Both of these approaches
have limitations in terms of dealing with accumulative quantization errors
that are propogated within the layers of a neural network during the forward
pass Zhao et al. (2019); Fan et al. (2020). To address this issue, we propose
a method called Self-Distilled Quantization (SDQ) that combines self-attention
and output distillation with quantization to compress large language models.
SDQ involves injecting quantization noise into the student network during
training and distilling knowledge from a fine-tuned teacher network from both
its final output and outputs of intermediate self-attention layers. By
distilling knowledge of the self-attention layers, as depicted in Figure 1, we
further reduce the compounding effect of quantization errors in the network.
We use SDQ for self-attention models and demonstrate its effectiveness in
compressing multilingual models XLM-R${}_{\text{Base}}$ and
InfoXLM${}_{\text{Base}}$, achieving high compression rates while maintaining
performance on the XGLUE benchmark. Lastly, we identify that quantization
error is largest at the output of self-attention modules.

Figure 1: Self-Attention Self-Distilled Quantization

##  2 Related Work

Combining quantization and distillation has been previously explored by Mishra
and Marr (2017), who used three different schemes to combine low bit precision
and knowledge distillation (KD) using a 4-bit ResNet network. Polino et al.
(2018) used a distillation loss with respect to a quantized teacher network to
train a student network, and also proposed differentiable quantization, which
optimizes the location of quantization points through SGD. Zhou et al. (2017)
used iterative quantization, supervised by a teacher network, to retrain an
FP-32 model with low precision convolution weights (binary, ternary, and 4
bits). Kim et al. (2019) used QAT and fine-tuning to mitigate the
regularization effect of KD on quantized models. Q8BERT Zafrir et al. (2019)
and fully Quantized Transformer Prato et al. (2019) applied QAT with the
Straight-Through Estimator to approximate non-differentiable quantization in
INT-8 format. We now describe the methodology of SDQ.

##  3 Methodology

We begin by defining a dataset $\mathcal{D}:=\\{(X_{i},y_{i})\\}^{D}_{i=1}$
with samples $s_{i}=(X_{i},\bm{y}_{i})$, where each
$X_{i}:=(\bm{x}_{1},\ldots,\bm{x}_{N})$ and $\bm{x}_{i}\in\mathbb{R}^{d}$ is
the $i$-th vector. For structured prediction $y_{i}\in\\{0,1\\}^{N\times
d_{y}}$ and for single and pairwise sentence classification,
$y_{i}\in\\{0,1\\}^{d_{y}}$, where $d_{y}$ is the number of classes. Let
$\bm{y}^{S}=f_{\theta}(X_{i})$ be the output prediction
($y^{S}\in\mathbb{R}^{d_{y}}$) from the student $f_{\theta}(\cdot)$ with
pretrained parameters $\theta:=\\{\text{\bf W}_{l},\bm{b}_{l}\\}_{l=1}^{L}$
for $L$ layers and the outputs of self-attention blocks are denoted as
$\bm{A}_{l}$. The loss function for standard classification fine-tuning is
defined as the cross-entropy loss $\ell_{CE}(\bm{y}^{S},\bm{y})$.

##### Self-Distilled Quantization

For self-distilled quantization, we also require a fine-tuned teacher network
$f_{\Theta}$, that has been tuned from the pretrained state $f_{\theta}$, to
retrieve the soft teacher labels $y^{T}:=f_{\Theta}(\bm{x})$, where
$y^{T}\in\mathbb{R}^{C}$ and $\sum_{c}^{C}y^{T}_{c}=1$. The soft label
$\bm{y}^{T}$ can be more informative than the one-hot targets $\bm{y}$ used
for standard classification as they implicitly approximate pairwise class
similarities through logit probabilities. The Kullbeck-Leibler divergence
(KLD) $\ell_{\mathrm{KLD}}$ is then used with the main task cross-entropy loss
$\ell_{CE}$ to express $\ell_{\mathrm{SDQ}_{\text{KLD}}}$ as shown in Equation
2,

$\begin{split}\ell_{\mathrm{SDQ}_{\text{KLD}}}\text{=}\,\ell_{\mathrm{CE}}(\bm{y}^{S},\bm{y})\text{+}\alpha\tau^{2}D_{\mathrm{KLD}}\big{(}\bm{y}^{S},\bm{y}^{T}\big{)}\end{split}$
(1)

where
$D_{\mathrm{KLD}}(\bm{y}^{S},\bm{y}^{T})=\mathbb{H}(\bm{y}^{T})-\bm{y}^{T}\log(\bm{y}^{S})$,
$\mathbb{H}(\bm{y}^{T})=\bm{y}^{T}\log(\bm{y}^{T})$ is the entropy of the
teacher distribution and $\tau$ is the softmax temperature. Following Hinton
et al. (2015), the weighted sum of the cross-entropy loss and the KLD loss
$\ell_{\mathrm{SDQ}_{\text{KLD}}}\text{=}\ell_{\mathrm{CE}}(\bm{y}^{S},\bm{y})\text{+}\alpha\tau^{2}D_{\mathrm{KLD}}\big{(}\bm{y}^{S},\bm{y}^{T}\big{)}$
is used as our main SDQ-based KD loss baseline, where $\alpha\in[0,1]$.
However, $D_{\mathrm{KLD}}$ only distils the knowledge from the soft targets
of the teacher but does not directly reduce accumulative quantization errors
of the outputs of successive self-attention layers. This brings us to our
proposed attention-based SDQ loss $\ell_{\mathrm{SDQ}_{\text{Att-KLD}}}$ shown
in Equation 2,

$\begin{split}\ell_{\mathrm{SDQ}_{\text{Att-
KLD}}}\text{=}\ell_{\mathrm{CE}}(\bm{y}^{S},\bm{y})\text{+}\alpha\tau^{2}D_{\mathrm{KLD}}\big{(}\bm{y}^{S},\bm{y}^{T}\big{)}\\\
+\beta\frac{1}{LH}\sum_{l=1}^{L}\sum_{h=1}^{H}\ell_{\mathrm{Attention}}\big{(}\text{\bf
A}^{S}_{lh},\text{\bf A}^{T}_{lh}\big{)}\end{split}$ (2)

where $\alpha$ and $\beta$ are regularization terms and
$\ell_{\text{Attention}}$ computes the loss between the student and teacher
outputs of each self-attention block in $L$ layers and $H$ attention heads per
layer. We also consider two baselines, $\ell_{\text{SDQ}_{\text{Att}}}$ which
is the same as Equation 2 without
$\alpha\tau^{2}D_{\mathrm{KLD}}(\bm{y}^{S},\bm{y}^{T})$ and
$\ell_{\text{SDQ}_{\text{Hid}}}$ which applies the Mean Squared Error (MSE)
loss between the hidden state outputs instead of the attention outputs. The
gradient of $D_{\text{KLD}}(\cdot,\cdot)$ is expressed as $\frac{\partial
D_{\text{KLD}}(\bm{y}^{S}_{i},\bm{y}^{T}_{i})}{\partial\bm{y}^{S}_{i}}=\tau(\bm{y}^{S}_{i}/\tau-\bm{y}^{T}_{i}/\tau)$
and as $\tau\to\infty$, the gradient is approximately
$1/(d_{y}\bm{y}^{S}_{i}-\bm{y}^{T}_{i})$. Similarly, the gradient of the MSE
loss on a single self-attention output in layer $l$ and head $h$ is
$1/n_{lh}(\bm{a}^{S}_{j}-\bm{a}^{T}_{j})$ for a single sample input $x$.
Hence, we see the connection between derivatives between the KLD loss and the
MSE loss when combining them in a single objective. We now move to desribing
how SDQ is used in two QAT methods.

##### Iterative Product Distilled Quantization

We first consider using SDQ with iPQ Stock et al. (2019). This is achieved by
quantizing $m$ subvectors for each $k$ columns of W where a codebook for each
$k$ subvectors is learned to map each subvector to its nearest neighbor in the
learned codebook $C\in\mathbb{R}^{k\times d}$ where $k$ is the number of
codewords. The codebook is updated by minimizing $||\text{\bf
W}-\tilde{\text{\bf W}}||_{2}^{2}=\sum_{i}^{d}||\text{\bf
W}_{[:,i]}-\phi(\text{\bf w}_{[:,i]})||^{2}_{2}$ where $\phi(\cdot)$ is the
quantization function. This objective can be efficiently minimized with the
k-means algorithm and the codewords of each layers are updated with SGD by
averaging the gradients of each assigned block of weights. This is done
iteratively from the bottom layers to the top layers throughout training where
the upper layers are finetuned while the lower layers are progressively being
quantized Stock et al. (2019). When using iPQ with SDQ, omitting the KLD loss
and cross-entropy loss, the objective is
$\ell_{\text{SDQ}_{\text{iPQ}}}=\sum_{l=1}^{L-F}\Big{[}||\text{\bf
W}_{l}-\tilde{\text{\bf
W}}_{l}||^{2}_{2}+\frac{\beta}{L\text{-}F}\sum_{i}^{d}(\text{\bf
A}^{S}_{l,i}-\text{\bf A}^{T}_{l,i})^{2}\Big{]}$ where $F$ is the number of
finetuned layers (non-quantized) at that point in training. Hence, SDQ
progressively quantizes the layers throughout training when used with iPQ.

##### Block-Wise Distilled Quantization Noise

For the majority of our QAT-based experiments we use Quant-Noise Fan et al.
(2020). Quant-Noise is a SoTA QAT method that applies (fake) block-wise
quantization noise at random to each weight matrix. Concretely, blocks of
weights $b_{kl}$ in $\text{\bf W}_{l}$ are chosen at random at a rate $p$ and
quantization noise is added to the chosen blocks.

(a) Quantization Error

(b) F-Norm Elementwise Errors

Figure 2: Dynamic Quantization of InfoXLMBase After Quantization Aware Fine-
Tuning on XNLI.

Student | Quant Method | Teacher | Quant Method | en | Avg.
---|---|---|---|---|---
XLM-R${}_{\text{Base}}$ Conneau et al. | - | - | - | 84.6 | 74.5
XLM-R${}_{\text{Base}}$ | - | - | - | 83.9 | 73.9
InfoXLM${}_{\text{Base}}$ | - | - | - | 84.1 | 74.6
InfoXLM${}_{\text{Base}}$ | PTQ${}_{\text{Dynamic}}$ | - | - | 81.7 | 71.4
XLM-R${}_{\text{Base}}$ | PTQ${}_{\text{Dynamic}}$ | - | - | 80.1 | 72.5
XLM-R${}_{\text{Base}}$ | QNAT | - | - | 82.1 | 70.5
InfoXLM${}_{\text{Base}}$ | QNAT | - | - | 83.7 | 73.0
XLM-R${}_{\text{Base}}$ | QNAT${}_{\text{KLD}}$ | XLM-R${}_{\text{Base}}$ | - | 83.4 | 72.5
XLM-R${}_{\text{Base}}$ | QNAT${}_{\text{KLD}}$ | InfoXLM${}_{\text{Base}}$ | - | 84.4 | 73.3
InfoXLM${}_{\text{Base}}$ | QNAT${}_{\text{KLD}}$ | InfoXLM${}_{\text{Base}}$ | - | 83.9 | 73.6
InfoXLM${}_{\text{Base}}$ | QNAT${}_{\text{Att}}$ | InfoXLM${}_{\text{Base}}$ | - | 84.1 | 73.2
InfoXLM${}_{\text{Base}}$ | QNAT${}_{\text{Att-KLD}}$ | InfoXLM${}_{\text{Base}}$ | - | 84.1 | 73.8
InfoXLM${}_{\text{Base}}$ | QNAT${}_{\text{Att}}$ | InfoXLM${}_{\text{Base}}$ | QNAT-PTQ | 83.3 | 72.1
InfoXLM${}_{\text{Base}}$ | QNAT${}_{\text{Hid}}$ | InfoXLM${}_{\text{Base}}$ | QNAT-PTQ | 81.1 | 70.7
InfoXLM${}_{\text{Base}}$ | QNAT${}_{\text{KLD}}$ | InfoXLM${}_{\text{Base}}$ | QNAT-PTQ | 83.7 | 73.1
InfoXLM${}_{\text{Base}}$ | QNAT${}_{\text{Att-KLD}}$ | InfoXLM${}_{\text{Base}}$ | QNAT-PTQ | 83.9 | 73.4
The best performance obtained are marked in bold.

Table 1: XNLI standard and zero-shot test accuracy using (Fake) Quantization-
Aware Training with INT-8 Post-Training (Real) Dynamic Quantization.

##  4 Empirical Results

(a) Query

(b) Key

(c) Value

(d) Intermediate Out

(e) Attention Output

(f) Output

Figure 3: Dynamic Quantization of InfoXLMBase After Quantization Aware Fine-
Tuning on XNLI.

Student | Teacher | Mem | XNLI | NC | NER | PAWSX | POS | QAM | QADSM | WPR | Avg.
---|---|---|---|---|---|---|---|---|---|---|---
X | - | 1.22 | 73.9 | 83.2 | 83.8 | 89.3 | 79.7 | 68.4 | 68.3 | 73.6 | 77.5
I | - | 1.22 | 74.6 | 83.6 | 85.9 | 89.6 | 79.8 | 68.6 | 68.9 | 73.8 | 78.1
X-PTQ${}_{\text{Dynamic}}$ | - | 0.52 | 71.4 | 81.5 | 82.9 | 87.1 | 76.1 | 66.3 | 65.8 | 68.2 | 74.9
I-PTQ${}_{\text{Dynamic}}$ | - | 0.52 | 72.5 | 81.8 | 83.0 | 87.8 | 75.8 | 66.6 | 66.1 | 68.7 | 75.3
X-QNAT | - | 0.52 | 70.5 | 81.8 | 83.0 | 87.4 | 78.4 | 66.8 | 66.9 | 70.4 | 75.7
I-QNAT | - | 0.52 | 73.0 | 82.1 | 83.1 | 87.8 | 78.0 | 67.2 | 67.2 | 70.8 | 76.2
X-QNAT${}_{\text{KLD}}$ | X | 0.52 | 72.5 | 82.0 | 83.2 | 88.1 | 78.8 | 67.1 | 67.2 | 70.7 | 75.8
X-QNAT${}_{\text{KLD}}$ | I | 0.52 | 73.3 | 82.1 | 82.8 | 88.2 | 78.3 | 67.3 | 67.5 | 70.5 | 75.9
I-QNAT${}_{\text{KLD}}$ | I | 0.52 | 73.6 | 82.6 | 83.1 | 88.4 | 79.5 | 67.6 | 67.9 | 71.8 | 76.8
I-QNAT${}_{\text{Att}}$ | I | 0.52 | 73.2 | 82.4 | 83.0 | 88.3 | 78.3 | 67.8 | 67.7 | 71.7 | 76.6
I-QNAT${}_{\text{Att-KLD}}$ | I | 0.52 | 73.8 | 82.8 | 83.4 | 88.8 | 79.5 | 67.9 | 68.0 | 72.4 | $\underline{\textbf{77.1}}$
I-QNAT${}_{\text{Att}}$ | I${}_{\text{QNAT-PTQ}}$ | 0.52 | 72.1 | 82.1 | 83.1 | 89.2 | 78.8 | 68.0 | 67.8 | 71.9 | 76.6
I-QNAT${}_{\text{Hid}}$ | I${}_{\text{QNAT-PTQ}}$ | 0.52 | 70.7 | 81.9 | 82.4 | 88.8 | 78.4 | 67.3 | 68.0 | 71.4 | 76.1
I-QNAT${}_{\text{KLD}}$ | I${}_{\text{QNAT-PTQ}}$ | 0.52 | 73.1 | 82.3 | 83.0 | 88.4 | 79.2 | 67.6 | 67.9 | 72.1 | 76.7
I-QNAT${}_{\text{Att-KLD}}$ | I${}_{\text{QNAT-PTQ}}$ | 0.52 | 73.4 | 82.5 | 83.3 | 88.9 | 79.6 | 67.9 | 68.2 | 72.6 | $\underline{\mathbf{77.1}}$

Table 2: XGLUE INT-8 Zero-Shot Quantization Results with Post-Training Dynamic
Quantization.

We begin by referring the reader to the supplementary material for the
experimental setup in subsection A.2 and subsection A.3. Before discussing the
main results on XGLUE, we first analyse the mean absolute quantization error
and the Frobenius norm of the elementwise difference in self-attention blocks
between an INT-8 dynamically quantized InfoXLM${}_{\text{Base}}$ and an
unquantized FP-32 InfoXLM${}_{\text{Base}}$ in Figure 2. We see in 2(a) that
the output layer contains the largest mean absolute error across each layer
and highest error variance. In contrast, query, key, value (QKV) parameters
have much smaller error. However, since most of the parameters are found in
the QKV layers, the sum of the quantization error is larger, as seen in 2(b).
This motivates us to focus on the output of the self-attention block when
minimizing quantization errors with our proposed loss in Equation 2 as the
mean error is higher near the output as it accumulates errors from previous
layers in the block. This is also reflected in the parameter distribution of
each layer type across all layers in Figure 3, where the x-axis is the mean
absolute quantization error and the y-axis is the layer indices. We see the
quantization noise is more apparent on the output layer as the Gaussian
distrbutions are non-smooth and have clear jitter effect.

##### XNLI Per Language Results

Table 1 shows the baselines and our SDQ methods applied to
XLM-R${}_{\text{Base}}$ and InfoXLM${}_{\text{Base}}$. Here, both models are
only trained on the English language and hence the remaining languages in the
evaluation set test the zero-shot performance after INT8 quantization (apart
from the first 3 rows that show FP-32 fine-tuned results). On average, we find
that best student networks results are found when distilling using
QNAT${}_{\text{Att-KLD}}$ SDQ with the outputs of an FP-32 teacher for
InfoXLM${}_{\text{Base}}$ at 73.8% test accuracy points, where the original
FP-32 InfoXLM${}_{\text{Base}}$ achieves 74.6%. Additionally we see that
QNAT${}_{\text{Att-KLD}}$ improves over QNAT${}_{\text{KLD}}$ distillation,
indicating that attention output distillation improves the generalization of
the INT-8 student model. We also found that largest performance drops
correspond to languages that less pretraining data and morphologically rich
(Swahili, Urdu, Arabic), while performance in English for the best INT-8
XLM-R${}_{\text{Base}}$ (84.4%) is within 0.2% of the original network (84.6%)
and the best InfoXLM${}_{\text{Base}}$ that uses QNAT${}_{\text{Att-KLD}}$ is
on par with the FP-32 results.

Figure 4: Accuracy versus Number of Bits for XGLUE Tasks with FP-32-16 and
INT-8-4-2 formats.

##### Quantization Results on XGLUE.

We show the per task test performance and the understanding score (i.e average
score) on XGLUE for quantization baselines and our proposed SDQ approaches in
Table 2 (for brevity we denote InfoXLM${}_{\text{Base}}$ as I and
XLM-R${}_{\text{Base}}$). Our proposed QNAT${}_{\text{Att-KLD}}$ achieves the
best average (Avg.) score and per task performance for all tasks, using a
fine-tuned InfoXLM${}_{\text{Base}}$ (XNLI, NC, NER and QAM) and a fine-tuned
InfoXLM${}_{\text{Base}}$ trained with QuantNoise and dynamically quantized
post-training (PAWSX, POS, QAM, QADSM and WPR). We also find that
QNAT${}_{\text{Att-KLD}}$ improves over QNAT${}_{\text{KLD}}$, highlighting
that the attention loss is improving quantized model performance.

##### Performance versus Compression Rate

Figure 4 shows how the performance changes for four approaches, including two
of our proposed objectives (QNAT${}_{\text{KLD}}$ and QNAT${}_{\text{Att-
KLD}}$), when training InfoXLM${}_{\text{Base}}$. As before,
PTQ${}_{\text{dynamic}}$ is a dynamically quantization fine-tuned
InfoXLM${}_{\text{Base}}$ and QNAT-PTQ${}_{\text{dynamic}}$ is the same as
PTQ${}_{\text{dynamic}}$ except fine-tuned also using QuantNoise. Unlike our
previous results, here we apply fake quantization at inference to achieve
compression lower than INT-8 and be comparable to previous work Fan et al.
(2019). We see that performance is generally well maintained up until 8 bits.
However, performance significantly degrades for all quantization methods for 4
and 2 bit weights. We find that QNAT${}_{\text{Att-KLD}}$ maintains higher
performance when compared to the baselines and directly quantizing with no QAT
(PTQ${}_{\text{dynamic}}$) leads to the poorest results, also reflected in
Table 2 results with real dynamic quantization at inference time.

Student | Teacher | XNLI | NC | NER | POS | Avg.
---|---|---|---|---|---|---
- | - | 74.6 | 83.6 | 85.9 | 79.7 | 81.0
iPQ${}_{\text{Scalar}}$ | - | 69.1 | 79.4 | 81.9 | 76.3 | 76.7
iPQ${}_{\text{Scalar-KLD}}$ | Standard | 70.4 | 80.1 | 82.3 | 76.9 | 77.4
iPQ${}_{\text{Scalar-KLD}}$ | iPQ${}_{\text{Scalar}}$ | 70.8 | 80.7 | 82.6 | 79.4 | 78.4
iPQ${}_{\text{Scalar-Att-KLD}}$ | Standard | 72.2 | 80.4 | 82.5 | 77.4 | 78.1
iPQ${}_{\text{Scalar-Att-KLD}}$ | iPQ${}_{\text{Scalar}}$ | 71.3 | 80.4 | 82.9 | 79.6 | 78.6
iPQ${}_{\text{EM}}$ | - | 69.1 | 79.4 | 81.9 | 76.3 | 76.7
iPQ${}_{\text{EM-KLD}}$ | Standard | 70.4 | 80.1 | 82.3 | 76.9 | 77.4
iPQ${}_{\text{EM-KLD}}$ | iPQ${}_{\text{EM}}$ | 72.8 | 81.6 | 82.8 | 79.8 | 79.3
iPQ${}_{\text{EM-Att-KLD}}$ | Standard | 73.2 | 82.3 | 82.7 | 79.1 | 79.3
iPQ${}_{\text{EM-Att-KLD}}$ | iPQ${}_{\text{EM}}$ | 73.1 | 82.5 | 83.0 | 79.2 | 79.5
QNAT | - | 70.5 | 81.8 | 83.3 | 78.4 | 78.5
QNAT${}_{\text{KLD}}$ | Standard | 73.2 | 82.6 | 83.1 | 79.5 | 79.6
QNAT${}_{\text{KLD}}$ | QNAT | 73.1 | 82.3 | 83.0 | 79.2 | 79.4
QNAT${}_{\text{Att-KLD}}$ | Standard | 73.8 | 82.8 | 83.4 | 79.5 | 79.9
QNAT${}_{\text{Att-KLD}}$ | QNAT | 73.4 | 82.5 | 83.3 | 79.6 | 79.7

Table 3: SoTA INT-8 Iterative Product Quantization methods with and without
SDQ for InfoXLM${}_{\text{Base}}$.

##### Ablation with Current QAT Methods

Table 3 shows the results from XGLUE tasks where the first two columns
describe how the student and teacher networks are trained and “Standard”
refers to standard FP-32 fine-tuning. This includes iPQ Stock et al. (2019)
with scalar quantization (iPQ${}_{\text{Scalar}}$), iPQ that uses expectation
maximization to create the codebook during training (iPQ${}_{\text{EM}}$) and
previous results of QuantNoise (QNAT) as a reference point. In this setup, we
only apply the attention loss, $\ell_{\text{Attention}}$, to the layers that
are quantized during iPQ. In all cases, adding that SDQ distillation of the
classification output and the self-attention outputs improves the average
performance.

##  5 Conclusion

In this paper we proposed an attention-based distillation that minimizes
accumulative quantization errors in fine-tuned masked language models. We
identified that most of the quantization errors accumulate at the output of
self-attention blocks and the parameter distribution of the output layer is
effected more by quantization noise. The proposed distillation loss
outperforms baseline distillation without the attention loss and the resulting
INT-8 models are within 1 understanding score points on the XGLUE benchmark
with _real_ quantization post-training. Moreover, fine-tuning the teacher
network with quantization-aware training can further improve student network
performance on some of the tasks. Further compression can be achieved up to
4-bit and 2-bit weights but performance steeply degrades as the network
capacity is drastically reduced coupled with the models having to generalize
to multiple languages it was not trained on.

##  6 Limitations

##### Dataset and Experimental Limitations.

The datasets and tasks we focus on are from the XGLUE benchmark Liang et al.
(2020). The structured prediction tasks, namely Named Entity Recognition (NER)
and Part of Speech (PoS) Tagging, both have a limited number of training
samples at 15k and 25.4k samples respectively. This is due to the difficulty
in annotating on the token level, however it can still be viewed as a
limitation when compared to the remaining sentence-level tasks the majority of
tasks have at least 100k samples.

##### Methodological Limitations.

Below are a list of the main methodological limitations we perceive of our
work:

  * •

Our method requires a teacher model that is already trained on the downstream
task which can then be used to perform knowledge distillation. This is
limiting when there are constraints on the computing resources required to
produce the quantized model.

  * •

We have focused on the problem of reducing accumulative qunatization errors
which become more apparent the deeper a network is. However, this problem is
intuitvely lessened when the model is shallow (e.g 3-4 layers) but perhaps
wider. Hence the results may be less significant if the model is shallower
than what we have experimented in this work.

  * •

By introducing the distillation loss we require an additional regualrization
term $\beta$ to be optimally set, relative to the main distillation loss
$\alpha$. This can be viewed as a potential limitation has it introduced an
additional hyperparameter to be searched to obtain best results on a given
task.

  * •

Lastly, since intermediate layer outputs of the teacher network are required
for self-attention distillation, we have to perform two forward passes during
training. Since standard KLD distillation only requires the output logits, it
is common to store the training data teacher logits, eliminating the need to
perform two forward passes at training data. However, this is not an option
with self-atttention outputs as the storage required offline scales with the
number of self-attention heads, number of layers and the size of the training
data.

##  7 Ethics Statement

Here we briefly discuss some ethical concerns of using such compressed models
in the real world, specifically the two techniques used in this work,
quantization and knowledge distillation. Hooker et al. (2020) have found that
compressed models can amplify existing algorithmic bias and perform very
poorly on a subset of samples while the average out-of-sample accuracy is
maintained close to the uncompressed model. This general finding for pruning
and quantization may be also extrapolated to our work (including
distillation), hence it is important to recognize that our work, much like the
remaining literature on compression, may have ethical concerns with regards to
algorithmic bias and how that effects downstream tasks. However, smaller
models are more cost-efficient and thus become more widely available to the
general public. To summarize, it is important to analyse any aforementioned
bias amplification for subsets of samples for downstream tasks compressed
models are used for.

## References

  * Banner et al. (2018) Ron Banner, Itay Hubara, Elad Hoffer, and Daniel Soudry. 2018.  Scalable methods for 8-bit training of neural networks.  _Advances in neural information processing systems_ , 31.
  * Bondarenko et al. (2021) Yelysei Bondarenko, Markus Nagel, and Tijmen Blankevoort. 2021.  Understanding and overcoming the challenges of efficient transformer quantization.  _arXiv preprint arXiv:2109.12948_.
  * Chi et al. (2020) Zewen Chi, Li Dong, Furu Wei, Nan Yang, Saksham Singhal, Wenhui Wang, Xia Song, Xian-Ling Mao, Heyan Huang, and Ming Zhou. 2020.  Infoxlm: An information-theoretic framework for cross-lingual language model pre-training.  _arXiv preprint arXiv:2007.07834_.
  * Chin et al. (2020) Ting-Wu Chin, Pierce I-Jen Chuang, Vikas Chandra, and Diana Marculescu. 2020.  One weight bitwidth to rule them all.  In _European Conference on Computer Vision_ , pages 85–103. Springer.
  * Conneau et al. (2019) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Unsupervised cross-lingual representation learning at scale.  _arXiv preprint arXiv:1911.02116_.
  * Dettmers et al. (2022) Tim Dettmers, Mike Lewis, Younes Belkada, and Luke Zettlemoyer. 2022.  Gpt3. int8 (): 8-bit matrix multiplication for transformers at scale.  _Advances in Neural Information Processing Systems_ , 35:30318–30332.
  * Faghri et al. (2020) Fartash Faghri, Iman Tabrizian, Ilia Markov, Dan Alistarh, Daniel M Roy, and Ali Ramezani-Kebrya. 2020.  Adaptive gradient quantization for data-parallel sgd.  _Advances in neural information processing systems_ , 33:3174–3185.
  * Fan et al. (2019) Angela Fan, Edouard Grave, and Armand Joulin. 2019.  Reducing transformer depth on demand with structured dropout.  _arXiv preprint arXiv:1909.11556_.
  * Fan et al. (2020) Angela Fan, Pierre Stock, Benjamin Graham, Edouard Grave, Rémi Gribonval, Herve Jegou, and Armand Joulin. 2020.  Training with quantization noise for extreme model compression.  _arXiv preprint arXiv:2004.07320_.
  * Hinton et al. (2015) Geoffrey Hinton, Oriol Vinyals, and Jeff Dean. 2015.  Distilling the knowledge in a neural network.  _arXiv preprint arXiv:1503.02531_.
  * Hooker et al. (2020) Sara Hooker, Nyalleng Moorosi, Gregory Clark, Samy Bengio, and Emily Denton. 2020\.  Characterising bias in compressed models.  _arXiv preprint arXiv:2010.03058_.
  * Jacob et al. (2018) Benoit Jacob, Skirmantas Kligys, Bo Chen, Menglong Zhu, Matthew Tang, Andrew Howard, Hartwig Adam, and Dmitry Kalenichenko. 2018.  Quantization and training of neural networks for efficient integer-arithmetic-only inference.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2704–2713.
  * Kim et al. (2019) Jangho Kim, Yash Bhalgat, Jinwon Lee, Chirag Patel, and Nojun Kwak. 2019.  Qkd: Quantization-aware knowledge distillation.  _arXiv preprint arXiv:1911.12491_.
  * Kim et al. (2020) Jangho Kim, KiYoon Yoo, and Nojun Kwak. 2020.  Position-based scaled gradient for model quantization and sparse training.  _arXiv preprint arXiv:2005.11035_.
  * Kim et al. (2021) Sehoon Kim, Amir Gholami, Zhewei Yao, Michael W Mahoney, and Kurt Keutzer. 2021\.  I-bert: Integer-only bert quantization.  In _International conference on machine learning_ , pages 5506–5518. PMLR.
  * Krishnamoorthi (2018) Raghuraman Krishnamoorthi. 2018.  Quantizing deep convolutional networks for efficient inference: A whitepaper.  _arXiv preprint arXiv:1806.08342_.
  * Liang et al. (2020) Yaobo Liang, Nan Duan, Yeyun Gong, Ning Wu, Fenfei Guo, Weizhen Qi, Ming Gong, Linjun Shou, Daxin Jiang, Guihong Cao, et al. 2020.  Xglue: A new benchmark datasetfor cross-lingual pre-training, understanding and generation.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 6008–6018.
  * Mishra and Marr (2017) Asit Mishra and Debbie Marr. 2017.  Apprentice: Using knowledge distillation techniques to improve low-precision network accuracy.  _arXiv preprint arXiv:1711.05852_.
  * Neill (2020) James O’ Neill. 2020.  An overview of neural network compression.  _arXiv preprint arXiv:2006.03669_.
  * NVIDIA (2017) Tesla NVIDIA. 2017.  Nvidia tesla v100 gpu architecture.  _Tesla NVIDIA_.
  * Polino et al. (2018) Antonio Polino, Razvan Pascanu, and Dan Alistarh. 2018.  Model compression via distillation and quantization.  _arXiv preprint arXiv:1802.05668_.
  * Prato et al. (2019) Gabriele Prato, Ella Charlaix, and Mehdi Rezagholizadeh. 2019.  Fully quantized transformer for machine translation.  _arXiv preprint arXiv:1910.10485_.
  * Stock et al. (2019) Pierre Stock, Armand Joulin, Rémi Gribonval, Benjamin Graham, and Hervé Jégou. 2019.  And the bit goes down: Revisiting the quantization of neural networks.  _arXiv preprint arXiv:1907.05686_.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  _arXiv preprint arXiv:1706.03762_.
  * Wang et al. (2018) Naigang Wang, Jungwook Choi, Daniel Brand, Chia-Yu Chen, and Kailash Gopalakrishnan. 2018.  Training deep neural networks with 8-bit floating point numbers.  _Advances in neural information processing systems_ , 31.
  * Zafrir et al. (2019) Ofir Zafrir, Guy Boudoukh, Peter Izsak, and Moshe Wasserblat. 2019.  Q8bert: Quantized 8bit bert.  In _2019 Fifth Workshop on Energy Efficient Machine Learning and Cognitive Computing-NeurIPS Edition (EMC2-NIPS)_ , pages 36–39. IEEE.
  * Zhao et al. (2019) Ritchie Zhao, Yuwei Hu, Jordan Dotzel, Chris De Sa, and Zhiru Zhang. 2019.  Improving neural network quantization without retraining using outlier channel splitting.  In _International conference on machine learning_ , pages 7543–7552. PMLR.
  * Zhou et al. (2017) Aojun Zhou, Anbang Yao, Yiwen Guo, Lin Xu, and Yurong Chen. 2017.  Incremental network quantization: Towards lossless cnns with low-precision weights.  _arXiv preprint arXiv:1702.03044_.

##  Appendix A Supplementary Material

###  A.1 Self-Attention in Transformers

Consider a dataset $D=\\{(X_{i},y_{i})\\}_{i=1}^{m}$ for $D\in\mathcal{D}$ and
a sample $s:=(X,y)$ where the sentence $X:=(x_{1},\ldots x_{n})$ with $n$
being the number of words $x\in X$. We can represent a word as an input
embedding $\bm{x}_{w}\in\mathbb{R}^{d}$, which has a corresponding target
vector $\bm{y}$. In the pre-trained transformer models we use, $X_{i}$ is
represented by 3 types of embeddings; word embeddings ($\text{\bf
X}_{w}\in\mathbb{R}^{n\times d}$), segment embeddings ($\text{\bf
X}_{s}\in\mathbb{R}^{n\times d}$) and position embeddings ($\text{\bf
X}_{p}\in\mathbb{R}^{n\times d}$), where $d$ is the dimensionality of each
embedding matrix. The self-attention block in a transformer mainly consists of
three sets of parameters: the query parameters $\text{\bf
Q}\in\mathbb{R}^{d\times l}$, the key parameters $\text{\bf
K}\in\mathbb{R}^{d\times l}$ and the value parameters $\text{\bf
V}\in\mathbb{R}^{d\times o}$. For 12 attention heads (as in
XLM-R${}_{\text{Base}}$ and InfoXLM${}_{\text{Base}}$), we express the forward
pass as follows:

$\displaystyle\overrightarrow{\text{\bf X}}=\text{\bf X}_{w}+\text{\bf
X}_{s}+\text{\bf X}_{p}$ (3) $\displaystyle\overrightarrow{\text{\bf
Z}}:=\bigoplus_{i=1}^{12}\text{softmax}\big{(}\overrightarrow{\text{\bf
X}}\text{\bf Q}_{(i)}\text{\bf K}^{T}_{{(i)}}\overrightarrow{\text{\bf
X}}^{T}\big{)}\overrightarrow{\text{\bf X}}\text{\bf V}_{(i)}$ (4)
$\displaystyle\overrightarrow{\mathbb{Z}}=\text{Feedforward}(\text{LayerNorm}(\overrightarrow{\text{\bf
Z}}+\overrightarrow{\text{\bf X}}))$ (5)
$\displaystyle\overleftarrow{\mathbb{Z}}=\text{Feedforward}(\text{LayerNorm}(\overleftarrow{\text{\bf
Z}}+\overleftarrow{\text{\bf X}}))$ (6)

The last hidden representations of both directions are then concatenated
$\mathbb{Z}^{\prime}:=\overleftarrow{\mathbb{Z}}\bigoplus\overrightarrow{\mathbb{Z}^{\prime}}$
and projected using a final linear layer $\text{\bf W}\in\mathbb{R}^{d}$
followed by a sigmoid function $\sigma(\cdot)$ to produce a probability
estimate $\hat{y}$, as shown in (7). Words from (step-3) that are used for
filtering the sentences are masked using a [PAD] token to ensure the model
does not simply learn to correctly classify some samples based on the
association of these tokens with counterfacts. A linear layer is then fine-
tuned on top of the hidden state, $\bm{h}_{X,\texttt{[CLS]}}$ emitted
corresponding to the [CLS] token. This fine-tunable linear layer is then used
to predict whether the sentence is counterfactual or not, as shown in Equation
7, where $\mathcal{B}\subset D$ is a mini-batch and $\mathcal{L}_{ce}$ is the
cross-entropy loss.

$\displaystyle\mathcal{L}_{ce}:=\frac{1}{|\mathcal{B}|}\sum_{(X,y)\in\mathcal{B}}\bm{y}\log\big{(}\sigma(\bm{h}_{X,\texttt{[CLS]}}\cdot\text{\bf
W})\big{)}$ (7)

Configurations We use XLM-R${}_{\text{Base}}$ and InfoXLM${}_{\text{Base}}$,
which uses 12 Transformer blocks, 12 self-attention heads with a hidden size
of 768. The default size of 512 is used for the sentence length and the
sentence representation is taken as the final hidden state of the first [CLS]
token.

###  A.2 Experimental Setup and Hardware Details

Below describes the experimental details, including model, hyperparameter and
quantization details. We choose modestly sized cross-lingual language models
as the basis of our experiments, namely XLM-R${}_{\text{Base}}$ Conneau et al.
(2019) and InfoXLM${}_{\text{Base}}$ Chi et al. (2020), both approximately
1.1GB in memory and these pretrained models are retrieved from the huggingface
model hub.

We choose both XLM-R${}_{\text{Base}}$ and InfoXLM${}_{\text{Base}}$ because
they are relatively small Transformers and are required to generalized to
languages other than the language used for fine-tuning. Hence, we begin from a
point that model are already relatively difficult to compress and are further
motivated by the findings that larger overparameterized networks suffer less
from PTQ to 8-bit integer format and lower (Jacob et al., 2018;
Krishnamoorthi, 2018).

For both XLM-R${}_{\text{Base}}$ and InfoXLM${}_{\text{Base}}$ the hyper-
parameters are set as follows: 768 hidden units, 12 heads, GELU activation, a
dropout rate of 0.1, 512 max input length, 12 layers in encoder. The Adam
Optimizer with a linear warm-up Vaswani et al. (2017) and set the learning
rate to 2e-5 for most tasks. For all sentence classification tasks the batch
size is set to 32 and we fine-tune with 10 epochs. For POS Tagging and NER, we
fine-tune with 20 epochs and set the learning rate to 2e-5. We select the
model with the best average results on the development sets of all languages.
For SDQ-based models, we report the best performing model for
$\alpha\in[0.1,0.2,0.5,0.8]$ and $\beta\in[10,100,200,500]$. All experiments
are carried out on Tesla V100-SXM2 32 Gigabyte GPUs NVIDIA (2017) with no
constraint on GPU hours used on these machines. In all reported results, we
report the best (max) result from 8-16 different runs when searching for
$\alpha$ and $\beta$ depending on each particular task.

###  A.3 Model Configuration and Hyperparameter Settings

XLM-R${}_{\text{Base}}$ and InfoXLM${}_{\text{Base}}$ uses 12 Transformer
blocks, 12 self-attention heads with a hidden size of 768. The default size of
512 is used for the sentence length and the sentence representation is taken
as the final hidden state of the first [CLS] token. A fine-tuned linear layer
W is used on top of both models, which is fed to through a softmax function
$\sigma$ as $p(c|h)=\sigma(\text{\bf W}\bm{h})$ where $c$ is used to calibrate
the class probability estimate and we maximize the log-probability of
correctly predicting the ground truth label.

Table 4 shows the pretrained model configurations that were already predefined
before our experiments. The number of (Num.) hidden groups here are the number
of groups for the hidden layers where parameters in the same group are shared.
The intermediate size is the dimensionality of the feed-forward layers of the
the Transformer encoder. The ‘Max Position Embeddings’ is the maximum sequence
length that the model can deal with.

Hyperparameters | XLM-R${}_{\text{Base}}$ | InfoXLM${}_{\text{Base}}$
---|---|---
Vocab Size | 250002 | 250002
Max Pos. Embeddings | 514 | 514
Hidden Size | 3072 | 3072
Encoder Size | 768 | 768
Num. Hidden Layers | 12 | 12
Num. Hidden Groups | 1 | 1
Num. Attention Heads | 12 | 12
Hidden Activations | GeLU | GeLU
Layer Norm. Epsilon | $10^{-12}$ | $10^{-12}$
Fully-Connected Dropout Prob. | 0.1 | 0.1
Attention Dropout Prob. | 0 | 0

Table 4: Model Hyperparameter Settings

We now detail the hyperparameter settings for transformer models and the
baselines. We note that all hyperparameter settings were performed using a
manual search over development data.

####  A.3.1 Transformer Model Hyperparameters

We did not change the original hyperparameter settings that were used for the
original pretraining of each transformer model. The hyperparameter settings
for these pretrained models can be found in the class arguments python
documentation in each configuration python file in the
https://github.com/huggingface/transformers/blob/master/src/transformers/ e.g
configuration_.py.For fine-tuning transformer models, we manually tested
different combinations of a subset of hyperparameters including the learning
rates $\\{50^{-4},10^{-5},50^{-5}\\}$, batch sizes $\\{16,32,128\\}$, warmup
proportion $\\{0,0.1\\}$ and $\epsilon$ which is a hyperparameter in the
adaptive momentum (adam) optimizer. Please refer to the huggingface
documentation at https://github.com/huggingface/transformers for further
details on each specific model e.g at
https://github.com/huggingface/transformers/blob/master/src/transformers/modeling_roberta.py,
and also for the details of the architecture that is used for sentence
classification and token classification.


# Solving optimal control problems with non-smooth solutions using an
integrated residual method and flexible mesh

Lucian Nita1, Eric C. Kerrigan1, Eduardo M. G. Vila1 and Yuanbo Nie2 This work
has received funding from the Engineering and Physical Sciences Research
Council under a Doctoral Training Grant (reference number: EP/T51780X/1)1
Department of Electrical & Electronic Engineering, Imperial College London,
SW7 2AZ London, UK, email: {lucian.nita16, e.kerrigan<EMAIL_ADDRESS>Department of Automatic Control and Systems Engineering, University of
Sheffield, S1 3JD Sheffield, UK, email<EMAIL_ADDRESS>
###### Abstract

Solutions to optimal control problems can be discontinuous, even if all the
functionals defining the problem are smooth. This can cause difficulties when
numerically computing solutions to these problems. While conventional
numerical methods assume state and input trajectories are continuous and
differentiable or smooth, our method is able to capture discontinuities in the
solution by introducing time-mesh nodes as decision variables. This allows one
to obtain a higher accuracy solution for the same number of mesh nodes
compared to a fixed time-mesh approach. Furthermore, we propose to first solve
a sequence of suitably-defined least-squares problems to ensure that the error
in the dynamic equation is below a given tolerance. The cost functional is
then minimized subject to an inequality constraint on the dynamic equation
residual. We demonstrate our implementation on an optimal control problem that
has a chattering solution. Solving such a problem is difficult, since the
solution involves infinitely many switches of decreasing duration. This
simulation shows how the flexible mesh is able to capture discontinuities
present in the solution and achieve superlinear convergence as the number of
mesh intervals is increased.

##  I INTRODUCTION

Solving a sequence of constrained optimal control problems (OCPs) in real-time
is a very powerful technique typically used in model predictive control (MPC).
However, it can also be challenging in practice to reliably compute a solution
since the continuous-time, infinite-dimensional OCP is solved using finite-
dimensional numerical solvers. Since we are solving a discretized version of
the original problem, the obtained solution is not guaranteed to be feasible
for the original continuous-time OCP. Additionally, it can also be difficult
to preserve the accuracy of the solution and obtain superlinear convergence as
the number of mesh nodes is increased, especially in problems with
discontinuous solutions.

The most commonly used direct transcription method for solving optimal control
problems is direct collocation, which is considered to be the current state of
the art [13]. While collocation has the advantage of being able to handle
complex dynamical models, collocation has the fundamental drawback of not
guaranteeing an acceptable accuracy in between the collocation points.

The idea of using integrated residuals as part of the transcription process
overcomes some of the limitations of collocation [9, 10]. Compared to
classical time-marching schemes (shooting methods) or point-wise residual
minimisation (collocation), integrated residual methods have the benefit of
producing a solution with a more uniform error over the whole time domain.

State-of-the-art mesh refinement methods are hp-adaptive methods [7, 12]. Some
advanced methods have discontinuity detection schemes, but most current
refinement strategies rely on knowing beforehand whether the solution will be
discontinuous and cannot provide an efficiency comparable to the continuous
case for general problems. By adding a flexible mesh, as in this paper, it is
possible to develop a method that has similar convergence properties for
problems with discontinuous or continuous solutions.

This paper extends the work in [11], which used the integrated residual method
for solving differential equations, feasibility problems and constraint
satisfaction problems with a flexible mesh. The central contribution of this
paper is to extend such integrated residual methods to the solution of OCPs.
In the transcription process, a flexible time-mesh will be introduced in order
to achieve superlinear convergence for discontinuous problems during the mesh
refinement phase. Moreover, the proposed algorithm is able to solve difficult
problems to a user-defined accuracy. A numerical example shows how our method
performs on a control problem with an optimal chattering solution.

In Section II, we introduce the optimal control problem formulation. In
Section III-A we present the integrated residual method for transcribing the
constrained optimal control problem. Section III-B describes the concept of a
flexible time mesh and how this method can improve convergence. Section III-C
presents an algorithm for solving an OCP to a user-specified accuracy. Section
IV demonstrates how the proposed algorithm from Section III-C can be used
jointly with a flexible mesh scheme to solve an optimal control problem and
construct a Pareto front between solution accuracy and a lower bound on the
optimal cost. Section V provides a summary of the main findings presented in
this paper and discusses potential improvements and future works.

##  II Problem definition

The objective functional of many optimal control and estimation problems can
be written in the general Bolza form

$\displaystyle\min_{x(\cdot),u(\cdot)}$
$\displaystyle\phi(x(t_{0}),x(t_{f}),t_{0},t_{f})+\int_{t_{0}}^{t_{f}}{L(x(t),u(t),t)}\mathop{}\\!\mathrm{d}t$
(1a) s.t. $\displaystyle F(\dot{x}(t),x(t),u(t),t)=0\quad\forall
t\in[t_{0},t_{f}],$ (1b) $\displaystyle G(\dot{x}(t),x(t),u(t),t)\leq
0\quad\forall t\in[t_{0},t_{f}],$ (1c) where
$x:\mathbb{R}\rightarrow\mathbb{R}^{N_{x}}$ are the state variables and
constrained to be continuous,
$\dot{x}:\mathbb{R}\rightarrow\mathbb{R}^{N_{x}}$ are the time derivatives of
the state $x$, and $u:\mathbb{R}\rightarrow\mathbb{R}^{N_{u}}$ are the control
inputs. The function
$F:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{u}}\times\mathbb{R}\rightarrow\mathbb{R}^{N_{F}}$,
which contains the dynamical model of the system, defines a set of $N_{F}$
equality constraints that have to be satisfied by the controlled system.
$G:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{u}}\times\mathbb{R}\rightarrow\mathbb{R}^{N_{g}}$
defines $N_{g}$ path inequality constraints.
$\phi:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$
is the Mayer cost functional, also called the boundary cost, with
$t_{0}\in\mathbb{R}$ and $t_{f}\in\mathbb{R}$ being the initial and final
times, respectively.
$L:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{u}}\times\mathbb{R}\rightarrow\mathbb{R}$
is the Lagrange cost functional, also known as the path cost. Additionally,
the problem may have one or more boundary constraints of the form
$\displaystyle\Psi_{E}(x(t_{0}),x(t_{f}),t_{0},t_{f})$ $\displaystyle=0,$ (1d)
$\displaystyle\Psi_{I}(x(t_{0}),x(t_{f}),t_{0},t_{f})$ $\displaystyle\leq 0,$
(1e)

where
$\Psi_{E}:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}^{N_{E}}$
are the boundary equality constraints, and
$\Psi_{I}:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}^{N_{I}}$
are the boundary inequality constraints.

##  III Solution method

In most real-time control applications we are heavily constrained by the
computational time. As a result, solving the entire problem (1) using direct
collocation has two fundamental drawbacks. Firstly, the designer is not able
to control the solution accuracy over the entire time interval without a
posteriori computing the error and conducting mesh refinement procedures.
Consequently, existing state-of-the-art methods may fail to ensure constraint
satisfaction. Secondly, existing schemes cannot terminate early and return the
best feasible solution that was achieved in the given amount of computational
time. Since one often wants to focus on fast constraint satisfaction, we
propose to initially solve a feasibility problem and refine the mesh until the
dynamic constraints are satisfied to a given accuracy. We will then use the
obtained solution as an initial guess to the optimal control solver, which
optimizes a transcribed version of the original problem (1).

###  III-A Integrated residual transcription

In the transcription process the infinite-dimensional OCP (1) has to be
converted into a finite-dimensional nonlinear programming problem (NLP). In
order to achieve this, the state $x(\cdot)$ and input $u(\cdot)$ trajectories
need to be parametrized by a finite number of decision variables $s_{i}^{j}$
and $c_{i}^{j}$ where the subscript $i$ denotes the interval number and $j$
denotes the index of the nodal point in interval $i$, as will be described
later. Using a linear combination of these decision variables, approximating
functions $\tilde{x}:\mathbb{R}\rightarrow\mathbb{R}^{N_{x}}$ and
$\tilde{u}:\mathbb{R}\rightarrow\mathbb{R}^{N_{u}}$ can be constructed.

Before aiming to minimize the objective, in most applications it is critical
to ensure the constraints are satisfied to a user-defined accuracy. For this
purpose we will introduce an error metric $\epsilon_{R}\in\mathbb{R}$ defined
as

$\epsilon_{R}:=\frac{1}{(t_{f}-t_{0})\cdot
N_{F}}\int_{t_{0}}^{t_{f}}{\left\lVert
F(\dot{\tilde{x}}(t),\tilde{x}(t),\tilde{u}(t),t)\right\rVert^{2}_{2}}\mathop{}\\!\mathrm{d}t$
(2)

based on the integral of the 2-norm squared of the dynamic equation residual.
The residual $\left\lVert
F(\dot{\tilde{x}}(t),\tilde{x}(t),\tilde{u}(t),t)\right\rVert^{2}_{2}$
indicates how well the numerical solution
$(\tilde{x}(\cdot),\tilde{u}(\cdot))$ satisfies the dynamic constraint (1b)
over the whole time interval. In contrast with direct collocation that
enforces constraint (1b) exactly, but only at a finite number of nodes called
collocation points, our method uses quadrature rules to integrate the residual
over the whole interval $[t_{0},t_{f}]$, thus guaranteeing a certain level of
accuracy in between the collocation points. Note also that the above error
metric is a scaled version of the integrated residual where the scaling factor
$\frac{1}{(t_{f}-t_{0})\cdot N_{F}}$ is introduced to average out the residual
over the interval $[t_{0},t_{f}]$ and over all components of the dynamics
function $F$.

Lagrange polynomial basis functions are often used to express the
approximating functions $\tilde{x}$ and $\tilde{u}$ [2, Sect. 1.17.1]. The
possible solution space is defined by the basis functions used to represent
approximation functions $(\tilde{x}(\cdot),\tilde{u}(\cdot))$. As a
consequence, the exact solution $(x(\cdot),u(\cdot))$ may not be representable
in that solution space, which implies that an exact representation of the
constraint (1b) can never be achieved in finite time (the integrated residual
$\epsilon_{R}$ can asymptotically converge to zero only in the limit as the
time-mesh is refined and the number of discretization points is increased).

To reduce the approximation error (as quantified by $\epsilon_{R}$) there are
two fundamental refinement strategies:

  * •

h-refinement involves splitting the entire time domain $[t_{0},t_{f}]$ into
$N$ subdomains, i.e. subintervals $[t_{i},t_{i+1}]$ such that

$\displaystyle\mathcal{T}_{i}:=[t_{i},t_{i+1}]\subset[t_{0},t_{f}],\ \forall
i\in\\{0,\dots,N-1\\},$ (3a)
$\displaystyle\cup_{i=0}^{N-1}[t_{i},t_{i+1}]=[t_{0},t_{f}],$ (3b)
$\displaystyle t_{i}<t_{i+1},\ \forall i\in\\{0,\dots,N-1\\}$ (3c)

where $t_{N}=t_{f}$. The refinement variable is therefore the number of
subdomains $N$.

  * •

p-refinement relies on constructing a polynomial approximation $\chi_{i}$ of
degree $a$ inside each subdomain $[t_{i},t_{i+1}]$ such that for all
$i\in\\{0,\dots,N-1\\}$:

$\displaystyle\tilde{x}(t)$ $\displaystyle:=\chi_{i}(t),\ \forall
t\in[t_{i},t_{i+1}]$ (4a) $\displaystyle\chi_{i}(t)$
$\displaystyle=\frac{\sum_{j=0}^{a}{\frac{w_{i}^{j}}{t-\tau_{i}^{j}}\cdot
s_{i}^{j}}}{\sum_{j=0}^{a}\frac{w_{i}^{j}}{t-\tau_{i}^{j}}},$ (4b)

where $s_{i}^{j}=\chi_{i}(\tau_{i}^{j})=\tilde{x}(\tau_{i}^{j})$ are NLP
decision variables, $w_{i}^{j}$ are polynomial weights and $\tau_{i}^{j}$ are
polynomial nodes [1]. In this case, polynomial refinement means increasing the
polynomial degree $a$. Note a similar expression for $\tilde{u}(\cdot)$ can be
derived with $b$ denoting the polynomial degree of $\tilde{u}(\cdot)$.

Note that these elementary methods can both be used during the mesh refinement
process leading to the so-called hp-type refinement method.

To enforce state continuity at mesh nodes $t_{i}$, the additional constraints

$\displaystyle\tilde{\chi}_{i}(t_{i+1})$
$\displaystyle=\tilde{\chi}_{i+1}(t_{i+1}),\quad\forall
i\in\\{0,\dots,N-2\\},$ (5)

are enforced by using the same variable $s_{i}^{a}=s_{i+1}^{0}$ to represent
both $\tilde{\chi}_{i}(t_{i+1})$ and $\tilde{\chi}_{i+1}(t_{i+1}),\forall
i\in\\{0,\dots,N-2\\}$.

###  III-B Residual minimization problem: Improving accuracy to ensure
feasibility

To efficiently solve feasibility and control problems with discontinuous
solutions, which are otherwise difficult to solve, we will use an integrated
residual method to tackle the dynamic constraints. The idea is similar to what
[3], [8] have proposed for solving differential equations and what has been
used in [11] for solving dynamic feasibility problems.

The first step of our approach is to solve a feasibility problem that aims to
satisfy constraints (1b)–(1e) to a given tolerance. This feasibility problem
is converted into a minimisation problem that minimizes the integrated
residual of the dynamics model $\epsilon_{R}$ below a user-specified value.

In numerical simulations, integrals from (1a) and (2) have to be approximated
using a $Q$-point Gaussian quadrature rule. Since $F$ is a general nonlinear
function, the approximation of the above integrals will not be exact. Apart
from the residual error $\epsilon_{R}$ appearing as a result of the
discretization, another numerical error is introduced, namely the quadrature
error

$\epsilon_{Q}\coloneqq\Big{|}\epsilon_{R}-\sum_{i=0}^{N-1}{\sum_{k=1}^{Q}{\sigma_{i}^{k}\cdot\left\lVert
F(\dot{\tilde{x}}(\rho_{i}^{k}),\tilde{x}(\rho_{i}^{k}),\tilde{u}(\rho_{i}^{k}),\rho_{i}^{k})\right\rVert^{2}_{2}\Big{|}}}$
(6)

where $\sigma_{i}^{k}$ for $k\in\\{1,\dots,Q\\}$ are the $Q$ quadrature
weights associated with the integration interval $[t_{i},t_{i+1}]$,
appropriately scaled by $N_{F}$ and interval length to include the initial
factors in (2), while $\rho_{i}^{k}$ are the quadrature nodes for the interval
$[t_{i},t_{i+1}]$.

In order to validate the obtained solution, we need to check whether
$\epsilon_{Q}$ is sufficiently small by recomputing the integrals with a
higher quadrature order. If the difference between the new value and the
solution obtained from the optimization problem is above a certain tolerance
$\varepsilon_{quad,tol}$, the problem needs to be resolved using a higher
value for $Q$.

We rely on mesh refinement to select appropriate values for $N$, $a$ and $b$.
Note however that conventional mesh refinement strategies applied to a fixed
time-mesh with nodes at predefined locations may not always achieve
superlinear converge to the solution as the number of nodes is increased.
Consider for example the case when a discontinuity in the solution $u(\cdot)$
is located in the interval $(t_{i},t_{i+1})$. In this case, a numerical
approximation of this discontinuous function is obtained using a continuous
polynomial basis (as described in Section III-A). In general, unless a mesh
node is located exactly at the point of discontinuity, a Gibbs phenomenon can
occur when interpolating a discontinuous function with a continuous one. This
leads to interpolation overshoots that cannot be eliminated in general by mesh
refinement schemes and will cause the error to plateau and not decrease beyond
a certain level.

In order to achieve superlinear convergence in cases where discontinuities are
present and produce an accurate solution, we propose including mesh points
$t_{i}$ as decision variables in the NLP formulation. As a result, time nodes
are allowed to move towards regions non-smoothness.

Recall that standard direct collocation methods compute an integral of the
residual only after the NLP has been solved [2]; they do not directly
constrain the integral of the residual while computing a solution to the NLP.
It follows that introducing mesh nodes as decision variables in standard
collocation methods can result in less accurate solutions than those with
fixed nodes, unless care is taken. This argument also motivates the
interdependence between the use of a flexible mesh and the integrated residual
transcription method proposed here.

In an ideal scenario where no quadrature error is present, nodes can be
allowed to move freely in the domain according to (3). However, since
quadrature error increases as the intervals expand, we still need to constrain
the allowed flexibility of the nodes. For a fixed parameter $\phi\in[0,1)$ we
impose upper and lower bounds on the interval length

$\displaystyle t_{i+1}-t_{i}$
$\displaystyle\leq(1+\phi)\cdot\frac{t_{f}-t_{0}}{N},\ \forall
i\in\\{1,\dots,N-1\\},$ (7a) $\displaystyle t_{i+1}-t_{i}$
$\displaystyle\geq(1-\phi)\cdot\frac{t_{f}-t_{0}}{N},\ \forall
i\in\\{1,\dots,N-1\\}.$ (7b)

Since $\phi\in[0,1)$ by definition, the order of nodes is preserved and
intervals do not overlap, as required by (3c).

For a better sparsity structure of the Hessian along with the ease of
implementing continuity constraints (5) and boundary constraints (1d), (1e)
the decision vector
$z\in\mathbb{R}^{N\cdot(N_{x}\cdot(a)+N_{u}\cdot(b+1)+1)+N_{x}+1}$ is ordered
as

$z:=(s_{0}^{0},t_{0},s_{0}^{1},\dots,s_{0}^{a-1},c_{0}^{0},\dots,c_{0}^{b},s_{0}^{a},t_{1},s_{1}^{1}\dots,s_{N-1}^{a},t_{N}).$
(8)

Hence, the minimum residual solution along with the optimal node locations can
be computed from the optimization problem

$\displaystyle\min_{z}\quad$
$\displaystyle\sum_{i=0}^{N-1}{\sum_{k=1}^{Q}{\sigma_{i}^{k}\cdot\left\lVert
F(\dot{\tilde{x}}(\rho_{i}^{k}),\tilde{x}(\rho_{i}^{k}),\tilde{u}(\rho_{i}^{k}),\rho_{i}^{k})\right\rVert^{2}_{2}}}$
(9a) s.t. $\displaystyle
G(\dot{\tilde{x}}(\tau_{i}^{j}),\tilde{x}(\tau_{i}^{j}),\tilde{u}(\tau_{i}^{j}),\tau_{i}^{j})\leq
0\quad\forall\tau_{i}^{j}\in\tau,$ (9b)
$\displaystyle\Psi_{E}(s_{0}^{0},s_{N-1}^{a},t_{0},t_{f})=0,$ (9c)
$\displaystyle\Psi_{I}(s_{0}^{0},s_{N-1}^{a},t_{0},t_{f})\leq 0,$ (9d)
$\displaystyle(1-\phi)\frac{t_{f}-t_{0}}{N}\leq
t_{i+1}-t_{i}\leq(1+\phi)\frac{t_{f}-t_{0}}{N},$ (9e)
$\displaystyle\tilde{\chi}_{i}(t_{i+1})=\tilde{\chi}_{i+1}(t_{i+1})\quad\forall
i\in\\{0,\dots,N-2\\},$ (9f)

where path inequality constraints (9b) are implemented at the support time
points $\tau_{i}^{j}$. Note that, since time-mesh nodes are added in the
decision vector, quadrature points $\rho_{i}^{k}$ and internal supports
$\tau_{i}^{j}$ become functions of $t_{i}$ and $t_{i+1}$. Hence, these values
need to be shifted and scaled accordingly.

###  III-C Cost minimization problem: From constrained control to optimal
control

In this paper we will focus on h-refinement. Starting from a coarse mesh
(small $N$) problem (9) is solved repeatedly until a desirable user-defined
tolerance $\epsilon_{tol}$ on (2) has been reached. Even if the convergence is
superlinear with respect to the number of subdomains $N$, the performance can
further be improved by providing a good initial guess obtained from
interpolating the solution obtained from the previous solution with a coarser
mesh.

After the desired tolerance has been reached, the mesh parameters $N$ and $a$
are fixed and the cost functional is minimized by solving

$\displaystyle\min_{z}\,$
$\displaystyle\phi(x(t_{0}),x(t_{f}),t_{0},t_{f})+\sum_{i=0}^{N-1}{\sum_{k=1}^{Q}{\sigma_{i}^{k}L(x(\rho_{i}^{k}),u(\rho_{i}^{k}),\rho_{i}^{k})}}$
(10a) s.t.
$\displaystyle\sum_{i=0}^{N-1}{\sum_{k=1}^{Q}{\sigma_{i}^{k}\cdot\left\lVert
F(\dot{\tilde{x}}(\rho_{i}^{k}),\tilde{x}(\rho_{i}^{k}),\tilde{u}(\rho_{i}^{k}),\rho_{i}^{k})\right\rVert^{2}_{2}}}\leq\epsilon_{tol}$
(10b)
$\displaystyle\eqref{eq:equation9B},\eqref{eq:equation9C},\eqref{eq:equation9D},\eqref{eq:equation9E},\eqref{eq:equation9F},$
(10c)

Since we have successfully solved (9), we know the dynamics constraint (10b)
should be feasible. We also have an upper bound on the cost and a sufficiently
good initial guess needed to efficiently warm-start problem (10).

The proposed algorithm for solving optimal control problems using integrated
residual transcription method is outlined in Algorithm 1.

0: $\epsilon_{tol}$, $\varepsilon_{quad,tol}$ and initial values for $N$, $a$,
$b$, $Q$ and $z^{*}$.

1: repeat

2: $z_{0}\leftarrow z^{*}$

3: $z^{*}\leftarrow\arg\min_{z}\eqref{eq:equation9A}$ s.t. (9b), (9c), (9d),
(9e), (9f)

4: if $\epsilon_{Q}\leq\varepsilon_{quad,tol}$ then

5: $\epsilon_{R}\leftarrow\min_{z}\eqref{eq:equation9A}$ s.t. (9b), (9c),
(9d), (9e), (9f)

6: increase $N$

7: else

8: increase $Q$

9: end if

10: until $\epsilon_{R}\leq\epsilon_{tol}$

11: $z_{0}\leftarrow z^{*}$

12: $z^{*}\leftarrow\min_{z}\eqref{eq:equation11A}$ s.t. (10b), (10c)

13: $\tilde{x},\tilde{u}\leftarrow\operatorname{interpolate}(z^{*})$

Algorithm 1 Algorithm for solving OCPs in form (1)

In the initialization phase, mesh variables $N$, $a$ and $b$ are set to small
values. The integrated residual minimization problem (9) is then solved using
warm starting and the solution checked if the obtained residual is below the
threshold tolerance $\epsilon_{tol}$. Note on line 2 of Algorithm 1 the
pseudo-code notation is simplistic, but the initial guess $z_{0}$ is not
directly set to $z^{*}$ since the size of $z$ increases as the mesh is
refined. However, $z_{0}$ will be an expanded and interpolated version of
$z^{*}$ as previously explained. Any suitable method can be used for
increasing $Q$ and $N$; we chose to double $N$ and $Q$ every time in our
example in Section IV. Once problem (9) has been solved and the residual
minimized, problem (10), which is a transcribed version of problem (1), is
solved using available NLP solvers.

##  IV Numerical Results

An optimal control problem solver based on the integrated residual method was
developed in the Julia v1.6 programming language. The package makes use of
barycentric interpolation routines as described in [1] to parametrize the
state and input variables. Numerical integration was performed using Gaussian
quadrature as detailed in [6]. Derivative information was obtained using
automatic differentiation (AD) tools [5] and supplied to the solver as the
gradient and Hessian of the Lagrangian function. The solver includes an
implementation of the flexible mesh scheme of Section III-A along with a fixed
mesh version, where time nodes $t_{i}$ are chosen to be in predefined
locations and not included as decision variables. In our implementation,
Chebyshev type 2 interpolation nodes and weights were used, since we want
internal supports $\tau_{i}^{0}$, $\tau_{i}^{a}$ to coincide with interval
boundaries $t_{i}$ and $t_{i+1}\,\forall i\in\\{0,\dots,N-1\\}$. The default
values for the number of intervals, state and input polynomial degree and
quadrature order are $N=5$, $a=2$, $b=1$ and $Q=3$.

In order to demonstrate the effectiveness of our method, we showcase the two
main features of our proposed method, namely superlinear convergence and the
ability to control the accuracy of state and input trajectories, on an optimal
control problem with a chattering solution. While solving optimal control
problems is the main focus of our work, the capabilities of the implemented
method can be used to solve feasibility problems and complex differential
equations as well (such as high index DAEs and differential inclusions).
Algorithm 1 was implemented using the interior point NLP solver Ipopt [14]
with a relative convergence tolerance set to $10^{-10}$. All tests were
performed on a laptop with an Intel® Core™ i7-4600U CPU at 2.10 GHz with 16 GB
of RAM.

###  IV-A Fuller problem description

We propose the numerical experiment to be an optimal control problem with a
discontinuous solution at non-trivial times in order to underline the
capability of our flexible-mesh optimal control solver to capture these
discontinuities. Additionally, we will analyse the impact of the desired
accuracy on the cost value for the numerically computed solution (which is a
lower bound for the exact optimal solution cost).

The chosen problem is a variation of the Fuller problem [4]

$\displaystyle\min_{p(\cdot),u(\cdot)}\quad$
$\displaystyle\int_{0}^{T}{p^{2}(t)}\mathop{}\\!\mathrm{d}t$ (11a) s.t.
$\displaystyle\ddot{p}(t)=u(t),\quad\quad\quad\,\,\,\,\forall t\in[0,T]$ (11b)
$\displaystyle u(t)\in[-0.01,0.01],\,\forall t\in[0,T]$ (11c) $\displaystyle
p(0)=p(T)=\dot{p}(T)=0,\,\dot{p}(0)=1,$ (11d)

with $T=300$ seconds, where $p(t)\in\mathbb{R}$ is the position and
$u(t)\in\mathbb{R}$ the control input at time $t$.

Since the chosen time $T$ was sufficiently large, the optimal control input
trajectory has a bang-bang structure with values alternating between $-0.01$
and $0.01$ and then reaching the steady state with $u(T)=0$. These switching
times are difficult to be captured by a numerical solver. Figure 1 displays
the state components $p$, $\dot{p}$ and the control input $u$ as obtained by
implementing Algorithm 1. Another relevant feature to observe is how the mesh
automatically becomes denser in the regions of sudden changes near the
switches and coarser where the solution is smoother.

Figure 1: Numerically computed control solution to Fuller problem using
flexible meshes with $N=20$ intervals, polynomial degrees $a=2$ , $b=1$,
flexibility parameter $\phi=0.5$ and desired accuracy
$\epsilon_{tol}=10^{-8}$. Blue dots indicate the location of mesh points.

Another aspect which motivates this choice of illustrative example is the
chattering phenomenon. As can be observed, instead of getting infinitely many
switches, we only capture a finite number of switches between values that are
not all on the input bounds. This behaviour is due to the specified
tolerances, as explained in Section IV-C below. The lower the tolerance, the
more accurate the numerical solution becomes.

###  IV-B Superlinear convergence for discontinuous solutions

Note that in this problem the solution can be represented by piecewise
polynomials of a sufficiently high degree. As a result, when setting $a\geq 2$
the minimum in (9) converges to $1.8729\cdot 10^{-15}$. In general, solutions
cannot be represented exactly using piecewise polynomials, hence we set the
polynomial degree to $a=b=1$ in order to reproduce the convergence behaviour
generally encountered for most practical problems.

Figure 2 presents the performance of our flexible mesh idea and compared to a
fixed mesh refinement procedure. The plot shows the variation of the minimum
integrated residual attainable for a certain number of mesh intervals $N$.
Note the scale is logarithmic and on the horizontal axis is plotted the
inverse of $N$. The slope of the blue line is approximately $2$ which is an
indicator of superlinear convergence. In contrast with our flexible mesh, the
red dots form a line of slope approximately equal to one for low values of $N$
and then start to plateau at a certain value.

Figure 2: Convergence of the Fuller problem as the mesh is refined with $N$
between 5 and 60 intervals, first order polynomial approximation $a=b=1$ and
flexibility parameter $\phi=0.5$.

###  IV-C Between accuracy and optimality

Figure 3 presents in green the solution for problem (9) using an increasing
number of mesh nodes $N$. Blue denotes the solution for (10) using an
increasing tolerance $\epsilon_{tol}$.

Figure 3: Impact of the desired accuracy on the cost value for the numerically
computed solution of the Fuller problem using flexible meshes with $\phi=0.5$.
Green dots denote solutions to problem (9) for different $N$ and when the cost
is not minimized. Blue dots represent the solutions to problem (10) for fixed
parameters $N=60$, $a=b=1$ as the tolerance $\epsilon_{tol}$ is varied.

The green dots represent what happens to the minimized integrated residual as
the number of intervals is changed. As can be observed, the cost does not
change significantly as the mesh is refined. However this is not always the
case and there can be situations where the solution of problem (9) is very far
away from the numerical solution computed from (10). Such an example is shown
in Figure 4. The blue dots are points on the Pareto front between numerical
lower bounds on the optimal cost value and discretization error for a fixed
number of intervals $N=20$.

Having an integrated residual based transcription allows the user to generate
an approximate numerical solution using a finite number of discretization
intervals $N$ and visualize the impact of different $\epsilon_{tol}$ tolerance
values on the computed solution. As one would expect, the objective value
increases as the tolerance is decreased, since the solution is captured more
accurately by our numerical scheme and a tighter lower bound to the exact
solution can be produced.

The Pareto front shown in Figure 4 uses a logarithmic scale for the y-axis and
a linear scale for the x-axis. In this example, decreasing the tolerance
$\epsilon_{tol}$ means that more full switches will be captured in the region
$t\approx 280$ seconds, thus improving the solution accuracy. However,
reducing $\epsilon_{tol}$ will increase the computational time. As can be
noticed, below a certain accuracy threshold the cost value no longer changes
significantly, meaning that the gap between the numerical and exact optimal
cost cannot be further reduced without increasing the number of mesh nodes
$N$. In general, we aim to find a value for $\epsilon_{tol}$ that can maintain
a good balance between solution accuracy and computational time.

Figure 4: Impact of the desired accuracy on the cost value for the numerically
computed solution for the Fuller problem using flexible meshes with
$\phi=0.5$. The green dot shows the solution to (9) for $N=20$, $a=2$, $b=1$.
Blue dots represent the solutions to problem (10) as the tolerance
$\epsilon_{tol}$ is varied.

##  V Conclusions and future works

Numerically solving optimal control problems for a given discretization mesh
involves a trade-off between computational time and solution accuracy. In our
example, we can use the change in cost value as a function of the residual
tolerance $\epsilon_{tol}$ as a metric to determine whether a satisfactory
solution has been found. While in most transcription methods this trade-off
cannot be easily ensured a priori, our proposed algorithm explicitly includes
the residual tolerance $\epsilon_{tol}$ as a parameter and is able to
construct the Pareto front between the cost value and residual tolerance
$\epsilon_{tol}$. This opens up the possibility for early termination in real-
time control applications. As discussed, it is important to use a
transcription method that relies on error measures that are integrated along
the entire solution trajectory to assess convergence, instead of measured
error at a finite number of sampled locations. In this way, the error between
mesh nodes can be accounted for.

When discontinuities are present in the solution, numerically approximating
state and input trajectories to an acceptable tolerance is especially
challenging. Our method proposes the use of a flexible mesh for capturing
discontinuities by including time mesh nodes in the decision vector. The
efficiency of this method was demonstrated through an illustrative example.

The implementation of the method is in an early development stage and many
improvements are possible in order to demonstrate its full capabilities.
Further research could be conducted on improving the mesh refinement process
by including early termination procedures, thus increasing the computational
efficiency of the overall solution process. Other future work could aim at
automatically increasing the quadrature order such that the constraints in (7)
can be removed. Theoretical convergence and performance guarantees also need
to be investigated.

## References

  * [1] Jean-Paul Berrut and Lloyd N Trefethen.  Barycentric lagrange interpolation.  SIAM review, 46(3):501–517, 2004.
  * [2] John T Betts.  Practical methods for optimal control and estimation using nonlinear programming.  SIAM, 2010.
  * [3] Ernest D. Eason.  A review of least-squares methods for solving partial differential equations.  International journal for numerical methods in engineering, 10(5):1021–1046, 1976.
  * [4] Anthony T. Fuller.  Study of an optimum non-linear control system.  International Journal of Electronics, 15(1):63–71, 1963.
  * [5] Michael Innes.  Don’t unroll adjoint: Differentiating ssa-form programs.  arXiv preprint arXiv:1810.07951, 2018.
  * [6] Steven G. Johnson.  QuadGK.jl: Gauss–Kronrod integration in Julia.  https://github.com/JuliaMath/QuadGK.jl, 2013.
  * [7] Fengjin Liu, William W. Hager, and Anil V. Rao.  Adaptive mesh refinement method for optimal control using nonsmoothness detection and mesh size reduction.  Journal of the Franklin Institute, 352(10):4081–4106, 2015.
  * [8] Daniele Mortari, Hunter Johnston, and Lidia Smith.  High accuracy least-squares solutions of nonlinear differential equations.  Journal of computational and applied mathematics, 352:293–307, 2019\.
  * [9] Yuanbo Nie and Eric C. Kerrigan.  Efficient and more accurate representation of solution trajectories in numerical optimal control.  IEEE Control Systems Letters, 4(1):61–66, 2019.
  * [10] Yuanbo Nie and Eric C. Kerrigan.  Solving dynamic optimization problems to a specified accuracy: An alternating approach using integrated residuals.  IEEE Transactions on Automatic Control, 2022.
  * [11] Lucian Nita, Eduardo M. G. Vila, Marta A. Zagorowska, Eric C. Kerrigan, Yuanbo Nie, Ian McInerney, and Paola Falugi.  Fast and accurate method for computing non-smooth solutions to constrained control problems.  In Proc. 20th European Control Conference (ECC), 2022.
  * [12] Luís Tiago Paiva and Fernando A.C.C. Fontes.  Adaptive time–mesh refinement in optimal control problems with state constraints.  Discrete & Continuous Dynamical Systems, 35(9):4553, 2015.
  * [13] James Blake Rawlings, David Q. Mayne, and Moritz Diehl.  Model predictive control: theory, computation, and design, volume 2.  Nob Hill Publishing Madison, 2017.
  * [14] Andreas Wächter and Lorenz T. Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Mathematical programming, 106(1):25–57, 2006.