# A geometric approach to approximating the limit set of eigenvalues for
banded Toeplitz matrices 111The work presented in this article is based on the
Master’s thesis of Teodor Bucht, available at http://lup.lub.lu.se/student-
papers/record/9109060.

Teodor Bucht222Centre for Mathematical Sciences, Lund University, Box 118,
22100 Lund, Sweden and Jacob S. Christiansen333Centre for Mathematical
Sciences, Lund University, Box 118, 22100 Lund, Sweden

###### Abstract

This article is about finding the limit set for banded Toeplitz matrices. Our
main result is a new approach to approximate the limit set $\Lambda(b)$ where
$b$ is the symbol of the banded Toeplitz matrix. The new approach is
geometrical and based on the formula
$\Lambda(b)=\cap_{\rho\in(0,\infty)}\operatorname{sp}T(b_{\rho})$. We show
that the full intersection can be approximated by the intersection for a
finite number of $\rho$’s, and that the intersection of polygon approximations
for $\operatorname{sp}T(b_{\rho})$ yields an approximating polygon for
$\Lambda(b)$ that converges to $\Lambda(b)$ in the Hausdorff metric. Further,
we show that one can slightly expand the polygon approximations for
$\operatorname{sp}T(b_{\rho})$ to ensure that they contain
$\operatorname{sp}T(b_{\rho})$. Then, taking the intersection yields an
approximating superset of $\Lambda(b)$ which converges to $\Lambda(b)$ in the
Hausdorff metric, and is guaranteed to contain $\Lambda(b)$. We implement the
algorithm in Python and test it. It performs on par to and better in some
cases than existing algorithms. We argue, but do not prove, that the average
time complexity of the algorithm is $O(n^{2}+mn\log m)$, where $n$ is the
number of $\rho$’s and $m$ is the number of vertices for the polygons
approximating $\operatorname{sp}T(b_{\rho})$. Further, we argue that the
distance from $\Lambda(b)$ to both the approximating polygon and the
approximating superset decreases as $O(1/\sqrt{k})$ for most of $\Lambda(b)$,
where $k$ is the number of elementary operations required by the algorithm.

Keywords. Banded Toeplitz matrices, limiting set of eigenvalues, polygon
approximations.

Mathematics Subject Classification. 15B05, 47B06, 47B35, 65F15.

##  1 Introduction

In this work we study the eigenvalues of $n\times n$ truncated infinite
Toeplitz matrices

$T_{n}(a):=\begin{pmatrix}a_{0}&a_{-1}&a_{-2}&\ldots&a_{-n+1}\\\
a_{1}&a_{0}&a_{-1}&\ldots&a_{-n+2}\\\ a_{2}&a_{1}&a_{0}&\ldots&\ldots\\\
\ldots&\ldots&\ldots&\ldots&a_{-1}\\\
a_{n-1}&a_{n-2}&\ldots&a_{1}&a_{0}\end{pmatrix},$ (1)

as $n\to\infty$. Here $a:\mathbb{T}\to\mathbb{C}$ is the associated symbol of
the matrix defined by

$a(t):=\sum_{n=-\infty}^{\infty}a_{n}t^{n},\quad t\in\mathbb{T}.$

We regard $T_{n}(a)$ as a truncation of the infinite Toeplitz matrix

$T(a):=(a_{j-k})_{j,k=0}^{\infty},$

which naturally can be regarded as an operator on $l^{2}(\mathbb{N})$.

A key idea for investigating $T_{n}(a)$ is that $T_{n}(a)$ should “mimic” the
properties of $T(a)$, at least for large $n$. In particular,
$\operatorname{sp}T_{n}(a)$ should mimic $\operatorname{sp}T(a)$, so that we
can gain insight about $\operatorname{sp}T_{n}(a)$ by computing
$\operatorname{sp}T(a)$. The exact result we will use that links the two
spectra will be stated below in (3).

Recall the following result due to Gohberg [10, 11] describing
$\operatorname{sp}T(a)$ for continuous symbols.

$\operatorname{sp}T(a)=a(\mathbb{T})\cup\\{\lambda\in\mathbb{C}\setminus
a(\mathbb{T}):\operatorname{wind}{(a-\lambda)}\neq 0\\},$

where $\operatorname{wind}(a-\lambda)$ denotes the winding number of $a$
around $\lambda$.

As for $T_{n}(a)$, denote its eigenvalues by
$\lambda_{j}^{(n)},\;j=1,2,\ldots,n$, and define the Borel measures

$\mu_{n}(E):=\frac{\\#\\{j:\lambda_{j}^{(n)}\in E\\}}{n},$

where $E\subset\mathbb{C}$ is any Borel subset. These measures describe which
fraction of the eigenvalues of $T_{n}(a)$ are located in $E$. When studying
eigenvalues of $T_{n}(a)$, a central question is whether there exists a
limiting measure $\mu$ that $\mu_{n}$ converge to, at least weakly, as
$n\to\infty$, and if it exists, we want to compute it. This limiting measure
is also referred to as the limiting eigenvalue distribution.

Another central problem when studying the eigenvalues of $T_{n}(a)$ is finding
the strong and weak limit sets defined by

$\begin{split}\Lambda_{s}&:=\bigl{\\{}\lambda\in\mathbb{C}:\lambda_{j_{n}}^{(n)}\to\lambda\text{
for some sequence
}j_{n}\to\infty\bigr{\\}}=\limsup\operatorname{sp}T_{n}(a),\\\
\Lambda_{w}&:=\bigl{\\{}\lambda\in\mathbb{C}:\lambda_{j_{k}}^{(n_{k})}\to\lambda\text{
for some sequences
}j_{k},n_{k}\to\infty\bigr{\\}}=\liminf\operatorname{sp}T_{n}(a).\end{split}$

Loosely speaking, $\Lambda_{s}$ and $\Lambda_{w}$ tell us where the
eigenvalues cluster. We will occasionally use the notation $\Lambda_{s}(a)$
and $\Lambda_{w}(a)$ to signify that we mean the limit sets with respect to
the symbol $a$.

When the symbol is a real valued $L^{\infty}$ function, the first Szegő limit
theorem (see, e.g., [6, Theorem 5.10]) describes the limit sets and the
limiting eigenvalue distribution. So the case of real-valued symbols is well
understood. However, a Toeplitz matrix is generally not hermitian, and so the
symbol is generally not real. Finding the limiting eigenvalue distribution and
limit sets for arbitrary complex-valued symbols remains an open problem.

The pioneers of dealing with such issues for complex-valued symbols were
Schmidt and Spitzer. In [16], they investigated the subclass of banded
Toeplitz matrices. A Toeplitz matrix is said to be banded it its symbol $b$ is
a Laurent polynomial, that is,

$b(t)=\sum_{n=-r}^{s}b_{n}t^{n},$

for some integers $r$ and $s$. This means that there only are a finite number
of non-zero elements in each row of $T_{n}(b)$.

A crucial insight of Schmidt and Spitzer is that if $\Lambda_{s}(b)$ and
$\Lambda_{w}(b)$ are to “mimic” the behavior of $\operatorname{sp}T(b)$, then
they are forced to also mimic the behavior of $\operatorname{sp}T(b_{\rho})$,
where $b_{\rho}(t):=b(\rho t)$. For $T_{n}(b_{\rho})$ and $T_{n}(b)$ only
differ by a similarity transform,

$T_{n}(b_{\rho})=\operatorname{diag}(\rho,\rho^{2},\ldots,\rho^{n})\;T_{n}(b)\;\operatorname{diag}(\rho^{-1},\rho^{-2},\ldots,\rho^{-n}).$

So $\Lambda_{s}(b_{\rho})=\Lambda_{s}(b)$ and
$\Lambda_{w}(b_{\rho})=\Lambda_{w}(b)$. Consequently, $\Lambda_{s}(b)$ and
$\Lambda_{w}(b)$ should also mimic

$\bigcap_{\rho\in(0,\infty)}\operatorname{sp}T(b_{\rho}).$ (2)

Schmidt and Spitzer also defined the polynomial

$\begin{multlined}Q(\lambda,z):=z^{r}\bigl{(}b(z)-\lambda\bigr{)}=\\\
b_{-r}+b_{-r+1}z+\ldots+(b_{0}-\lambda)z^{r}+\ldots+b_{s-1}z^{r+s-1}+b_{s}z^{r+s}\end{multlined}Q(\lambda,z):=z^{r}\bigl{(}b(z)-\lambda\bigr{)}=\\\
b_{-r}+b_{-r+1}z+\ldots+(b_{0}-\lambda)z^{r}+\ldots+b_{s-1}z^{r+s-1}+b_{s}z^{r+s}$

and showed that there is a connection between the roots of $Q(\lambda,z)$ and
the set in (2). In particular, if we for fixed $\lambda$ denote the roots of
$Q(\lambda,z)$ as $z_{j}(\lambda),\;j=1,2,\ldots,r+s$, and order them so that
$\left\lvert z_{j}(\lambda)\right\rvert\leq\left\lvert
z_{j+1}(\lambda)\right\rvert,\;j=1,2,\ldots,r+s-1$, Schmidt and Spitzer showed
that

$\Lambda(b):=\bigl{\\{}\lambda\in\mathbb{C}:\left\lvert
z_{r}(\lambda)\right\rvert=\left\lvert
z_{r+1}(\lambda)\right\rvert\bigr{\\}}=\bigcap_{\rho\in(0,\infty)}\operatorname{sp}T(b_{\rho}).$
(3)

The main result of Schmidt and Spitzer was proving that

$\Lambda_{s}(b)=\Lambda_{w}(b)=\Lambda(b),$

this also explains why $\Lambda(b)$ is simply called the limit set.
Additionally, they showed that $\Lambda(b)$ is a finite union of analytic
arcs. A few years later, Ullman [17] showed that $\Lambda(b)$ is a connected
set. Further significant results on the asymptotic eigenvalue behavior for
banded Toeplitz matrices include Hirschman [12], Widom [19, 20] and Duits and
Kuijlaars [9]. Hirschman derived the limiting eigenvalue distribution for
banded matrices, Widom simplified the proofs of both Schmidt and Spitzer and
Hirschman, and Duits and Kuijlaars characterized the limiting eigenvalue
distribution in terms of an equilibrium problem.

Also, noteworthy monographs on the subject of Toeplitz matrices include
Böttcher and Silbermann [6] and Böttcher and Grudsky [4]. In [6], the general
theory of Toeplitz matrices is presented, whilst [4] focuses on spectral
properties of banded Toeplitz matrices.

As explained, there is for banded Toeplitz matrices a well developed theory
that describes the structure of the limit set $\Lambda(b)$ and how dense the
eigenvalues cluster on $\Lambda(b)$. But how do we find $\Lambda(b)$ for
practical examples? A natural approach to finding $\Lambda(b)$ is to calculate
$\operatorname{sp}T_{n}(b)$ for $n$ of increasing size and simply observe
where the eigenvalues seem to cluster. However, in general this approach does
not work. Consider for instance $b(t)=t^{-4}+t^{1}$. Then one can show that

$\Lambda(b)=\bigl{\\{}\lambda\in\mathbb{C}:\lambda=re^{\frac{2\pi
i}{5}},\,0\leq r\leq 5\cdot 4^{-4/5}\bigr{\\}}.$

Additionally, one can show that $\operatorname{sp}T_{n}(b)\subset\Lambda(b)$
for all $n$. We compare $\Lambda(b)$ with the approximated eigenvalues of
$T_{200}(b)$ and $T_{400}(b)$ using NumPy. The results are shown in Figure 1,
where we see that the eigenvalues do not cluster around $\Lambda(b)$ and even
seem to stray away from $\Lambda(b)$ which we know that the true eigenvalues
do not. So what we see in Figure 1 is purely a result of numerical errors.

(a) $\Lambda(b)$ and the eigenvalues of $T_{200}(b)$ for $b(t)=t^{-4}+t$.

(b) $\Lambda(b)$ and the eigenvalues of $T_{400}(b)$ for $b(t)=t^{-4}+t$.

Figure 1: Examples of poor numerical behavior for computed eigenvalues.

In the recent article [3], it is stated that “Finding this limiting set
[refering to $\Lambda(b)$] nevertheless remains a challenge”. Then the authors
propose an algorithm that is based on [2], which essentially looks at the
roots of $Q(z,\lambda)$. This algorithm produces some points contained in the
limit set, and the more sample points used, the better the representation will
be. However, this algorithm requires finding zeros of polynomials, which is a
hard problem if the degree of the polynomial is big, and for some cases a
large number of sampling points are needed. Instead of searching for
$\Lambda(b)$ algebraically, we present a novel approach to finding it
geometrically.

###  1.1 Previous and new work on computing $\Lambda(b)$

A naive approach to calculating $\Lambda(b)$ is to sample a grid in
$\mathbb{C}$ and check what points almost are in the set, by solving
$b(z)=\lambda$ and looking at the root sizes. However, this approach is time
consuming and it is not obvious how a point “almost” being in $\Lambda(b)$
should be interpreted.

In [3], a better algorithm for calculating $\Lambda(b)$ is presented. The idea
is: instead of sampling $\lambda$, sample $\varphi\in(0,2\pi)$ and solve

$b(z)-b(ze^{i\varphi})=0.$ (4)

This approach works since if $\lambda\in\Lambda(b)$ then two of the roots of
$b(z)-\lambda=0$, $z_{r}(\lambda)$ and $z_{r+1}(\lambda)$ have the same
modulus, so $z_{r+1}(\lambda)/z_{r}(\lambda)=e^{i\varphi}$, which implies that
$z_{r}(\lambda)$ is a root of (4). Note that in the case
$z_{r}(\lambda)=z_{r+1}(\lambda)$, (4) is the equation $0=0$. But for double
roots, the derivative also has a zero, so to find the cases where
$z_{r}(\lambda)=z_{r+1}(\lambda)$, we solve

$b^{\prime}(z)=0.$ (5)

From solving (5), and (4) for the sampled $\varphi$’s, we get some candidates
$z_{k}$. For each of these candidates, we can calculate the corresponding
$\lambda_{k}:=b(z_{k})$ and check each of these $\lambda_{k}$ for membership
in $\Lambda(b)$ by solving the equation

$b(z)=\lambda_{k}$

and checking if the $r$’th and $(r+1)$’th smallest roots (when ordered by
modulus) have the same absolute value.

The output from this algorithm is a set of points belonging to $\Lambda(b)$,
and it involves numerically finding the roots of multiple polynomials that
have degree $r+s$. For fixed $r$ and $s$, the computational time is at most
$O(N)$, where $N$ is the amount of sampled $\varphi$’s.

######  Example 1.1.

Consider the symbol

$b(t)=-2t^{-1}+4(1-i)+7it-3(1+i)t^{2}+t^{3}.$ (6)

This symbol is presented in [3] as an example of where the numerically
computed eigenvalues for $T_{n}(b_{\rho})$ do not converge to the true
$\Lambda(b)$, and where the computed eigenvalues vary with $\rho$, which they
theoretically do not. We estimated $\Lambda(b)$ using the algebraic approach
just presented. The results can be seen in Figure 2. If we zoom in on the top
part of the limit set, as we do in Figure 2(b), it is not clear whether
$\Lambda(b)$ forms a “cross” and separates the complex plane or not. It took
roughly 3 minutes to run the algorithm.

(a) $\Lambda(b)$ estimated using the algebraic approach (from section 1.1)
with $2\cdot 10^{5}$ $\varphi$’s sampled uniformly in $[0,\pi]$.

(b) A zoomed in version of 2(a).

Figure 2: The results from estimating $\Lambda(b)$ using the algebraic
algorithm presented above with $b$ as in (6).

The idea behind our new approach presented in this article is to approximate
$\Lambda(b)$ using (2). To do this, we first show that only a compact set of
$\rho$’s is needed in the intersection, and then we introduce polygons
approximating $\operatorname{sp}T(b_{\rho})$, which we denote by
$\operatorname{sp}T(b_{\rho}^{D})$. We show that the intersection of
$\operatorname{sp}T(b_{\rho_{i}}^{D})$ for a finite number of sampled
$\rho_{i}$’s converges to $\Lambda(b)$ in the Hausdorff metric as the number
of sampled $\rho_{i}$’s tend to infinity. This gives us a way to approximate
$\Lambda(b)$; we first find the compact interval of $\rho$’s needed in (2),
then sample $\rho$’s in this interval, construct polygon approximations for
$\operatorname{sp}T(b_{\rho_{i}})$, and finally compute the intersection which
yields an approximating polygon for $\Lambda(b)$.

Further, we show that one can slightly expand the polygon approximations
$\operatorname{sp}T(b_{\rho_{i}}^{D})$ to make sure that the expanded polygons
contain $\operatorname{sp}T(b_{\rho})$. Computing the intersection of the
expanded polygons then yields an approximating superset of $\Lambda(b)$ that
(contains $\Lambda(b)$ and) approaches $\Lambda(b)$ in the Hausdorff metric.
Additionally, we argue but do not prove that the Hausdorff distance between
$\Lambda(b)$ and the approximating superset should decrease as
$O(1/\sqrt{k})$, where $k$ is the number of basic operations needed to run the
algorithm.

###  1.2 Outline

This article is organized as follows. In Section 2 we present the theoretic
background of our algorithm, and show that the approximating polygon and
approximating superset being produced both converge to $\Lambda(b)$. In
Section 3 we discuss how the algorithm can be implemented. We also suggest a
better way of sampling $\rho$’s for the intersection, and compare the naive
and improved way of sampling $\rho$’s. In Section 4 we reason about the time
complexity and convergence speed of the algorithm. Lastly, we provide some
concluding remarks in Section 5.

##  2 A geometric approach

One way to look at the problem of finding $\Lambda(b)$ is the result of
Schmidt and Spitzer [16] presented in (3), which says

$\bigcap_{\rho\in(0,\infty)}\operatorname{sp}T(b_{\rho})=\Lambda(b).$ (7)

We can either work with the left hand side and try to compute it
geometrically, or look at the right hand side and compute it algebraically,
which is what previous works have done, see subsection 1.1. We will now show a
geometric approach to approximating $\Lambda(b)$, which involves approximating
the left-hand side in (7).

The first observation that speaks in favor of
$\cap_{\rho\in(0,\infty)}\operatorname{sp}T(b_{\rho})$ even being possible to
compute numerically is that not all $\rho$’s are needed.

######  Proposition 2.1.

Given a Laurent polynomial $b$, there exists $\rho_{l},\rho_{h}$ with
$0<\rho_{l}\leq\rho_{h}<\infty$ such that

$\Lambda(b)=\bigcap_{\rho\in(0,\infty)}\operatorname{sp}T(b_{\rho})=\bigcap_{\rho\in[\rho_{l},\rho_{h}]}\operatorname{sp}T(b_{\rho}).$
(8)

Furthermore, examples of such $\rho_{l}$ and $\rho_{h}$ can be found by
solving a polynomial equation with real coefficients.

###### Proof.

The idea behind the proof is that for small enough $\rho$ the lowest order
term in $b$ will dominate, so then $b_{\rho}(\mathbb{T})$ will roughly trace
out $r$ circles of radius $b_{-r}\rho^{-r}$ winding in negative direction.
Thus, intersecting with $\operatorname{sp}T(b_{\rho})$ for small enough $\rho$
does not erase any new points. The same argument holds for large enough
$\rho$.

Now, choose $K:=\left\lVert b_{1}\right\rVert_{\infty}$ and $0<\rho_{l}\leq
1\leq\rho_{h}$ such that $\operatorname{sp}T(b_{\rho})\supset B(0,K)$ for
$\rho<\rho_{l}$ and $\rho_{h}<\rho$, where $B(0,K)$ is the closed ball
centered at $0$ with radius $K$. Then we have

$\bigcap_{\rho\in(0,\infty)}\operatorname{sp}T(b_{\rho})\supset
B(0,K)\cap\Biggl{(}\,\bigcap_{\rho\in[\rho_{l},\rho_{h}]}\operatorname{sp}T(b_{\rho})\Biggr{)}=\bigcap_{\rho\in[\rho_{l},\rho_{h}]}\operatorname{sp}T(b_{\rho})\supset\Lambda(b),$

where the equality holds since $\operatorname{sp}T(b_{1})\subset B(0,K)$ and
$\rho_{l}\leq\rho_{h}$. As the left-most and right-most sets are equal in the
chain of inclusions, all sets must be equal, showing our desired result.

Now we motivate how $\rho_{l}$ and $\rho_{h}$ with the correct properties can
be found. As before, $K=\left\lVert b_{1}\right\rVert_{\infty}$. Let
$\rho_{l}$ be the minimum of 1 and the smallest real positive solution to
$\left\lvert b_{-r}\right\rvert\rho^{-r}-\sum_{n=-r+1}^{s}\left\lvert
b_{n}\right\rvert\rho^{n}=K$. Hence $\left\lvert
b_{-r}\right\rvert\rho^{-r}-\sum_{n=-r+1}^{s}\left\lvert
b_{n}\right\rvert\rho^{n}>K$ for all $\rho<\rho_{l}$. We now estimate the size
of $b_{\rho}(t)$. For $\rho$ in $0<\rho<\rho_{l}$, we have

$\begin{split}\left\lvert
b_{\rho}(t)\right\rvert=\left\lvert\sum_{n=-r}^{s}b_{n}\rho^{n}t^{n}\right\rvert&\geq\left\lvert
b_{-r}\right\rvert\rho^{-r}-\left\lvert\sum_{n=-r+1}^{s}b_{n}\rho^{n}t^{n}\right\rvert\\\
&\geq\left\lvert b_{-r}\right\rvert\rho^{-r}-\sum_{n=-r+1}^{s}\left\lvert
b_{n}\right\rvert\rho^{n}>K.\end{split}$

Hence $\operatorname{sp}T(b_{\rho})\supset B(0,K)$ for $\rho$ in
$0<\rho<\rho_{l}$. Analogously we can define $\rho_{h}$ as the maximum of 1
and the largest positive root to $\left\lvert
b_{s}\right\rvert\rho^{s}-\sum_{n=-r}^{s-1}\left\lvert
b_{n}\right\rvert\rho^{n}=K$ and do similar calculations to get that
$\operatorname{sp}T(b_{\rho})\supset B(0,K)$ for $\rho>\rho_{h}$, which
concludes our proof. ∎

So, only a compact set of $\rho$’s is needed to obtain $\Lambda(b)$. A natural
approach to approximating $\Lambda(b)$ is to select a partition
$\rho_{l}=\rho_{0}<\rho_{1}<\ldots<\rho_{n-1}<\rho_{n}=\rho_{h}$ and compute

$\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}).$ (9)

The granularity of the partition $(\rho_{j})_{j=0}^{n}$ is defined as

$\Delta:=\max_{0\leq j\leq n-1}(\rho_{j+1}-\rho_{j}).$

The question is now if (9) will approach $\Lambda(b)$, and in what way? In
what follows, we will show that (9) approaches $\Lambda(b)$ in the Hausdorff
metric as the granularity of $(\rho_{j})_{j=0}^{n}$ tends to $0$. Recall that
the Hausdorff metric, denoted $d_{H}$, operates on the set $\mathcal{C}$
consisting of all the bounded non-empty closed subsets of $\mathbb{C}$. It is
defined by

$d_{H}(A,B):=\inf\left\\{\epsilon\,:\,A\subset(B)_{\epsilon},B\subset(A)_{\epsilon}\right\\},$

where $(X)_{\epsilon}$ denotes the $\epsilon$–fattening of $X$, that is,

$(X)_{\epsilon}:=\bigcup_{x\in X}\left\\{z\in\mathbb{C}\,:\,\left\lvert
z-x\right\rvert\leq\epsilon\right\\}.$

As is well known, $(\mathcal{C},d_{H})$ forms a complete metric space, see
e.g. [7, p. 121–124].

To show that (9) approaches $\Lambda(b)$, we need the following lemma.

######  Lemma 2.2.

For a Laurent polynomial $b$, the map

$\rho\mapsto\operatorname{sp}T(b_{\rho})$

from $(\mathbb{R}^{+},|\cdot|)\to(\mathcal{C},d_{H})$ is continuous.

###### Proof.

Let $(\mathcal{T}_{C},\left\lVert\,\cdot\,\right\rVert)$ denote the space of
Toeplitz operators with continuous symbol endowed with the operator norm. We
see that our map is the composition of three maps,

$\begin{split}\rho\mapsto
b_{\rho},&\quad(\mathbb{R}^{+},|\cdot|)\to(C(\mathbb{T}),||\cdot||_{\infty}),\\\
a\mapsto
T(a),&\quad(C(\mathbb{T}),||\cdot||_{\infty})\to(\mathcal{T}_{C},||\cdot||),\\\
T(a)\mapsto\operatorname{sp}T(a),&\quad(\mathcal{T}_{C},||\cdot||)\to(\mathcal{C},d_{H}).\end{split}$
(10)

The first and second maps in (10) are clearly continuous, and by [13, Theorem
10] the third map is also continuous. ∎

######  Remark 2.3.

Recall that the map $A\mapsto\operatorname{sp}A$, from the set of bounded
operators on a Hilbert space to $(\mathcal{C},d_{H})$ is not continuous is
general, see e.g. [15, Example 3.10]. According to [13], it was conjectured
however that $A\mapsto\operatorname{sp}A$ is continuous on the manifold of
Toeplitz operators. But this conjecture was proven false in [5], where the
authors constructed symbols $a_{n},a\in L^{\infty}(\mathbb{T})$ such that
$\left\lVert a_{n}-a\right\rVert_{\infty}\to 0$ but
$d_{H}(\operatorname{sp}T(a_{n}),\operatorname{sp}T(a))\not\rightarrow 0$. It
should be noted that continuity is proven for a much larger class of symbols
than $C(\mathbb{T})$ in [13].

Now we can use Lemma 2.2 to show the promised result, that (9) approaches
$\Lambda(b)$ in the limit.

######  Theorem 2.4.

Let $\rho_{l}$ and $\rho_{h}$ be bounds as in Theorem 2.1 and let
$(\rho_{j})_{j=0}^{n}$ be a partition of $[\rho_{l},\rho_{h}]$ of granularity
$\Delta$. Then we have

$\lim_{\Delta\to
0}\Biggl{[}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}})\Biggr{]}=\Lambda(b)$

in the Hausdorff metric.

###### Proof.

From Lemma 2.2 we have that the map $\rho\mapsto\operatorname{sp}T(b_{\rho})$
from $[\rho_{l},\rho_{h}]\to\mathcal{C}$ is uniformly continuous. So for a
given $\epsilon>0$ there is a $\delta>0$ such that for all
$\rho,\rho^{*}\in[\rho_{l},\rho_{h}]$ with
$\left\lvert\rho-\rho^{*}\right\rvert<\delta$ we have
$d_{H}(\operatorname{sp}T(b_{\rho^{*}}),\operatorname{sp}T(b_{\rho}))<\epsilon$.
It thus suffices to show that for all partitions $(\rho_{j})_{j=0}^{n}$ with
granularity $\Delta<2\delta$ we have

$d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\bigcap_{\rho\in[\rho_{l},\rho_{h}]}\operatorname{sp}T(b_{\rho})\Biggr{)}\leq\epsilon.$

Naturally we always have
$\bigcap_{\rho\in[\rho_{l},\rho_{h}]}\operatorname{sp}T(b_{\rho})\subset\left(\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}})\right)_{\epsilon}$.
So we need only to prove the opposite inclusion. To do this, take an arbitrary

$\lambda\notin\Biggl{(}\,\bigcap_{\rho\in[\rho_{l},\rho_{h}]}\operatorname{sp}T(b_{\rho})\Biggr{)}_{\epsilon}=\bigcap_{\rho\in[\rho_{l},\rho_{h}]}\bigl{(}\operatorname{sp}T(b_{\rho})\bigr{)}_{\epsilon}.$

So there is a $\rho^{*}\in[\rho_{l},\rho_{h}]$ such that
$\lambda\notin(\operatorname{sp}T(b_{\rho^{*}}))_{\epsilon}$. Since
$\Delta<2\delta$ there is a $\rho_{k}$ in the partition with
$\left\lvert\rho^{*}-\rho_{k}\right\rvert<\delta$. This implies
$d_{H}(\operatorname{sp}T(b_{\rho^{*}}),\operatorname{sp}T(b_{\rho}))<\epsilon$,
which means
$(\operatorname{sp}T(b_{\rho^{*}}))_{\epsilon}\supset\operatorname{sp}T(b_{\rho_{k}})$.
So $\lambda\notin\operatorname{sp}T(b_{\rho_{k}})$, and we are done. ∎

Theorem 2.4 hints at great possibilities for computing $\Lambda(b)$
numerically. But it is difficult to represent $\operatorname{sp}T(b_{\rho})$
numerically. A possible solution is to approximate
$\operatorname{sp}T(b_{\rho})$ by a polygon with vertices sampled along
$b_{\rho}(\mathbb{T})$. In what follows we shall show that such a polygon
approximation indeed gives the desired result. We will partition $[0,2\pi]$
with $0=v_{0}<v_{1}<\ldots<v_{n-1}<v_{n}=2\pi$. The granularity for
$(v_{j})_{j=0}^{n}$ is defined as before. A partition of $[0,2\pi]$ defines a
polygon discretization $b_{\rho}^{D}$ approximating $b_{\rho}$, on the form
$b_{\rho}^{D}:\mathbb{T}\to\mathbb{C}$,

$b_{\rho}^{D}(e^{iv})=\frac{v-v_{j}}{v_{j+1}-v_{j}}\left(b_{\rho}(e^{iv_{j+1}})-b_{\rho}(e^{iv_{j}})\right)+b_{\rho}(e^{iv_{j}}),\;\;v\in[v_{j},v_{j+1}).$
(11)

Since $b_{\rho}^{D}$ is a continuous symbol, we have

$\operatorname{sp}T(b_{\rho}^{D})=b_{\rho}^{D}(\mathbb{T})\cup\left\\{\lambda\in\mathbb{C}:\operatorname{wind}(b_{\rho}^{D}-\lambda)\neq
0\right\\}.$

Naturally, we wish to bound the Hausdorff distance between the spectrum of
$T(b_{\rho})$ and $T(b_{\rho}^{D})$.

######  Proposition 2.5.

Let $b$ be a Laurent polynomial and $b_{\rho}^{D}$ a polygon discretization of
$b_{\rho}$ as is (11). Then

$d_{H}\bigl{(}\operatorname{sp}T(b_{\rho}),\operatorname{sp}T(b_{\rho}^{D})\bigr{)}\leq\left\lVert
b_{\rho}-b_{\rho}^{D}\right\rVert_{\infty}.$

###### Proof.

Take an arbitrary $\delta_{\rho}>\left\lVert
b_{\rho}-b_{\rho}^{D}\right\rVert_{\infty}$. Our goal is to show
$(\operatorname{sp}T(b_{\rho}))_{\delta_{\rho}}\supset\operatorname{sp}T(b_{\rho}^{D})$
and
$(\operatorname{sp}T(b_{\rho}^{D}))_{\delta_{\rho}}\supset\operatorname{sp}T(b_{\rho})$,
which we can do through a geometric argument. We will proceed by taking
$\lambda$ either $\notin(\operatorname{sp}T(b_{\rho}))_{\delta_{\rho}}$ or
$\notin(\operatorname{sp}T(b_{\rho}^{D}))_{\delta_{\rho}}$. Since
$\mathcal{R}(a)\subset\operatorname{sp}T(a)$ for any continuous symbol $a$, we
have by choice of $\delta_{\rho}$ that $\inf_{t\in\mathbb{T}}\left\lvert
b_{\rho}(t)-\lambda\right\rvert>0$ and $\inf_{t\in\mathbb{T}}\left\lvert
b_{\rho}^{D}(t)-\lambda\right\rvert>0$. Due to this we can parameterize

$b_{\rho}(e^{iv})-\lambda=r(v)e^{i\theta(v)}\;\text{ and
}\;b_{\rho}^{D}(v)-\lambda=r_{D}(v)e^{i\theta_{D}(v)},$

with continuous $r,r_{D},\theta,\theta_{D}$. Recall that
$b_{\rho}^{D}(e^{iv_{j}})=b_{\rho}(e^{iv_{j}})$ where $(v_{j})_{j=0}^{n}$ is
the partition of $[0,2\pi]$ defining $b_{\rho}^{D}$. The winding numbers can
thus be written as

$\begin{split}\operatorname{wind}(b_{\rho}^{D}-\lambda)&=\frac{1}{2\pi}\bigl{(}\theta_{D}(v_{n})-\theta_{D}(v_{0})\bigr{)}=\frac{1}{2\pi}\sum_{j=1}^{n}\theta_{D}(v_{j})-\theta_{D}(v_{j-1}),\\\
\operatorname{wind}(b_{\rho}-\lambda)&=\frac{1}{2\pi}\bigl{(}\theta(v_{n})-\theta(v_{0})\bigr{)}=\frac{1}{2\pi}\sum_{j=1}^{n}\theta(v_{j})-\theta(v_{j-1}).\end{split}$
(12)

Our strategy is to show that
$\theta_{D}(v_{j})-\theta_{D}(v_{j-1})=\theta(v_{j})-\theta(v_{j-1})$ for all
$j$. Begin with $j=1$, and assume that $\theta(0)=\theta_{D}(0)$. Denote
$b_{\rho}^{D}(v_{0})$ by $A$ and $b_{\rho}^{D}(v_{1})$ by $B$ and the line
through $A$ and $B$ by $l$. There are two cases to consider: the first case is
when $A\neq B$ and the projection $P$ of $\lambda$ onto $l$ lies on the
segment $AB$, the second case is when $A=B$ or the projection $P$ of $\lambda$
onto $l$ does not lie on the segment $AB$.

$B$$A$$\lambda$$P$$r_{\lambda}$$C$$b_{\rho}([v_{0},v_{1}])$ (a) Illustration
of the first case.

$B$$A$$\lambda$$C$$r_{\lambda}$$b_{\rho}([v_{0},v_{1}])$ (b) Illustration of
the second case.

Figure 3: Illustrations of the construction in the proof of Proposition 2.5.

In the first case consider the line through $P$ and $\lambda$. Note that
$P\lambda\perp AB$. Cut the line through $P$ and $\lambda$ at $\lambda$ so we
end up with a ray from $\lambda$ not intersecting $AB$. Denote this ray by
$r_{\lambda}$. Assume for contradiction that $r_{\lambda}$ intersects
$b_{\rho}(e^{iv})$ for some $v\in[v_{0},v_{1}]$, say at the point $C$, see
Figure 3(a) for an illustration. As $\left\lvert CP\right\rvert\leq\left\lVert
b_{\rho}-b_{\rho}^{D}\right\rVert_{\infty}<\delta_{\rho}$, this is a
contradiction since $\left\lvert P\lambda\right\rvert\leq\left\lvert
CP\right\rvert$ and $\left\lvert\lambda C\right\rvert\leq\left\lvert
CP\right\rvert$ (since $\lambda$ lies on the segment $CP$). Hence neither
$b_{\rho}(e^{iv})$ nor $b_{\rho}^{D}(e^{iv})$ intersect $r_{\lambda}$ for
$v\in[v_{0},v_{1}]$, so all angles $\theta(v)$ and $\theta_{D}(v)$ stay in the
same cut plane. Hence $\theta(v_{1})=\theta_{D}(v_{1})$ since
$\theta(v_{0})=\theta_{D}(v_{0})$.

We can handle the second case in a similar manner. Assume without loss of
generality that $A$ is closer to $\lambda$ than $B$, and let $r_{\lambda}$ be
the ray originating in $\lambda$ that lies on the line through $A$ and
$\lambda$ and does not contain $A$. Assume as before that $r_{\lambda}$
intersects $b_{\rho}(e^{iv})$ at some point $C$, see Figure 3(b) for an
illustration. We see that $\left\lvert AC\right\rvert<\delta_{\rho}$ which
using a similar argument again gives that neither $b_{\rho}(e^{iv})$ nor
$b_{\rho}^{D}(e^{iv})$ can intersect $r_{\lambda}$ for $v\in[v_{0},v_{1}]$.
This implies that $\theta(v_{1})=\theta_{D}(v_{1})$. Note that the case $A=B$
is also handled by this construction.

All in all we have shown that
$\theta_{D}(v_{j})-\theta_{D}(v_{j-1})=\theta(v_{j})-\theta(v_{j-1})$ for
$j=1$ and the argument can be continued inductively which shows it for all
$1\leq j\leq n$. From (12) we now get that
$\operatorname{wind}(b_{\rho}-\lambda)=\operatorname{wind}(b_{\rho}^{D}-\lambda)$
which means that
$\lambda\notin(\operatorname{sp}T(b_{\rho}))_{\delta_{\rho}}\Rightarrow\lambda\notin\operatorname{sp}T(b_{\rho}^{D})$
and
$\lambda\notin(\operatorname{sp}T(b_{\rho}^{D}))_{\delta_{\rho}}\Rightarrow\lambda\notin\operatorname{sp}T(b_{\rho})$,
which is what we wanted to show. ∎

Proposition 2.5 gives us a good bound on the error introduced by the polygon
approximation. Next we bound $\left\lVert
b_{\rho}-b_{\rho}^{D}\right\rVert_{\infty}$ and use this to get a uniform
bound for the error.

######  Lemma 2.6.

Let $\rho_{l}$ and $\rho_{h}$ be bounds as in Proposition 2.1 and let
$(v_{j})_{j=0}^{n}$ be a partition of $[0,2\pi]$ of granularity $\Delta$. Then
there exists $C>0$ not depending on $\Delta$ such that

$d_{H}\left(\operatorname{sp}T(b_{\rho}),\operatorname{sp}T(b_{\rho}^{D})\right)\leq
C\Delta^{2}$

for all $\rho\in[\rho_{l},\rho_{h}]$.

###### Proof.

When comparing $b_{\rho}^{D}$ with $b_{\rho}$, it is convenient to write

$b(\rho e^{iv})=x_{\rho}(v)+iy_{\rho}(v).$

Further, we see that $x_{\rho}(v)$ and $y_{\rho}(v)$ are smooth in both $\rho$
and $v$. In particular this implies that $x_{\rho}^{\prime\prime}(v)$ and
$y_{\rho}^{\prime\prime}(v)$ are continuous on
$[\rho_{l},\rho_{h}]\times[0,2\pi]$, so

$C:=\sup_{\rho\in[\rho_{l},\rho_{h}],v\in[0,2\pi]}\left\lvert
x_{\rho}^{\prime\prime}(v)\right\rvert+\left\lvert
y_{\rho}^{\prime\prime}(v)\right\rvert<\infty.$

By the mean value theorem for $x_{\rho}$ and $y_{\rho}$ separately we have for
$v\in(v_{j},v_{j+1})$

$\frac{b(\rho e^{iv_{j+1}})-b(\rho
e^{iv_{j}})}{v_{j+1}-v_{j}}=x_{\rho}^{\prime}(v_{jx})+i{y}_{\rho}^{\prime}(v_{jy}),$

for some $v_{jx},v_{jy}\in(v_{j},v_{j+1})$. Using (11) and the mean value
theorem repeatedly, we get for $v\in(v_{j},v_{j+1})$

$\begin{split}\left\lvert b_{\rho}^{D}(e^{iv})-b(\rho
e^{iv})\right\rvert&=\left\lvert(v-v_{j})({x}_{\rho}^{\prime}(v_{jx})+{y}_{\rho}^{\prime}(v_{jy}))+b(\rho
e^{iv_{j}})-b(\rho e^{iv})\right\rvert\\\
&=\left\lvert(v-v_{j})({x}_{\rho}^{\prime}(v_{jx})+{y}_{\rho}^{\prime}(v_{jy}))-(v-v_{j})({x}_{\rho}^{\prime}(v_{jx}^{*})+{y}_{\rho}^{\prime}(v_{jy}^{*}))\right\rvert\\\
&\leq\left\lvert
v-v_{j}\right\rvert\left(\left\lvert{x}_{\rho}^{\prime}(v_{jx})-{x}_{\rho}^{\prime}(v_{jx}^{*})\right\rvert+\left\lvert{y}_{\rho}^{\prime}(v_{jy})-{y}_{\rho}^{\prime}(v_{jy}^{*})\right\rvert\right),\\\
&\leq\Delta^{2}\bigl{(}\left\lvert{x}_{\rho}^{\prime\prime}(v_{jx}^{**})\right\rvert+\left\lvert{y}_{\rho}^{\prime\prime}(v_{jy}^{**})\right\rvert\bigr{)}\\\
&\leq C\Delta^{2},\end{split}$ (13)

where $v_{jx}^{*},v_{jy}^{*}\in(v_{j},v)$ and $v_{jx}^{**},v_{jy}^{**}$
between $v_{jx}$ and $v_{jx}^{*}$, and $v_{jy}$ and $v_{jy}^{*}$ respectively.
Now we have $\left\lVert b_{\rho}-b_{\rho}^{D}\right\rVert_{\infty}\leq
C\Delta^{2}$, which using Proposition 2.5 gives the desired result. ∎

With this result we are ready to prove a stronger version of Theorem 2.4.

######  Theorem 2.7.

Let $\rho_{l}$ and $\rho_{h}$ be bounds as in Theorem 2.1, and let
$(\rho_{j})_{j=0}^{n}$ denote partitions of $[\rho_{l},\rho_{h}]$ with
granularity $\Delta_{\rho}$. Moreover, let $(v_{j})_{j=0}^{m}$ denote
partitions of $[0,2\pi]$ with granularity $\Delta_{v}$. Then

$\lim_{\Delta_{\rho},\Delta_{v}\to
0}\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D})=\Lambda(b)$ (14)

in the Hausdorff metric.

###### Proof.

Fix an $\epsilon>0$. We can use the triangle inequality for $d_{H}$ to
estimate

$\begin{multlined}d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\Lambda(b)\Biggr{)}\leq
d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}})\Biggr{)}\\\
+d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\Lambda(b)\Biggr{)}.\end{multlined}d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\Lambda(b)\Biggr{)}\leq
d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}})\Biggr{)}\\\
+d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\Lambda(b)\Biggr{)}.$
(15)

Let $\delta_{\rho}$ be such that
$d_{H}\left(\cap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\Lambda(b)\right)<\epsilon/2$
for all partitions of $[\rho_{l},\rho_{h}]$ with granularity less than
$\delta_{\rho}$, which exists by Theorem 2.4. Let $C>0$ be as in Lemma 2.6 and
choose $\delta_{v}$ small enough so that
$d_{H}\left(\operatorname{sp}T(b_{\rho}),\operatorname{sp}T(b_{\rho}^{D})\right)\leq
C\delta_{v}^{2}<\epsilon/2$ for all partitions of $[0,2\pi]$ with granularity
less than $\delta_{v}$ and all $\rho\in[\rho_{l},\rho_{h}]$. Using this we can
prove that
$d_{H}\bigl{(}\cap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\cap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}})\bigr{)}\leq\epsilon/2$,
by noting that

$\begin{split}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D})\Biggr{)}_{\epsilon/2}=\bigcap_{j=0}^{n}\left(\operatorname{sp}T(b_{\rho_{j}}^{D})\right)_{\epsilon/2}&\supset\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\\\
\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}})\Biggr{)}_{\epsilon/2}=\bigcap_{j=0}^{n}\left(\operatorname{sp}T(b_{\rho_{j}})\right)_{\epsilon/2}&\supset\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}).\end{split}$

This concludes our proof since we now have
$d_{H}\bigl{(}\cap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\Lambda(b)\bigr{)}\leq\epsilon/2+\epsilon/2=\epsilon$.
∎

Thanks to Theorem 2.7 we have the outline of an algorithm that we know will
converge to $\Lambda(b)$. Given $b$ we sample $\rho$’s in
$[\rho_{l},\rho_{h}]$, and for each sampled $\rho$ construct an approximating
polygon according to (11). Then, taking the intersection with respect to non-
zero winding number of all the sampled polygons yields an approximating
polygon for $\Lambda(b)$. By Theorem 2.7 we know that the intersection of the
polygons will approach $\Lambda(b)$ as we decrease the granularities.
Pseudocode for the algorithm can be seen in Algorithm 1 below.

Algorithm 1 Basic approach to calculating $\Lambda(b)$ geometrically

procedure CalcLimitSet($b,n,m$)

$b$: the symbol,

$n$: the number of sampled $\rho$’s,

$m$: number of sampled $v$’s.

$\rho_{l},\rho_{h}\leftarrow$ bounds from Proposition 2.1

rhos $\leftarrow$ sample $n+1$ $\rho$’s in $[\rho_{l},\rho_{h}]$

vs $\leftarrow$ sample $m+1$ $v$’s in $[0,2\pi]$

$\Lambda\leftarrow$ $b_{\text{rhos[1]}}^{D}(vs)$

for $i\leftarrow 2$ to $n+1$ do

$\Lambda\leftarrow\Lambda\cap b_{\text{rhos[i]}}^{D}(vs)$

end for

return $\Lambda$

end procedure

######  Remark 2.8.

It would be useful to quantify the error introduced by the polygon
approximation. In fact, we see from the proof of Theorem 2.7 that there exists
$C>0$ such that

$d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\Lambda(b)\Biggr{)}\leq
C\Delta_{v}^{2}+d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\Lambda(b)\Biggr{)}.$
(16)

In other words, the polygon approximation gives a term that converges very
fast to 0, so the main error to worry about is the second term on the right
hand side of (16).

Now we have a way of producing an approximating polygon, but from just looking
at the approximating polygon we do not get a lot of information about the size
of the error, but this is likely something that is very useful to know in
practical applications. For example, we could imagine an application in
control theory were we want to check if the eigenvalues will lie in a specific
region, thus making the system stable or unstable. If we just have the
information from an approximating polygon in this application we would not be
able to say for certain if the limit set lies in a specific region or not,
since we do not know the size of the error. However, we could try to resolve
this issue by using the following idea: if we slightly expand every polygon in
(14) by an amount $\delta_{\rho}$, such that
$(\operatorname{sp}T(b_{\rho_{j}}^{D}))_{\delta_{\rho_{j}}}\supset\operatorname{sp}T(b_{\rho_{j}})$,
then the intersection of the expanded polygons would contain $\Lambda(b)$. The
intersection of the expanded polygons would then yield an approximating
superset of $\Lambda(b)$, which we know that $\Lambda(b)$ is a subset of. And
so if the approximating superset does not intersect a specific region, we know
for certain that $\Lambda(b)$ does not do that either.

######  Corollary 2.9.

Let $\rho_{l}$ and $\rho_{h}$ be bounds as in Theorem 2.1, let
$(\rho_{j})_{j=0}^{n}$ denote partitions of $[\rho_{l},\rho_{h}]$ with
granularity $\Delta_{\rho}$, and let $(v_{j})_{j=0}^{m}$ denote partitions of
$[0,2\pi]$ with granularity $\Delta_{v}$. Then there exists a family of
constants $C_{\rho}$ such that

$\begin{split}\bigcap_{j=0}^{n}\left(\operatorname{sp}T(b_{\rho_{j}}^{D})\right)_{C_{\rho_{j}}\Delta_{v}^{2}}&\supset\Lambda(b)\end{split}$

holds for all partitions of $[0,2\pi]$ with granularity less than
$\Delta_{v}$. Further, we have that

$\lim_{\Delta_{\rho},\Delta_{v}\to
0}\bigcap_{j=0}^{n}\left(\operatorname{sp}T(b_{\rho_{j}}^{D})\right)_{C_{\rho_{j}}\Delta_{v}^{2}}=\Lambda(b).$

###### Proof.

Let
$C_{\rho}:=\left\lVert{x}_{\rho}^{\prime\prime}(v)\right\rVert_{\infty}+\left\lVert{y}_{\rho}^{\prime\prime}(v)\right\rVert_{\infty}+\eta$
for arbitrary $\eta\in(0,1)$. From Proposition 2.5 and (13), we see that

$\left(\operatorname{sp}T(b_{\rho_{j}}^{D})\right)_{C_{\rho_{j}}\Delta_{v}^{2}}\supset\operatorname{sp}T(b_{\rho_{j}}),$

which shows the first statement about the inclusion.

For the second statement, let $C:=\sup_{\rho\in[\rho_{l},\rho_{h}]}C_{\rho}$
and note that

$\bigcap_{j=0}^{n}\left(\operatorname{sp}T(b_{\rho_{j}}^{D})\right)_{C_{\rho_{j}}\Delta_{v}^{2}}\subset\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D})\Biggr{)}_{C\Delta_{v}^{2}}.$

Now let $\Delta_{\rho}$ and $\Delta_{v}$ be such that
$d_{H}\bigl{(}\cap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\Lambda(b)\bigr{)}<\epsilon/2$
for partitions of granularities less than $\Delta_{\rho}$ and $\Delta_{v}$,
respectively. Further, suppose $\Delta_{v}$ fulfills
$C\Delta_{v}^{2}<\epsilon/2$. Then

$d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D}),\biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}^{D})\biggr{)}_{C\Delta_{v}^{2}}\,\Biggr{)}<\epsilon/2,$

and the triangle inequality for $d_{H}$ now gives what we want. ∎

######  Remark 2.10.

Using Corollary 2.9 we can easily modify Algorithm 1 to get a bound for
$\Lambda(b)$. Simply use the constants given in the proof of Corollary 2.9,
expand the polygons we intersect with by the given amount, and compute the
intersection in the same way as before.

##  3 Implementation details

The main operation we need in Algorithm 1 is polygon intersection, which is a
common problem in the field of computer graphics. Hence, it has been
thoroughly researched, and there exist efficient implementations. For this
work, we use the Python library pyclipper [1], which is a Python wrapper
library for the C$++$ library Clipper2 [14]. The library is based on Vatti’s
clipping algorithm, which is described in his article [18]. The processing
time for intersecting two polygons using Vatti’s algorithm scales linearly in
the number of vertices of the polygons being intersected.

To produce the bound for $\Lambda(b)$ as described in Remark 2.10, we also
need to be able to expand polygons. This problem is called polygon offsetting
in the literature, and it is a well researched problem. In fact, Clipper2
contains an implementation of the offsetting algorithm as described in [8]. It
should be noted that we do not want to compute
$(\operatorname{sp}T(b_{\rho}^{D}))_{C_{\rho}\Delta_{v}^{2}}$ exactly, since
the boundary then contains circular arcs, which are more difficult to
intersect with than polygons. Instead, we bound
$(\operatorname{sp}T(b_{\rho}^{D}))_{C_{\rho}\Delta_{v}^{2}}$ with a polygon
${P}_{\rho}$ such that

$\left((\operatorname{sp}T(b_{\rho}^{D})\right)_{C_{\rho}\Delta_{v}^{2}}\subset{P}_{\rho}\subset\left(\operatorname{sp}T(b_{\rho}^{D})\right)_{2C_{\rho}\Delta_{v}^{2}}$
(17)

The constant 2 here is chosen quite arbitrarily, and any constant would work
in theory. For practical reasons we want the chosen constant to be “large
enough” so that $P_{\rho}$ does not contain too many vertices, and “small
enough” so that $P_{\rho}$ is close to $\Lambda(b)$. Because of the way
$C_{\rho}$ is chosen in the proof of Corollary 2.9, the conclusions hold also
for $2C_{\rho}$, and because of (17), they then also hold for ${P}_{\rho}$.
Clipper2 contains an implementation for choosing such ${P}_{\rho}$. The time
complexity for these offset polygons is a bit tricky to analyze, as the
offsetting algorithm has time complexity $O((m+k)\log m)$ where $m$ is the
number of vertices of the input polygon and $k$ the number of self-
intersections in the raw offset curve [8]. Intuition tells us that for
sufficiently regular polygons, such as the ones we work with, $k$ should be of
the same order as $m$, so in the following analysis of time complexity we
assume that producing ${P}_{\rho}$ takes $O(m\log m)$ time and that the number
of vertices for ${P}_{\rho}$ is $O(m)$.

Using the above implementations of offsetting and polygon intersection we can
compute the approximating polygon and approximating superset of $\Lambda(b)$
given by Algorithm 1. We do this for the symbol presented in Example 2. The
results can be seen in Figure 4. If we look at the bigger picture as in Figure
4(a), we see that the results are very similar to those of the algebraic
approach. And if we zoom in, we see that our geometric approach struggles at
the same places as the algebraic one. However, Figure 4(d) shows clearly that
$\Lambda(b)$ actually separates the complex plane, which was not clear from
the results of the algebraic approach. It took roughly 3 minutes to run the
algorithm for finding the approximating polygon and the approximating
superset.

(a) An approximating polygon for $\Lambda(b)$ computed using Algorithm 1 with
$n=1250$ and $m=1000$.

(b) A zoomed in version of 4(a).

(c) An approximating superset of $\Lambda(b)$ computed using Algorithm 1 with
$n=1250$ and $m=1000$.

(d) A zoomed in version of 4(c). Blue polygons are positively oriented and the
red ones are negatively oriented. That is, the approximating superset is the
region outside the red polygons and inside the blue polygon.

Figure 4: A comparison between Algorithm 1 and the algebraic approach
presented in 1.1. The symbol $b$ being investigated is presented in (6). The
lime dots in all of the subfigures are the results from the algebraic
approach, these are the same results as presented in Figure 2.

Besides finding efficient polygon intersection and offsetting algorithms, the
most important issue in the implementation is how to choose
$(v_{j})_{j=0}^{m}$ and $(\rho_{j})_{j=0}^{n}$ (i.e., where to choose the
vertices of $b_{\rho_{j}}^{D}$ on $b_{\rho_{j}}(\mathbb{T})$ and what $\rho$’s
to choose for the intersections). We know from Theorem 2.7 that if we just
choose the granularities to be finer and finer, then the intersection will
approach $\Lambda(b)$. But we should be able to get a more computationally
efficient algorithm if the partitions are chosen in a good way.

For simplicity, the only investigated partition for $v$ is the uniform one,
that is, $v_{j}=2\pi j/m$. It is quite difficult to imagine some strategy
working much better. Choosing $\rho$’s on the other hand is more interesting.
A natural approach is to first calculate $[\rho_{l},\rho_{h}]$ using the
estimate from Proposition 2.1. Then, a naive solution is to sample $\rho$’s
uniformly in this interval.

There are two major points of critique for this method of sampling $\rho$’s.
First of all, the bounds $[\rho_{l},\rho_{h}]$ are not strict. For the
$\rho$’s that are close to the end points, $\operatorname{sp}T(b_{\rho})$
probably contains the entirety of the currently intersected polygon. Secondly,
from the proof of Theorem 2.4 we see that ideally, we want to choose more
$\rho$’s where “the Hausdorff metric of $\operatorname{sp}T(b_{\rho})$ varies
rapidly with $\rho$”.

To get an intuition of good heuristics for choosing stricter bounds for
$\rho$, it is helpful to plot $b_{\rho}(\mathbb{T})$ and let $\rho$ vary. An
insight from experimenting is that $\Lambda(b)$ is formed by the points where
$b_{\rho}(\mathbb{T})$ intersects itself, and when $\rho$ is varied, this
intersection moves along $\Lambda(b)$ in a continuous manner. This would
suggest that we only want to find the intervals for which
$b_{\rho}(\mathbb{T})$ self-intersect, and where these intersections are a
part of $\Lambda(b)$. As these intersections move continuously with $\rho$,
there should be a finite number of interesting intervals for $\rho$. Actually,
one can see from (4) that all $\lambda\in\Lambda(b)$, except for a finite
number, must arise from a self-intersection of $b_{\rho}(\mathbb{T})$, since
(4) exactly describes this type of intersection.

To find the interesting intervals of $\rho$’s analytically is not obvious. But
we can use the fact that when we intersect our partial polygon with
$\operatorname{sp}T(b_{\rho})$ for a good $\rho$, the area of the partial
polygon decreases. So we could sweep over $[\rho_{l},\rho_{h}]$ and calculate
the decrease in area when intersecting with different
$\operatorname{sp}T(b_{\rho})$. Also, note that the area of a polygon can be
calculated in linear time with respect to the number of vertices. Pseudocode
for calculating the approximating polygon and approximating superset of
$\Lambda(b)$ with the sweep improvement can be seen in Algorithm 2. How to
choose the intervals of $\rho$’s where the area is reduced the most deserves
some attention. If some $\rho$ is in both rhos and sweeprhos in Algorithm 2,
the area reduced will be 0. To avoid missing out on good $\rho$’s because of
this, a moving average filter is used for areareduce, and then all the good
intervals are selected as the ones where areareduce[$\rho$]
$\geq(\max_{\rho}\text{{areareduce}}[\rho])/10^{6}$.

The problem of sampling the same $\rho$ also applies when we sample rhos for
the next intersection-run. Typically, we will find the exact same interval of
good $\rho$’s that reduce area the most in each iteration. Therefore, a
randomly chosen $\epsilon$ with
$0\leq\epsilon\leq\text{sweeprhos}[i+1]-\text{sweeprhos}[i]$ is added to the
left boundaries of the good intervals and subtracted from the right
boundaries.

A significant optimization can be done when we calculate areareduce. The
important observation is that areareduce$[\rho]$ is non-increasing. When
updating areareduce we begin by calculating $\Lambda_{s}$ for the intervals of
$\rho$ that we just intersected with. Then, if areareduce for the $\rho$ that
did not reduce area the most in the previous iteration still are smaller than
the threshold, we don’t need to update those values. Practically, this
drastically reduces the number of intersections done in total. Most polygon
intersections used for calculating $\Lambda_{s}$ are done in the first
iteration, when we have no prior values in areareduce.

The second major point of critique to the naive method means that we ideally
want to sample $\rho$’s denser where “the Hausdorff metric of
$\operatorname{sp}T(b_{\rho})$ varies rapidly with $\rho$”. We tried a
strategy of doing this. The idea behind our strategy is the observation that
$b$ can be seen as a sum of two polynomials, one in $z$ and one in $1/z$. The
standard approach if nothing better is known is to sample uniformly, but
because of $1/z$, it would make sense to sample $\rho_{j}$ uniformly for
$\rho\geq 1$ and sampling $\rho<1$ in such a way that $1/\rho_{j}$ becomes
uniformly distributed.

We note that besides the above mentioned improvements, Algorithm 2 also
contains the details on how to produce the approximating superset, which is
denoted $\Lambda_{b}$ in the pseudocode.

Algorithm 2 Improved $\Lambda(b)$–calculator using area sweeps

procedure CalcLimitSet($b,n,m,l,\\#sweeps$)

$b$: the symbol,

$n$: the number of sampled $\rho$’s,

$m$: number of sampled $v$’s,

$l$: the number of sampled $\rho$’s for area sweeping,

$\\#sweeps$: how many times new $\rho$’s should be generated.

$\rho_{l},\rho_{h}\leftarrow$ bounds from Theorem 2.1

sweeprhos $\leftarrow$ sample $l$ points for area sweeping in
$[\rho_{l},\rho_{h}]$.

rhos $\leftarrow$ sample $n/\\#sweeps$ points in $[\rho_{l},\rho_{h}]$

vs $\leftarrow$ sample $m+1$ points in $[0,2\pi]$

$\Lambda\leftarrow$ $b_{\text{rhos[1]}}^{D}(vs)$

$\Lambda_{b}\leftarrow$
$(b_{\text{rhos[1]}}^{D}(vs))_{C_{\rho}\Delta_{v}^{2}}$

for $i\leftarrow 1$ to $\\#sweeps$ do

for each $\rho$ in rhos do

$\Lambda\leftarrow\Lambda\cap b_{\rho}^{D}(vs)$

$\Lambda_{b}\leftarrow\Lambda_{b}\cap(b_{\rho}^{D}(vs))_{C_{\rho}\Delta_{v}^{2}}$

end for

for each $\rho$ in sweeprhos do

$\Lambda_{s}\leftarrow\Lambda\cap b_{\rho}^{D}(vs)$

areareduce[$\rho$] $\leftarrow$ Area($\Lambda$) - Area($\Lambda_{s}$)

end for

rhos $\leftarrow$ sample $n/\\#sweeps$ points in the intervals that reduces
area the most

end for

return $\Lambda,\Lambda_{b}$

end procedure

###  3.1 Testing suggested improvements

Regarding the suggested improvement about sampling $\rho$’s denser where the
Hausdorff metric varies rapidly, we found that the ideal strategy for sampling
$\rho$ is highly dependent on the symbol. So in practical applications where
high performance is required, the user should fine tune the sampling strategy
for the specific symbol.

To investigate the other suggested improvements (i.e., using Algorithm 2
instead of Algorithm 1), we again compute an approximating polygon and
approximating superset of $\Lambda(b)$ for the symbol given in Example 2. To
get a fair comparison, we run Algorithm 2 with $n=1000$, $m=1000$, $l=250$ and
$\\#sweeps=2$, which means that we use the same number of polygon
intersections for both algorithms.

The approximating polygon and approximating superset of $\Lambda(b)$ computed
with Algorithm 2 can be seen in Figure 5. We cannot really see any big
difference compared to Figure 4 in the zoomed out pictures. But if we look at
Figures 5(b) and 5(d), we see that the approximating polygon looks more like
the type of set we are searching for (a connected finite union of analytic
arcs), and the approximating superset in 5(d) is tighter than the one in 4(d).
These properties should generalize since it is never good to “waste” polygon
intersections. It took roughly 3 minutes to run the algorithm for finding the
approximating polygon and approximating superset.

(a) An approximating polygon for $\Lambda(b)$ computed using Algorithm 2 with
$n=1000$, $m=1000$, $l=250$, $\\#sweeps=2$.

(b) A zoomed in version of 5(a).

(c) An approximating superset of $\Lambda(b)$ computed using Algorithm 2 with
$n=1000$, $m=1000$, $l=250$, $\\#sweeps=2$.

(d) A zoomed in version of 5(c). Blue polygons are positively oriented and the
red ones are negatively oriented. That is, the approximating superset is the
region outside the red polygons and inside the blue polygon.

Figure 5: The results of Algorithm 2. The symbol $b$ being investigated is
presented in (6). The lime dots in all of the subfigures are the results from
the algebraic approach, these are the same results presented in Figure 2.

##  4 Time complexity and convergence speed

Analyzing the time complexity and convergence speed of Algorithm 2 rigorously
seems difficult, but we will provide an argument and an example to back up
that the running time for Algorithm 2 on average probably is $O(n^{2}+nm\log
m)$. Using that bound, we then argue that the distance from the approximating
polygon to $\Lambda(b)$ for most of $\Lambda(b)$ decreases as $O(1/\sqrt{k})$
where $k$ is the number of elementary operations.

###  4.1 Time complexity

We wish to analyze the time complexity of Algorithm 2, and to simplify that
analysis we begin by looking at only the part of the Algorithm that produces
the approximating polygon for $\Lambda(b)$. The most computationally complex
parts of Algorithm 2 are the polygon intersections which we know are linear in
the number of vertices of the incoming polygons. So what we really want to
analyze is how the number of vertices in the approximating polygon $\Lambda$
grows. Let $w(i)$ be the number of vertices of $\Lambda$ after $i$
intersections. We will now argue that $w(i)$ should increase at most linearly.
To see this, consider Figure 6. When we intersect $\Lambda$ with
$b_{\rho}^{D}(vs)$, the typical situation is described there. Specifically,
Figure 6 is zoomed in where $b_{\rho}^{D}(vs)$ intersects itself and the
intersection is part of $\Lambda(b)$, which happens a bounded number of times.
In Figure 6 we see that 3 vertices are added to $\Lambda$ for the specific
self-intersection of $b_{\rho}^{D}(vs)$. Since the number of self-
intersections of $b_{\rho}^{D}(vs)$ is bounded, $w(i)$ should increase at most
linearly.

Limit set Figure 6: A probable illustration of how it looks when $\Lambda$ is
intersected with $b_{\rho}^{D}(vs)$ in Algorithm 2. The red line is the true
$\Lambda(b)$, the black lines are $\Lambda$ prior to being intersected, and
the blue lines are $b_{\rho}^{D}(vs)$. The figure is zoomed in where
$b_{\rho}^{D}(vs)$ intersects itself. The green points mark the new vertices
in $\Lambda$, and the red ones are the vertices being removed from $\Lambda$.

To test our hypothesis about $w(i)$ growing linearly, we generate a random
Laurent polynomial with $r=s=5$. We get the polynomial

$\begin{multlined}b(t)=(9.563-6.844i)t^{5}+(-2.246-4.182i)t^{4}+(0.8478+0.4667i)t^{3}\\\
+(-4.125+1.539i)t^{2}+(9.114+2.07i)t+(6.679+7.744i)\\\
+(5.21-3.819i)t^{-1}+(-4.62-3.149i)t^{-2}+(-6.911-1.115i)t^{-3}\\\
+(4.486+4.44i)t^{-4}+(4.279+0.9657i)t^{-5}.\end{multlined}b(t)=(9.563-6.844i)t^{5}+(-2.246-4.182i)t^{4}+(0.8478+0.4667i)t^{3}\\\
+(-4.125+1.539i)t^{2}+(9.114+2.07i)t+(6.679+7.744i)\\\
+(5.21-3.819i)t^{-1}+(-4.62-3.149i)t^{-2}+(-6.911-1.115i)t^{-3}\\\
+(4.486+4.44i)t^{-4}+(4.279+0.9657i)t^{-5}.$ (18)

We then run Algorithm 2 with $m=500$, $n=3000$, $l=500$, $\\#sweeps=6$, $\rho$
sampled uniformly, and keep track of $w(i)$, and plot it in Figure 7. The
resulting approximating polygon can be seen in Figure 8. Figure 7 looks like
it should according to our argument: $w(i)$ appears to grow piecewise
linearly, with the slope changing. These changes correspond to the alterations
of the amount of self-intersections. The algorithm finds the same interval of
good $\rho$’s after each intersection-run, which explains the repetition of
the pattern. Also note that the repetitive pattern we see is exactly what
should happen if our reasoning in Figure 6 is correct, since we intersect over
the same interval of $\rho$’s each time.

Figure 7: The number of vertices of $\Lambda$ in Algorithm 2 as a function of
the number of intersections. Figure 8: The approximating polygon for
$\Lambda(b)$ given by Algorithm 2 for $b$ as in (18), $m=500$, $n=3000$,
$l=500$, $\\#sweeps=6$.

Assuming $w(i)=O(i+m)$, the total complexity is

$\begin{split}&\sum_{i=1}^{n}O\bigl{(}w(i)+m\bigr{)}+\sum_{s=1}^{\\#sweeps}\sum_{i=1}^{l}O\bigr{(}w(s\cdot
n/\\#sweeps)+m\bigr{)}\\\
&=O\left(\frac{n(1+2m+n+2m)}{2}\right)+O\left(\frac{l\\#sweeps(n/\\#sweeps+2m+n+2m)}{2}\right)\\\
&=O(n^{2}+nm+nl\\#sweeps+ml\\#sweeps).\end{split}$

However, as we remarked earlier, we can practically optimize away the
$\\#sweeps$-factor if it is of moderate size, and in regular use cases it
should never be set higher than about 10. Also, typically $l$ is chosen as
smaller than $n$, else we would “waste” too many intersections on area
sweeping. So with these reasonable restrictions on $\\#sweeps$ and $l$, the
time complexity for computing the approximating polygon in Algorithm 2 should
on average be $O(n^{2}+nm)$.

We now analyze the time complexity of Algorithm 2 when we also produce the
approximating superset. Remember the earlier assumptions, that the number of
vertices for ${P}_{\rho}$ is $O(m)$ and producing ${P}_{\rho}$ takes $O(m\log
m)$ time. Using a similar argument as above, we could argue that the number of
vertices for $\Lambda_{b}$ in Algorithm 2 grows linearly. So the time
complexity for computing the approximating superset is almost the same as
computing the approximating polygon, the only operation that adds additional
time complexity is computing ${P}_{\rho}$. We compute $n$ ${P}_{\rho}$’s in
Algorithm 2 so this step adds an additional $O(nm\log m)$ term to the overall
time complexity. Hence the total time complexity for computing the
approximating superset using Algorithm 2 is $O(n^{2}+mn)+O(nm\log
m)=O(n^{2}+mn\log m)$. Since computing the approximating polygon took
$O(n^{2}+mn)$ time, the total time complexity for running Algorithm 2 is
$O(n^{2}+mn\log m)$.

###  4.2 Convergence speed

To argue about the convergence speed we first remind ourselves about Remark
2.8 which bounds the part of the error originating in the polygon
approximation. Because of (16), we have a good idea of how big we must choose
$m$ to get a certain error level. Therefore we will consider $m$ to be
constant and reason about the convergence speed as $n$ increases. This
approach also makes sense in terms of the algorithm, $m$ is set in the
beginning and after that we can do as many intersections as we want.

From (16) we see that what we want to investigate is

$d(n):=d_{H}\Biggl{(}\,\bigcap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\Lambda(b)\Biggr{)}.$

So $d(n)$ describes how the error changes when $n$ increases. Consider Figure
6 but think of the black lines describing
$\cap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}})$, which makes sense since if
$b^{\prime}_{\rho}(t)\neq 0$, then locally $b_{\rho}(t)$ is almost linear.
Note that the distances between the red curve and the red dots in Figure 6 is
a lower bound of $d(n)$. If we increase $n$ by one, and intersect with the
spectrum corresponding to the blue lines, we argue that the now relevant
distances, the ones between the green dots and the red curve, should decrease
geometrically. Assuming this is true for all of $\Lambda(b)$, we see that as
$n$ is doubled, the maximum distance $d(n)$ should have decreased
geometrically, since doubling $n$ corresponds to intersecting with a “blue
spectrum” between each pair of black self-intersections, as in Figure 6, and
each black self-intersection corresponds to at most one $\rho$. So we have
argued that

$d(2n)\leq d(n)/r,$ (19)

for some factor $r>1$. If we assume the new intersections to be placed in the
“middle of the old ones”, we should have $r\approx 2$.

However, ultimately we are interested in how $d(n)$ decreases as a function of
the number of elementary operations. Let $k$ denote the number of elementary
operations. As we argued in the previous subsection, $k=O(n^{2}+mn\log m)$,
and we assumed $m$ to be constant, so $k<(c^{*})^{2}n^{2}$, which means
$n>\sqrt{k}/c^{*}$. Let $d^{*}(k)$ denote
$d_{H}(\cap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\Lambda(b))$ as a
function of the number of elementary operations. Further, we see that $d(n)$
and $d^{*}(k)$ are non-increasing, and get

$d^{*}(k)\leq d\left(\left\lfloor\sqrt{k}/c^{*}\right\rfloor\right)\leq
d\left(n_{0}2^{\left\lfloor\log_{2}\left(\left\lfloor\frac{\sqrt{k}}{c^{*}}\right\rfloor/n_{0}\right)\right\rfloor}\right),$
(20)

for some $n_{0}$, which is the $n$ at which we assume we begin our doubling
argument. Using the bounds in (20) and (19) recursively, we obtain

$\begin{multlined}d^{*}(k)\leq\frac{d(n_{0})}{r^{\left\lfloor\log_{2}\left(\left\lfloor\frac{\sqrt{k}}{c^{*}}\right\rfloor/n_{0}\right)\right\rfloor}}\leq\frac{d(n_{0})}{2^{(\log_{2}r)\left(\log_{2}\left(\left\lfloor\frac{\sqrt{k}}{c^{*}}\right\rfloor/n_{0}\right)-1\right)}}\\\
=\frac{d(n_{0})}{\left(\left\lfloor\frac{\sqrt{k}}{c^{*}}\right\rfloor/2n_{0}\right)^{\log_{2}r}}\leq\frac{c}{k^{\frac{1}{2}\log_{2}r}}.\end{multlined}d^{*}(k)\leq\frac{d(n_{0})}{r^{\left\lfloor\log_{2}\left(\left\lfloor\frac{\sqrt{k}}{c^{*}}\right\rfloor/n_{0}\right)\right\rfloor}}\leq\frac{d(n_{0})}{2^{(\log_{2}r)\left(\log_{2}\left(\left\lfloor\frac{\sqrt{k}}{c^{*}}\right\rfloor/n_{0}\right)-1\right)}}\\\
=\frac{d(n_{0})}{\left(\left\lfloor\frac{\sqrt{k}}{c^{*}}\right\rfloor/2n_{0}\right)^{\log_{2}r}}\leq\frac{c}{k^{\frac{1}{2}\log_{2}r}}.$

for some constant $c$. The conclusion of the above argument is that for a
fixed $m$ we should have $d^{*}(k)\leq O(1/k^{\frac{1}{2}\log_{2}r})$. And if
we – as commented on earlier – choose the $\rho$’s in a good way, $r$ should
be $\approx 2$, which means $d^{*}(k)\leq O(1/\sqrt{k})$.

So we have argued for
$d_{H}(\cap_{j=0}^{n}\operatorname{sp}T(b_{\rho_{j}}),\Lambda(b))$ decreasing
as $O(1/\sqrt{k})$, and assuming that $\Delta_{v}$ is small enough, this
argument also carries over to the distances between the approximating polygon
and $\Lambda(b)$, and between the approximating superset and $\Lambda(b)$.
This can be seen from (16) and the proof of Corollary 2.9.

However, this entire argument is based on the assumption that Figure 6 is a
good approximation of the local behavior of intersected spectra. An implicit
assumption in Figure 6 is for example that $b_{\rho}(\mathbb{T})$ is not
tangent to itself (which is true for all but a finite number of points on
$\Lambda(b)$). If one examines the symbol in Example 2 more closely, one sees
that this is exactly what happens in the center of “the cross”, since

$b_{1}(e^{i0})=b_{1}(e^{i\frac{\pi}{2}})=0,\;\;\frac{\mathop{}\\!{d}}{\mathop{}\\!{d}v}b_{1}(e^{iv})\Bigr{|}_{v=0}=\frac{\mathop{}\\!{d}}{\mathop{}\\!{d}v}b_{1}(e^{iv})\Bigr{|}_{v=\frac{\pi}{2}}=-(1+i).$

##  5 Conclusion

A new approach to approximating $\Lambda(b)$ geometrically has been presented.
Pseudocode for the approach can be found in Algorithm 2. Compared to the
previous approaches that are algebraically based, it has the advantage of not
requiring root finding of arbitrary complex polynomials. Furthermore, one
could argue that the output from the geometric algorithm is more natural than
the output for the algebraic one, since you get a connected subset of the
complex plane instead of points sampled on $\Lambda(b)$. In some cases, the
geometric algorithm yields better information on $\Lambda(b)$. An example of
this can be seen in Figure 5(d) where it is not clear from the lime dots
whether or not $\Lambda(b)$ separates the plane, but from the approximating
superset we see that $\Lambda(b)$ does split the plane.

## Acknowledgments

The authors wish to thank Mikael Persson Sundqvist for fruitful and insightful
discussions and Frank Wikström for suggesting the approach of constructing an
approximating superset.

J. S. Christiansen’s research was supported in part by the Swedish Research
Council (VR) Grant No. 2018-03500.

## References

  * [1] Pyclipper.  https://github.com/fonttools/pyclipper.  Accessed: 2022-11-24.
  * [2] R. M. Beam and R. F. Warming.  The asymptotic spectra of banded Toeplitz and quasi-Toeplitz matrices.  SIAM J. Sci. Comput. 14 (1993), pp. 971–1006.
  * [3] A. Böttcher, J. Gasca, S. M. Grudsky, and A. V. Kozak.  Eigenvalue Clusters of Large Tetradiagonal Toeplitz Matrices.  Integral Equations Operator Theory 93 (2021).
  * [4] A. Böttcher and S. M. Grudsky.  Spectral properties of banded Toeplitz matrices. Society for Industrial and Applied Mathematics, Philadelphia, 2005.
  * [5] A. Böttcher, S. M. Grudsky, and I. M. Spitkovsky.  The spectrum is discontinuous on the manifold of Toeplitz operators. Arch. Math. (Basel) 75 (2000), pp. 46–52.
  * [6] A. Böttcher and B. Silbermann.  Introduction to large truncated Toeplitz matrices.  Springer New York, NY, 1999.
  * [7] E. Čech and M. Katětov.  Point Sets.  Academia, Publishing House of the Czechoslovak Academy of Sciences, Prague, 1969.
  * [8] X. Chen and S. McMains.  Polygon offsetting by computing winding numbers,  in International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, volume 4739, pp. 565–575, 2005.
  * [9] M. Duits and A. B. J. Kuijlaars.  An equilibrium problem for the limiting eigenvalue distribution of banded Toeplitz matrices.  SIAM J. Matrix Anal. Appl. 30 (2008), pp. 173–196.
  * [10] I. Gohberg.  On an application of the theory of normed rings to singular integral equations. Uspekhi Mat. Nauk 7 (1952), pp. 149–156.
  * [11] I. Gohberg.  On the number of solutions of homogeneous singular equations with continuous coefficients. Dokl. Akad. Nauk 112 (1958), pp. 327–330.
  * [12] I. I. Hirschman.  The spectra of certain Toeplitz matrices. Illinois J. Math. 11 (1967), pp. 145–159.
  * [13] I. S. Hwang and W. Y. Lee.  On the continuity of spectra of Toeplitz operators. Arch. Math. (Basel) 70 (1998), pp. 66–73.
  * [14] A. Johnson.  Clipper2 - Polygon Clipping and Offsetting Library.  http://www.angusj.com/clipper2/Docs/Overview.htm.  Accessed: 2022-11-24.
  * [15] T. Katō.  Perturbation theory for linear operators. Springer Berlin, Heidelberg, 1995.
  * [16] P. Schmidt and F. Spitzer.  The Toeplitz matrices of an arbitrary Laurent polynomial. Math. Scand. 8 (1960), pp. 15–38.
  * [17] J. L. Ullman.  A problem of Schmidt and Spitzer. Bulletin of the American Mathematical Society 73 (1967), pp. 883–885.
  * [18] B. R. Vatti.  A generic solution to polygon clipping. Communications of the ACM 35 (1992), pp. 56–63.
  * [19] H. Widom.  Eigenvalue distribution of nonselfadjoint Toeplitz matrices and the asymptotics of Toeplitz determinants in the case of nonvanishing index, in  Operator Theory: Advances and Applications 48 (1990), pp. 387–421.
  * [20] H. Widom.  Eigenvalue Distribution for Nonselfadjoint Toeplitz Matrices,  in Toeplitz Operators and Related Topics, E. L. Basor and I. Gohberg, eds., Birkhäuser Basel., Basel, 1994, pp. 1–8.