# Fast Toeplitz eigenvalue computations, joining interpolation-extrapolation
matrix-less algorithms and simple-loop theory

M. Bogoya<EMAIL_ADDRESS>University of Insubria, Como,
Italy.
Dipartimento di Scienza e Alta Tecnologia. S.E. Ekström sven-
<EMAIL_ADDRESS>Uppsala University, Uppsala, Sweden.
Department of Information Technology, Division of Scientific Computing.  S.
Serra–Capizzano<EMAIL_ADDRESS>University of Insubria, Como,
Italy.
Dipartimento di Scienza e Alta Tecnologia.

###### Abstract

Under appropriate technical assumptions, the simple-loop theory allows to
deduce various types of asymptotic expansions for the eigenvalues of Toeplitz
matrices generated by a function $f$. Independently and under the milder
hypothesis that $f$ is even and monotonic over $[0,\pi]$, matrix-less
algorithms have been developed for the fast eigenvalue computation of large
Toeplitz matrices, within a linear complexity in the matrix order: behind the
high efficiency of such algorithms there are the expansions predicted by the
simple-loop theory, combined with the extrapolation idea.

Here we focus our attention on a change of variable, followed by the
asymptotic expansion of the new variable, and we adapt the matrix-less
algorithm to the considered new setting.

Numerical experiments show a higher precision (till machine precision) and the
same linear computation cost, when compared with the matrix-less procedures
already presented in the relevant literature. Among the advantages, we
concisely mention the following: a) when the coefficients of the simple-loop
function are analytically known, the algorithm computes them perfectly; b)
while the proposed algorithm is better or at worst comparable to the previous
ones for computing the inner eigenvalues, it is extremely better for the
computation of the extreme eigenvalues.
Keywords: Eigenvalue computation, Toeplitz matrix, Matrix-less method,
Asymptotic expansion.

MSC Classes: Primary 15B05, 65F15, 65D05, 47B35. Secondary 15A18, 47A38.

##  1 Introduction

The target of this note is to design fast procedures for the computation of
all the spectra of large Toeplitz matrices having an even generating function
which is monotone in the interval $(-\pi,\pi)$. For the formal definition of
Toeplitz matrix generated by a Lebesgue integrable function over the basic
interval $Q\equiv(-\pi,\pi)$ see the first lines of Section 2.

This topic has been studied in the recent years by several researchers.
Indeed, taking into account a clear numerical evidence developed in a
systematic series of the numerical tests, in [15] the second author formulated
the following conjecture.

######  Conjecture 1.1.

Let $l,g$ be two real-valued even functions with $g>0$ on $(0,\pi)$, and
suppose that $f\equiv\frac{l}{g}$ is monotone increasing over $(0,\pi)$. Set
$X_{n}\equiv T_{n}^{-1}(g)T_{n}(l)$ for all $n$. Then, for every integer
$\alpha\geqslant 0$, every $n$, and every $j=1,\ldots,n$, the following
asymptotic expansion holds:

$\lambda_{j}(X_{n})=f(\theta_{j,n})+\sum_{k=1}^{\alpha-1}c_{k}(\theta_{j,n})h^{k}+E_{j,n,\alpha},$

where:

  * •

the eigenvalues of $X_{n}$ are arranged in non-decreasing order,
$\lambda_{1}(X_{n})\leqslant\cdots\leqslant\lambda_{n}(X_{n})$; 111Note that
the eigenvalues of $X_{n}$ are real, because $T_{n}(g)$ is symmetric positive
definite and $X_{n}$ is similar to the symmetric matrix
$T_{n}^{-\frac{1}{2}}(g)T_{n}(l)T_{n}^{-\frac{1}{2}}(g)$.

  * •

$\\{c_{k}\\}_{k=1}^{\infty}$ is a sequence of functions from $(0,\pi)$ to
$\mathbb{R}$ which depends only on $l,g$;

  * •

$h\equiv\frac{1}{n+1}$ and $\theta_{j,n}\equiv\frac{j\pi}{n+1}=j\pi h$;

  * •

$E_{j,n,\alpha}=O(h^{\alpha})$ is the remainder (the error), which satisfies
the inequality $|E_{j,n,\alpha}|\leqslant\kappa_{\alpha}h^{\alpha}$ for some
constant $\kappa_{\alpha}$ depending only on $\alpha,l,g$.

In the case where $g=1$ identically, Conjecture 1.1 was originally formulated
and supported through numerical experiments in [15]. Then the algorithmic
proposal was extended and refined in [1, 12, 13, 14]. When $g\equiv 1$ and $l$
satisfies further technical additional assumptions, those of the simple-loop
method, Conjecture 1.1 was formally proved by Bogoya, Böttcher, Grudsky, and
Maximenko in a series of papers [5, 6, 9].

However, again in the case of $g\equiv 1$, the power of the simple-loop method
has been not exploited completely, since in the case of a continuous
generating function, the distribution results reported in Theorem 2.2 imply
that $\lambda_{j}(T_{n}(f))=f(s_{j,n})$ with $s_{j,n}$ belonging to $(0,\pi)$
and distributed as the identity function. However, more is known and indeed
also in the independent variable $\theta$ there exists an asymptotic expansion
regarding exactly the points $s_{j,n}$.

This note deals with the adaptation of the interpolation-extrapolation
algorithms to the previous change of variable, joined with a trick at the end
points introduced in [8].

The numerical results are extremely precise, even compared with the already
good performances described in [1, 12, 13, 14, 15], since it is not difficult
to reach machine precision, and the complexity is still linear.

The present note is organized as follows. Preliminary definitions, tools, and
results are concisely reported in Section 2. Section 3 presents the new
adapted algorithm for computing the Toeplitz eigenvalues: as in [14], our
technique combines the extrapolation procedure proposed in [1, 15] – which
allows the computation of some of the eigenvalues of $X_{n}$ – with an
appropriate interpolation process, designed for the simultaneous computation
of all the eigenvalues of $X_{n}$, with the additional end point trick in [8].
In Section 4 we present the numerical experiments, while in Section 5 we draw
conclusions and we list few open problems for research lines to be
investigated in the next future.

##  2 Preliminaries and Tools

For a real or complex valued function $f$ in $L^{1}[-\pi,\pi]$, let
$\mathfrak{a}_{j}(f)$ be its $j$th Fourier coefficient, i.e.

$\mathfrak{a}_{j}(f)\equiv\frac{1}{2\pi}\int_{-\pi}^{\pi}f(\theta)\mathrm{e}^{-\mathrm{i}j\theta}\operatorname{d\\!}\theta,\quad
j\in\mathbb{Z},$

and consider the sequence $\\{T_{n}(f)\\}_{n=1}^{\infty}$ of the $n\times n$
Toeplitz matrices defined by
$T_{n}(f)\equiv(\mathfrak{a}_{j-k}(f))_{j,k=0}^{n-1}$. The function $f$ is
customarily referred to as the generating function of this sequence.

As a second step, we introduce some notations and definitions concerning
general sequences of matrices. For any function $F$ defined on the complex
field and for any matrix $A_{n}$ of size $d_{n}$, by the symbol
$\Sigma_{\lambda}(F,A_{n})$, we denote the mean

$\Sigma_{\lambda}(F,A_{n})=\frac{1}{d_{n}}\sum_{j=1}^{d_{n}}F[\lambda_{j}(A_{n})],$

while by the symbol $\Sigma_{\sigma}(F,A_{n})$, we denote the mean

$\Sigma_{\sigma}(F,A_{n})=\frac{1}{d_{n}}\sum_{j=1}^{d_{n}}F[\sigma_{j}(A_{n})].$

######  Definition 2.1.

Given a sequence $\\{A_{n}\\}$ of matrices of size $d_{n}$ with
$d_{n}<d_{n+1}$ and given a Lebesgue-measurable function $\psi$ defined over a
measurable set $K\subset{\mathbb{R}}^{\nu}$, $\nu\in\mathbb{N}^{+}$, of finite
and positive Lebesgue measure $\mu(K)$, we say that $\\{A_{n}\\}$ is
distributed as $(\psi,K)$ in the sense of the eigenvalues if for any
continuous function $F$ with bounded support, the following limit relation
holds

$\lim_{n\rightarrow\infty}\Sigma_{\lambda}(F,A_{n})=\frac{1}{\mu(K)}\int_{K}F(\psi)\,\operatorname{d\\!}\mu.$

In this case, we write in short $\\{A_{n}\\}\sim_{\lambda}(\psi,K)$.
Furthermore we say that $\\{A_{n}\\}$ is distributed as $(\psi,K)$ in the
sense of the singular values if for any continuous function $F$ with bounded
support, the following limit relation holds

$\lim_{n\rightarrow\infty}\Sigma_{\sigma}(F,A_{n})=\frac{1}{\mu(K)}\int_{K}F(|\psi|)\,\operatorname{d\\!}\mu.$

In this case, we write in short $\\{A_{n}\\}\sim_{\sigma}(\psi,K)$, which is
equivalent to $\\{A_{n}^{*}A_{n}\\}\sim_{\lambda}(|\psi|^{2},K)$.

In Remark 2.1 we provide an informal meaning of the notion of eigenvalue
distribution. For the singular value distribution similar statements can be
written.

######  Remark 2.1.

The informal meaning behind the above definition is the following. If $\psi$
is continuous, $n$ is large enough, and

$\Big{\\{}{\bf x}_{j}^{(d_{n})},\ j=1,\ldots,d_{n}\Big{\\}}$

is an equispaced grid on $K$, then a suitable ordering $\lambda_{j}(A_{n})$,
$j=1,\ldots,d_{n}$, of the eigenvalues of $A_{n}$ is such that the pairs
$\big{\\{}\big{(}{\bf x}_{j}^{(d_{n})},\lambda_{j}(A_{n})\big{)},\
j=1,\ldots,d_{n}\big{\\}}$ reconstruct approximately the hypersurface

$\big{\\{}({\bf x},\psi({\bf x})),\ {\bf x}\in K\big{\\}}.$

In other words, the spectrum of $A_{n}$ ‘behaves’ like a uniform sampling of
$\psi$ over $K$. For instance, if $\nu=1$, $d_{n}=n$, and $K=[a,b]$, then the
eigenvalues of $A_{n}$ are approximately equal to
$\psi\big{(}a+\frac{j}{n}(b-a)\big{)}$, $j=1,\ldots,n$, for $n$ large enough
and up to at most $o(n)$ outliers. Analogously, if we have $\nu=2$,
$d_{n}=n^{2}$, and $K=[a_{1},b_{1}]\times[a_{2},b_{2}]$, then the eigenvalues
of the matrix $A_{n}$ are approximately equal to
$\psi\big{(}a_{1}+\frac{j}{n}(b_{1}-a_{1}),a_{2}+\frac{k}{n}(b_{2}-a_{2})\big{)}$,
$j,k=1,\ldots,n$, for $n$ large enough and up to at most $o(n^{2})$ outliers.

The asymptotic distribution of eigenvalues and singular values of Toeplitz
matrix sequences has been studied deeply and continuously in the last century
(for example see [2, 3, 10, 16, 17] and references therein). The starting
point of this theory, which contains many extensions and other results, is a
famous theorem of Szegő [18], which we report in the version due to
Tyrtyshnikov and Zamarashkin [25].

######  Theorem 2.2.

If $f$ is integrable over $Q\equiv(-\pi,\pi)$, and if $\\{T_{n}(f)\\}$ is the
sequence of Toeplitz matrices generated by $f$, then

$\\{T_{n}(f)\\}\sim_{\sigma}(f,Q).$

Moreover, if $f$ is also real-valued, then each matrix $T_{n}(f)$ is Hermitian
and

$\\{T_{n}(f)\\}\sim_{\lambda}(f,Q).$

Furthermore, strong localization results are known in the case where the
generating function is real-valued, as stated in [20, Th.2.2] which we partly
report below.

######  Theorem 2.3.

If $f$ is integrable and real-valued almost everywhere over
$Q\equiv(-\pi,\pi)$, $m$ is the essential infimum of $f$, $M$ is the essential
supremum of $f$, and if $\\{T_{n}(f)\\}$ is the sequence of Toeplitz matrices
generated by $f$, then

$\lambda_{j}(T_{n}(f))\in(m,M)$

for every $j=1,\ldots,n$, for every positive integer $n$, under the assumption
that $m<M$. Moreover, if $m=M$ then the generating function $f$ is constant
almost everywhere and trivially we conclude that $T_{n}(f)$ coincides with $m$
times the identity matrix.

###  2.1 The Simple-Loop case

For $\alpha>0$, the well-known weighted Wiener algebra $W^{\alpha}$ is the
collection of all functions $f\colon\mathbb{T}\to\mathbb{C}$ whose Fourier
coefficients satisfy

$\|f\|_{\alpha}\equiv\sum_{j=-\infty}^{\infty}|\mathfrak{a}_{j}(f)|(|j|+1)^{\alpha}<\infty.$

It is easy to see that if $f\in W^{\alpha}$ then $f\in
C^{\lfloor\alpha\rfloor}[-\pi,\pi]$, hence the constant $\alpha$ is measuring
the smoothness of the symbol $f$. In what follows we extend every symbol $f$
to the whole real line in the natural way turning it into a $2\pi$-periodic
function, and we denote this extension by $f$ as well.

The simple-loop class, denoted by $\operatorname{SL}^{\alpha}$, consists in
the collection of all the real-valued symbols in $W^{\alpha}$ tracing out a
simple-loop over the interval $[-\pi,\pi]$ with the following properties:

  1. (i)

the range of $f$ is a segment $[0,\mu]$ with $\mu>0$;

  2. (ii)

$f(0)=f(2\pi)=0$, $f^{\prime\prime}(0)=f^{\prime\prime}(2\pi)>0$;

  3. (iii)

there is a unique $\theta_{0}\in[0,2\pi]$ such that $f(\theta_{0})=\mu$,
$f^{\prime}(\theta)>0$ for $\theta\in(0,\theta_{0})$ and
$f^{\prime}(\theta)<0$ for $\theta\in(\theta_{0},2\pi)$.

Note that if $f\in\operatorname{SL}^{\alpha}$ then
$f^{\prime}(0)=f^{\prime}(2\pi)=0$, and that $\theta_{0}=\pi$ for every even
symbol, i.e. a symbol satisfying $f(\theta)=f(-\theta)$ for
$\theta\in\mathbb{R}$.

Consider an even symbol $f\in\operatorname{SL}^{\alpha}$ with $\alpha>2$. For
a $n\times n$ matrix $A$ let $\lambda_{j}(A)$ $(j=1,\ldots,n)$ be its
eigenvalues. In the works [5, 6, 7, 8], for example, the authors state that

  1. (i)

the eigenvalues of $T_{n}(f)$ are all distinct, i.e.

$\lambda_{1}(T_{n}(f))<\lambda_{2}(T_{n}(f))<\cdots<\lambda_{n}(T_{n}(f));$

  2. (ii)

the numbers
$s_{j,n}\equiv\big{[}f\big{|}_{[0,\pi]}\big{]}^{-1}(\lambda_{j}(T_{n}(f)))$
$(j=1,\ldots,n)$ satisfy

$(n+1)s_{j,n}+\eta(s_{j,n})=\pi j+E_{j,n,\alpha},$ (2.1)

where $E_{j,n,\alpha}$ is an error term satisfying certain bounding condition,
and $\eta$ is a function depending only on $f$ with certain smoothness
depending on $\alpha$.

  3. (iii)

the previous equation (2.1) has exactly one solution $s_{j,n}\in[0,\pi]$ for
each $j=1,\ldots,n$.

In the bulk of the so-called simple-loop method, the authors use the Banach
fixed-point theorem to iterate over (2.1) and solve it for $s_{j,n}$,
obtaining an expansion of the kind (see [5, Th.2.2] for example)

$\lambda_{j}(T_{n}(f))=f(s_{j,n}),\qquad
s_{j,n}=\theta_{j,n}+\sum_{k=1}^{\lfloor\alpha\rfloor}r_{k}(\theta_{j,n})h^{k}+E_{j,n,\alpha},$
(2.2)

where

  * •

the numbers $s_{j,n}$ are arranged in nondecreasing order;

  * •

$h\equiv\frac{1}{n+1}$ and $\theta_{j,n}\equiv\pi jh$;

  * •

the coefficients $r_{k}$ depend only on $f$ and can be found explicitly, for
example

$r_{1}=-\eta$, |  $r_{3}=-\eta(\eta^{\prime})^{2}-\frac{1}{2}\eta^{2}\eta^{\prime\prime}$,
---|---
$r_{2}=\eta\eta^{\prime}$, |  $r_{4}=\eta(\eta^{\prime})^{3}+\frac{3}{2}\eta^{2}\eta^{\prime}\eta^{\prime\prime}+\frac{1}{6}\eta^{3}\eta^{\prime\prime\prime}$.

  * •

$E_{j,n,\alpha}=O(h^{\alpha})$ is the remainder (error) term, which satisfies
the bounding $|E_{j,n,\alpha}|\leqslant\kappa_{\alpha}h^{\alpha}$ for some
constant $\kappa_{\alpha}$ depending only on $\alpha$ and $f$.

Using the expansion (2.2) and the smoothness of the symbol $f$, the authors
apply $f$ to both sides obtaining an expansion of the kind

$\lambda_{j}(T_{n}(f))=f(\theta_{j,n})+\sum_{k=1}^{\lfloor\alpha\rfloor}c_{k}(\theta_{j,n})h^{k}+E_{j,n,\alpha},$
(2.3)

with similar characteristics but where the coefficients $c_{k}$ involve the
symbol and its derivatives, for instance $c_{1}=-f^{\prime}\eta$ and
$c_{2}=f^{\prime}\eta\eta^{\prime}+\frac{1}{2}f^{\prime\prime}\eta^{2}$. The
previous works [8, 13, 14, 15] used (2.3) as the basic expansion. We noticed
that (2.3) is absorbing all the almost non-increasining consequences of $f$ at
$\\{0,\pi\\}$ and all the related troubles that the derivatives of $f$ can
produce. Hence, we decided to work with (2.2) instead.

Indeed, in the light of Theorem 2.2, since $f$ is even and real-valued, we
find $\\{T_{n}(f)\\}\sim_{\lambda}(f,[0,\pi])$. Consequently we deduce

$\\{{\operatorname{diag}}_{j=1,\ldots,n}(s_{j,n})\\}\sim_{\lambda}({\rm
id},[0,\pi]),\ \ \ \ {\rm id}(\theta)\equiv\theta,$

with $s_{j,n}\in(0,\pi)$, by virtue of Theorem 2.3. Notice that the function
id$(\cdot)$ is very basic, and, as already claimed, when compared with the
study and the proposals in [5, 6, 9, 8, 12, 13, 14, 15], the troubles that the
derivatives of $f$ produce are completely removed: these nice features are
clearly evident by looking at the high precision of the numerical
computations, reported in Section 4, containing the numerical experiments.

##  3 The Algorithm

Our algorithm is based on the expansion (2.2) and is an evolution of the
algorithms proposed in [8, 13, 14, 15]. As in the mentioned works it is suited
for parallel implementation and can be called matrix-less since it does not
require to calculate or even, to store the matrix entries. For every
$n\in\mathbb{N}$ let $h\equiv\frac{1}{n+1}$ and $\theta_{j,n}\equiv\pi jh$,
thus the collection $\\{\theta_{j,n}\\}_{j=0}^{n+1}$ is a regular grid for the
interval $[0,\pi]$ with step size $\pi h$. We will use similar notations for
$n$ and $h$, for example, $h_{k}$ means $\frac{1}{n_{k}+1}$, and so on. We
assume that

  * •

the symbol $f$ is even and real-valued, strictly increasing in the interval
$[0,\pi]$, and $f(0)=0$;

  * •

$n_{1}$ and $\alpha$ are fixed natural numbers and $n\gg n_{1}$;

  * •

for $k=1,\ldots,\alpha$ let $n_{k}\equiv 2^{k-1}(n_{1}+1)-1$;

  * •

for $j_{1}=1,\ldots,n$ and $k=1,\ldots,\alpha$, let $j_{k}\equiv
2^{k-1}j_{1}$.

Note that $j_{k}$ depends on $j_{1}$, and similarly, $n_{k}$ depends on
$n_{1}$, but for notation simplicity, we suppressed those dependencies. The
indexes $j_{k}$ and the matrix sizes $n_{k}$ were calculated in such a way
that

$\theta_{j_{1},n_{1}}=\theta_{j_{2},n_{2}}=\cdots=\theta_{j_{\alpha},n_{\alpha}},$

which is the key idea of the following extrapolation phase, see Figure 1. We
also want to emphasize that $f$ is not necessarily simple-loop.

$\theta_{1,n_{1}}$$\theta_{2,n_{1}}$$\theta_{n_{1}-1,n_{1}}$$\theta_{n_{1},n_{1}}$$n_{1}$$n_{2}$$n_{\alpha-1}$$n_{\alpha}$
Figure 1: The regular grids $\\{\theta_{j,n_{k}}\\}$ for $j=1,\ldots,n_{k}$
and $k=1,\ldots,\alpha$. In the precomputing phase we will need to calculate
the eigenvalues of $T_{n_{k}}(f)$ for $k=1,\ldots,\alpha$ corresponding to the
blue and orange dots combined, but in the interpolation phase, we will use
only the eigenvalues corresponding to the orange dots, that is
$\lambda_{j_{k}}(T_{n_{k}}(f))$ for $j_{1}=1,\ldots,n_{1}$ and
$k=1,\ldots,\alpha$.

As in [8, 13, 14, 15], our algorithm is designed to calculate eigenvalues for
“big” matrix sizes $n$ with respect to $n_{1},\ldots,n_{\alpha}$, meaning
that, from a computational viewpoint, the calculation of the eigenvalues of
$T_{n}(f)$ is hard while for $T_{n_{k}}(f)$ can be easily done with any
standard eigensolver (i.e. Eigenvalues in Mathematica or eig in Matlab). But
our proposal is able to reach machine precision accuracy easily. The algorithm
has two phases, the first one involves an extrapolation procedure, and the
second one consists in a local interpolation technique. As a precomputing
phase we need to calculate the eigenvalues of $T_{n}(f)$ for
$n=n_{1},\ldots,n_{\alpha}$.

Extrapolation For each fixed $j_{1}=1,\ldots,n_{1}$ let
$\sigma_{j_{1}}\equiv\theta_{j_{1},n_{1}}=\cdots=\theta_{j_{\alpha},n_{\alpha}}$
(see the orange dots in Figure 1), and apply $\alpha$ times the expansion
(2.2) obtaining

$\displaystyle\lambda_{j_{1}}(T_{n_{1}}(f))-f(\sigma_{j_{1}})$
$\displaystyle=$ $\displaystyle
r_{1}(\sigma_{j_{1}})h_{1}+r_{2}(\sigma_{j_{1}})h_{1}^{2}+\cdots+r_{\alpha}(\sigma_{j_{1}})h_{1}^{\alpha}+E_{j_{1},n_{1},\alpha},$
$\displaystyle\lambda_{j_{2}}(T_{n_{2}}(f))-f(\sigma_{j_{1}})$
$\displaystyle=$ $\displaystyle
r_{1}(\sigma_{j_{1}})h_{2}+r_{2}(\sigma_{j_{1}})h_{2}^{2}+\cdots+r_{\alpha}(\sigma_{j_{1}})h_{2}^{\alpha}+E_{j_{2},n_{2},\alpha},$
$\displaystyle\vdots$
$\displaystyle\lambda_{j_{\alpha}}(T_{n_{\alpha}}(f))-f(\sigma_{j_{1}})$
$\displaystyle=$ $\displaystyle
r_{1}(\sigma_{j_{1}})h_{\alpha}+r_{2}(\sigma_{j_{1}})h_{\alpha}^{2}+\cdots+r_{\alpha}(\sigma_{j_{1}})h_{\alpha}^{\alpha}+E_{j_{\alpha},n_{\alpha},\alpha}.$

Let $\hat{r}_{k}(\sigma_{j_{1}})$ be the approximation of
$r_{k}(\sigma_{j_{1}})$ obtained by removing all the error terms
$E_{j_{k},n_{k},\alpha}$ and solving the resulting linear system:

$\begin{bmatrix}h_{1}&h_{1}^{2}&\cdots&h_{1}^{\alpha}\\\
h_{2}&h_{2}^{2}&\cdots&h_{2}^{\alpha}\\\ \vdots&\vdots&\ddots&\vdots\\\
h_{\alpha}&h_{\alpha}^{2}&\cdots&h_{\alpha}^{\alpha}\end{bmatrix}\begin{bmatrix}\hat{r}_{1}(\sigma_{j_{1}})\\\
\hat{r}_{2}(\sigma_{j_{1}})\\\ \vdots\\\
\hat{r}_{\alpha}(\sigma_{j_{1}})\end{bmatrix}=\begin{bmatrix}\lambda_{j_{1}}(T_{n_{1}}(f))\\\
\lambda_{j_{2}}(T_{n_{2}}(f))\\\ \vdots\\\
\lambda_{j_{\alpha}}(T_{n_{\alpha}}(f))\end{bmatrix}-f(\sigma_{j_{1}})\begin{bmatrix}1\\\
1\\\ \vdots\\\ 1\end{bmatrix}.$ (3.1)

As mentioned in previous works, a variant of this extrapolation strategy was
first suggested by Albrecht Böttcher in [5, §7] and is analogous to the
Richardson extrapolation employed in the context of Romberg integration [22,
§3.4].

Interpolation For any $j\in\\{1,\ldots,n\\}$ we will estimate
$r_{k}(\theta_{j,n})$. If $\theta_{j,n}$ coincides with one of the points in
the grid $\\{\theta_{0,n_{1}},\ldots,\theta_{n_{1}+1,n_{1}}\\}$, then we have
the approximations $\hat{r}(\theta_{j,n})$ from the extrapolation phase for
free. In any other case, we will do it by interpolating the data

$(\theta_{0,n_{k}},\hat{r}_{k}(\theta_{0,n_{k}})),(\theta_{1,n_{k}},\hat{r}_{k}(\theta_{1,n_{k}})),\ldots,(\theta_{n_{k}+1,n_{k}},\hat{r}_{k}(\theta_{n_{k}+1,n_{k}})),$

for $k=1,\ldots,\alpha$, and then evaluating the resulting polynomial at
$\theta_{j,n}$. This interpolation can be done in many ways, but to avoid
spurious oscillations explained by the Runge phenomenon [11, p.78], and
following the strategy of the previous works, we decided to do it considering
only the $\alpha-k+5$ points in the grid
$\\{\theta_{0,n_{1}},\ldots,\theta_{n_{1}+1,n_{1}}\\}$ which are closest to
$\theta_{j,n}$. Those points can be determined uniquely unless $\theta_{j,n}$
is the mid point of two consecutive points in the grid, in which case we can
take any of the possible two choices.

Finally, our eigenvalue approximation with $k$ terms, is given by

$\hat{\lambda}_{j,k}^{\textrm{NAS}}(T_{n}(f))\equiv
f\Big{(}\theta_{j,n}+\sum_{\ell=0}^{k-1}\hat{r}_{\ell}(\theta_{j,n})h^{\ell}\Big{)},$
(3.2)

where $k=1,\ldots,\alpha$, and NAS stands for “Numerical Algorithm in the
variable $s_{j,k}$”.

######  Remark 3.1.

To get the best results possible, in the precomputing phase, we advise to use
the proposed algorithm (3.2), calculating the eigenvalues of $T_{n}(f)$ for
$n=n_{1},\ldots,n_{\alpha}$ with a significant number of precision digits,
let’s say $60$. For $\alpha=5$ and $n_{1}=100$, for instance, this can be done
in a standard computer in a few minutes, and it only needs to be done once.
While for the extrapolation phase (3.1), we advice to do ti with the largest
$\alpha$ possible.

##  4 Numerical Experiments

Let $\lambda_{j,k}^{\operatorname{SL}}(T_{n}(f))$ be the $k$th term
approximation of $\lambda_{j}(T_{n}(f))$ obtained with the simple-loop method
in the variable $s_{j,n}$, that is

$\lambda_{j,k}^{\operatorname{SL}}(T_{n}(f))\equiv
f\Big{(}\theta_{j,n}+\sum_{k=1}^{k-1}r_{k}(\theta_{j,n})h^{k}\Big{)},$

where $\theta_{j,n}=\pi jh$. Let also $\lambda_{j,k}^{\textrm{NA}}(T_{n}(f))$
be the $k$th term approximation of $\lambda_{j}(T_{n}(f))$ given by the
Numerical Algorithm in [14] and finally, let
$\lambda_{j,k}^{\textrm{MNA}}(T_{n}(f))$ be the respective approximation given
by the Modified Numerical Algorithm [8, §4]. In order to compare the results
of the different methods we use the following notation for the absolute
individual errors

$\displaystyle\varepsilon^{\operatorname{SL}}_{j,n,k}$
$\displaystyle\equiv|\lambda_{j}(T_{n}(f))-\lambda_{j,k}^{\operatorname{SL}}(T_{n}(f))|,$
$\displaystyle\varepsilon^{\textrm{NA}}_{j,n,k}$
$\displaystyle\equiv|\lambda_{j}(T_{n}(f))-\lambda_{j,k}^{\textrm{NA}}(T_{n}(f))|,$
$\displaystyle\varepsilon^{\textrm{MNA}}_{j,n,k}$
$\displaystyle\equiv|\lambda_{j}(T_{n}(f))-\lambda_{j,k}^{\textrm{MNA}}(T_{n}(f))|,$
$\displaystyle\varepsilon^{\textrm{NAS}}_{j,n,k}$
$\displaystyle\equiv|\lambda_{j}(T_{n}(f))-\lambda_{j,k}^{\textrm{NAS}}(T_{n}(f))|,$

and the respective maximum absolute errors

$\displaystyle\varepsilon^{\operatorname{SL}}_{n,k}$
$\displaystyle\equiv\max\\{\varepsilon^{\operatorname{SL}}_{j,n,k}\colon
j=1,\ldots,n\\},$ $\displaystyle\varepsilon^{\textrm{NA}}_{n,k}$
$\displaystyle\equiv\max\\{\varepsilon^{\textrm{NA}}_{j,n,k}\colon
j=1,\ldots,n\\}],$ $\displaystyle\varepsilon^{\textrm{MNA}}_{n,k}$
$\displaystyle\equiv\max\\{\varepsilon^{\textrm{MNA}}_{j,n,k}\colon
j=1,\ldots,n\\},$ $\displaystyle\varepsilon^{\textrm{NAS}}_{n,k}$
$\displaystyle\equiv\max\\{\varepsilon^{\textrm{NAS}}_{j,n,k}\colon
j=1,\ldots,n\\}.$

We start with an example involving a well-known simple-loop symbol for which
we can exactly calculate the coefficients $r_{k}$ in (2.2) easily, thus we
will be able to compare the accuracy of the different eigenvalue
approximations.

######  Example 4.1 (A simple-loop symbol).

Consider the even simple-loop symbol given by

$f(\theta)\coloneqq\frac{(1+\rho)^{2}}{2}\cdot\frac{1-\cos(\theta)}{1-2\rho\cos(\theta)+\rho^{2}}\qquad(0\leqslant\theta\leqslant
2\pi),$ (4.1)

for a constant $0<\rho<1$, see Figure 2. This symbol was inspired in the
Kac–Murdock–Szegő Toeplitz matrices introduced in [19] and subsequently
studied in [23, 24], which usually are present in important physics models.
The respective Fourier coefficients can be explicitly calculated as
$\mathfrak{a}_{k}(f)=\frac{1}{4}(\rho^{2}-1)\rho^{|k|-1}$ for $k\neq 0$ and
$\frac{1}{2}(1+\rho)$ for $k=0$.

$f$$f^{\prime}$
$f^{\prime\prime}$

Figure 2: The symbol $f$ in (4.1) and its first two derivatives for
$\rho=\frac{1}{2}$.

We then have

$\|f\|_{\alpha}=\frac{1+\rho}{2}+\frac{\rho^{2}-1}{2\rho}\sum_{k=1}^{\infty}\rho^{k}(k+1)^{\alpha},$

which is finite for every $\alpha>0$, the remaining simple-loop conditions are
easily verified in Figure 2. Then $f\in\operatorname{SL}^{\alpha}$ for any
$\alpha>0$. According to [8, §4], the function $\eta$ in (2.1) is nicely given
by

$\eta(s)=2\arctan\Big{(}\frac{\rho\sin(s)}{1-\rho\cos(s)}\Big{)}.$

The Figure 3 and the Table 1 shows that, in this case, the approximation
$\lambda_{j,k}^{\textrm{MNA}}(T_{n}(f))$ can produce good results until the
level $k=2$ but it becomes unstable from this point on. While our algorithm
(3.2) is still producing fine results in the $4$th level and, for a matrix of
size $4096$, it reaches machine-precision from level $3$. The Figure 3 reveals
also, that for the symbol (4.1), our proposed algorithm (3.2) can match the
exact asymptotic simple-loop expansion until level $4$.

$n$ | $256$ | $512$ | $1024$ | $2048$ | $4096$
---|---|---|---|---|---
$\varepsilon_{n,1}^{\textrm{MNA}}$ | $3.0897\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $1.5494\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $7.7577\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $3.8816\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $1.9415\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$
$\varepsilon_{n,1}^{\textrm{NAS}}$ | $3.0897\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $1.5494\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $7.7577\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $3.8816\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $1.9415\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$
$(n+1)\varepsilon_{n,1}^{\textrm{NAS}}$ | $7.9405\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $7.9482\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $7.9517\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $7.9534\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $7.9542\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$
$\varepsilon_{n,2}^{\textrm{MNA}}$ | $1.5098\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $3.7949\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $9.5024\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $2.3819\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $5.9529\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$
$\varepsilon_{n,2}^{\textrm{NAS}}$ | $1.3575\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $3.4113\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $8.5515\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $2.1407\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $5.3553\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$
$(n+1)^{2}\varepsilon_{n,2}^{\textrm{NAS}}$ | $8.9661\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $8.9775\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $8.9844\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $8.9875\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $8.9890\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$
$\varepsilon_{n,3}^{\textrm{MNA}}$ | $8.8167\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $1.4780\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $4.7324\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $2.4238\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.2270\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$
$\varepsilon_{n,3}^{\textrm{NAS}}$ | $5.4356\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $6.8619\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $8.6153\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $1.0794\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $1.3507\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$
$(n+1)^{3}\varepsilon_{n,3}^{\textrm{NAS}}$ | $9.2267\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $9.2640\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $9.2778\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $9.2852\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $9.2887\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$
$\varepsilon_{n,4}^{\textrm{MNA}}$ | $6.6888\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $3.1948\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.5778\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $7.8461\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $3.8983\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$
$\varepsilon_{n,4}^{\textrm{NAS}}$ | $3.4700\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $2.1887\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$ | $1.3740\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-12}$ | $8.6077\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-14}$ | $5.4131\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-15}$
$(n+1)^{4}\varepsilon_{n,4}^{\textrm{NAS}}$ | $1.5138\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5158\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5166\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5172\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5252\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$
Table 1: The maximum errors $\varepsilon_{n,k}^{\textrm{MNA}}$,
$\varepsilon_{n,k}^{\textrm{NAS}}$, and maximum normalized errors
$(n+1)^{k}\varepsilon_{n,k}^{\textrm{NAS}}$ for the levels $k=1,2,3,4$ and
different matrix sizes $n$, corresponding to the symbol (4.1) with
$\rho=\frac{1}{2}$. We used a grid of size $n_{1}=100$.

$k=2$
$k=3$
$k=4$

Figure 3: The base-10 logarithm for the individual errors
$\varepsilon_{j,n,k}^{\operatorname{SL}}$ (thick gray),
$\varepsilon_{j,n,k}^{\textrm{MNA}}$ (green), and
$\varepsilon_{j,n,k}^{\textrm{NAS}}$ (blue) for the symbol $f$ in (4.1) with
$\rho=\frac{1}{2}$, a matrix size $n=4096$, a grid size $n_{1}=100$, and
different levels $k$. The proposed numerical algorithm (3.2) is represented by
the blue curve, while the exact asymptotic simple-loop expansion, by the gray
curve.

######  Example 4.2 (A non-simple-loop symbol).

We now test our algorithm with a Real Cosine Trigonometric Polynomial (RCTP),
see [15, §1]. For $\ell\in\mathbb{Z}_{+}$, consider the symbol

$f_{\ell}(\theta)\equiv(2-2\cos(\theta))^{\ell},\quad\theta\in Q.$ (4.2)

The respective Fourier coefficients can be exactly calculated as
$\mathfrak{a}_{k}(f_{\ell})=(-1)^{k}{{2\ell}\choose{\ell+k}}$ for
$|k|\leqslant\ell$ and $\mathfrak{a}_{k}(f_{\ell})=0$ in any other case, then
the respective Toeplitz matrices $T_{n}(f_{\ell})$ are banded with a band of
size $2\ell+1$. The case $\ell=2$ was carefully studied by Barrera and Grudsky
in [4], where they formally deduced that

$s_{j,n}=\frac{\pi(j+1)}{n+2}+\frac{u_{1,j}}{n+2}+\frac{u_{2,j}}{(n+2)^{2}}+O(h^{3}),$

with some bounded and continuous coefficients $u_{1,j},u_{2,j}$. See Theorem
2.5 there. The previous expansion is slightly different from (2.2) but we was
able to show that our algorithm is producing fine results in this case.

It is clear that $f_{\ell}\in W^{\alpha}$ for any $\ell\in\mathbb{Z}_{+}$ and
any $\alpha>0$, but $f_{\ell}$ is simple-loop only when $\ell=1$ because
$f^{\prime\prime}_{\ell}(0)=0$ for $\ell\neq 1$. The Figures 5, 6, and the
Tables 2, 3, show the data for the cases $\ell=2,3$. Since our method is based
on the simple-loop expansion (2.2), we expected difficulties for the very
first eigenvalues, corresponding to the point $\theta=0$, nevertheless, the
numerical approximations are good enough for machine precision purposes.

The Figure 4 shows a comparison between the individual errors
$\varepsilon_{j,n,k}^{\textrm{NA}}$, given by the eigenvalue approximation of
the numerical algorithm [15], and $\varepsilon_{j,n,k}^{\textrm{MNA}}$, which
corresponds to its boundary modification given by [8, §4]. Then it is clear
that the modified version works better.

Figure 4: The base-10 logarithm for the individual errors $\varepsilon_{j,n,k}^{\textrm{NA}}$ (thick cyan curve) and $\varepsilon_{j,n,k}^{\textrm{MNA}}$ (thin blue curve), corresponding to the algorithms [15] and [8, §4], respectively. We worked with the symbol (4.2) where $\ell=3$, level $k=4$, a matrix size $n=4096$, and a grid size $n_{1}=100$. $n$ | $256$ | $512$ | $1024$ | $2048$ | $4096$
---|---|---|---|---|---
$\varepsilon_{n,1}^{\textrm{MNA}}$ | $1.6269\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $8.1578\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $4.0848\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $2.0439\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $1.0223\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$
$\varepsilon_{n,1}^{\textrm{NAS}}$ | $1.6269\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $8.1578\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $4.0848\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $2.0439\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $1.0223\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$
$(n+1)\varepsilon_{n,1}^{\textrm{NAS}}$ | $4.1811\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $4.1850\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $4.1869\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $4.1878\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $4.1883\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$
$\varepsilon_{n,2}^{\textrm{MNA}}$ | $2.2934\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $5.7602\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $1.4434\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $3.6126\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $9.0367\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$
$\varepsilon_{n,2}^{\textrm{NAS}}$ | $2.7270\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $6.8421\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $1.7136\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $4.2880\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $1.0725\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$
$(n+1)^{2}\varepsilon_{n,2}^{\textrm{NAS}}$ | $1.8011\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.8006\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.8004\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.8003\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.8002\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$
$\varepsilon_{n,3}^{\textrm{MNA}}$ | $7.5158\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $9.4903\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.2050\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.5758\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $2.5206\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$
$\varepsilon_{n,3}^{\textrm{NAS}}$ | $6.9024\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $8.6696\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.0863\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.3595\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $1.7004\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$
$(n+1)^{3}\varepsilon_{n,3}^{\textrm{NAS}}$ | $1.1717\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.1704\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.1698\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.1695\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.1694\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$
$\varepsilon_{n,4}^{\textrm{MNA}}$ | $3.7838\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $2.2540\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $1.6954\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$ | $5.9426\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-12}$ | $3.2321\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-12}$
$\varepsilon_{n,4}^{\textrm{NAS}}$ | $2.7800\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.3631\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $7.4328\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-12}$ | $4.5503\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-13}$ | $5.4968\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-14}$
$(n+1)^{4}\varepsilon_{n,4}^{\textrm{NAS}}$ | $1.2128\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $9.4408\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $8.2044\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $8.2044\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5487\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$
Table 2: The maximum errors $\varepsilon_{n,k}^{\textrm{MNA}}$, $\varepsilon_{n,k}^{\textrm{NAS}}$, and maximum normalized errors $(n+1)^{k}\varepsilon_{n,k}^{\textrm{NAS}}$ for the levels $k=1,2,3,4$ and different matrix sizes $n$, corresponding to the symbol (4.2) with $\ell=2$. We used a grid of size $n_{1}=100$. $n$ | $256$ | $512$ | $1024$ | $2048$ | $4096$
---|---|---|---|---|---
$\varepsilon_{n,1}^{\textrm{MNA}}$ | $9.1868\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $4.6172\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $2.3146\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $1.1588\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $5.7978\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$
$\varepsilon_{n,1}^{\textrm{NAS}}$ | $9.1868\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $4.6172\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $2.3146\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $1.1588\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-2}$ | $5.7978\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$
$(n+1)\varepsilon_{n,1}^{\textrm{NAS}}$ | $2.3610\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.3686\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.3725\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.3744\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.3753\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$
$\varepsilon_{n,2}^{\textrm{MNA}}$ | $2.2588\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $5.6759\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $1.4227\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $3.5615\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $8.9091\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$
$\varepsilon_{n,2}^{\textrm{NAS}}$ | $3.0497\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $7.6550\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $1.9176\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $4.7989\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $1.2003\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$
$(n+1)^{2}\varepsilon_{n,2}^{\textrm{NAS}}$ | $2.0143\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.0146\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.0147\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.0148\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.0148\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$
$\varepsilon_{n,3}^{\textrm{MNA}}$ | $9.5995\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $1.2158\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $1.5664\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $2.1683\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $3.7673\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$
$\varepsilon_{n,3}^{\textrm{NAS}}$ | $1.3355\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $1.6765\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $2.1002\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $2.6281\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $3.2868\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$
$(n+1)^{3}\varepsilon_{n,3}^{\textrm{NAS}}$ | $2.2669\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.2634\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.2617\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.2608\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $2.2604\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$
$\varepsilon_{n,4}^{\textrm{MNA}}$ | $5.7200\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $3.9578\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $1.7663\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $9.3710\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$ | $4.8060\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$
$\varepsilon_{n,4}^{\textrm{NAS}}$ | $7.6467\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $4.8020\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $3.0083\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$ | $1.8824\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-12}$ | $1.1772\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-13}$
$(n+1)^{4}\varepsilon_{n,4}^{\textrm{NAS}}$ | $3.3358\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $3.3258\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $3.3206\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $3.3181\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $3.3168\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$
Table 3: The same as Table 2 but this time with $\ell=3$.

$k=2$
$k=3$
$k=4$

Figure 5: The base-10 logarithm for the individual errors
$\varepsilon_{j,n,k}^{\operatorname{SL}}$ (thick gray),
$\varepsilon_{j,n,k}^{\textrm{MNA}}$ (green), and
$\varepsilon_{j,n,k}^{\textrm{NAS}}$ (blue) for the RCTP symbol $f_{\ell}$ in
(4.2) with $\ell=2$, a matrix size $n=4096$, a grid size $n_{1}=100$, and
different levels $k$.

$k=2$
$k=3$
$k=4$

Figure 6: The same as Figure 5 but this time with $\ell=3$.

######  Example 4.3 (A matrix order dependent symbol).

We now test our algorithm with a symbol which is a linear combination of RCTPs
with coefficients depending on the matrix order $n$. Consider the symbol

$F_{n}(\theta)\equiv
f_{2}(\theta)+\alpha_{1}f_{1}(\theta)h^{2}+\alpha_{0}f_{0}(\theta)h^{4},$
(4.3)

where $f_{\ell}$ is given by (4.2) and $\alpha_{\ell}$ are real constants. We
previously studied this symbol in [8, §4.1] where we proposed an improvement
to the numerical algorithm [14]. This symbol commonly arises when discretizing
differential equations with the Finite Differences method.

The respective Fourier coefficients can be exactly calculated using the
previous example and the linearity of the Fourier transform, as
$\mathfrak{a}_{k}(F_{n})=\alpha_{2}\mathfrak{a}_{k}(f_{2})+\alpha_{1}\mathfrak{a}_{k}(f_{1})h^{2}+\alpha_{0}\mathfrak{a}_{k}(f_{0})h^{4}$.
Therefore, the Toeplitz matrices $T_{n}(F_{n})$ are banded and penta-diagonal.
For implementing the numerical algorithm in [14], we need to assume an
eigenvalue expansion with the form

$\lambda_{j}(T_{n}(F_{n}))=f_{2}(\theta_{j,n})+\sum_{\ell=1}^{\alpha-1}c_{\ell}(\theta_{j,n})h^{\ell}+E_{\alpha,j,n},$

where the coefficients $c_{\ell}$ are continuous and bounded functions from
$[0,\pi]$ to $\mathbb{R}$, and the remainder (error) term $E_{\alpha,j,n}$
satisfies the inequality $|E_{\alpha,j,n}|\leqslant\kappa_{\alpha}h^{\alpha}$.
For the respective boundary modification proposed in [8], we need to note that

$\displaystyle c_{2}(0)$ $\displaystyle=\alpha_{1}f_{1}(0)=0,$ $\displaystyle
c_{2}(\pi)$ $\displaystyle=\alpha_{1}f_{1}(\pi)=4\alpha_{1},$ $\displaystyle
c_{4}(0)$ $\displaystyle=\alpha_{0}f_{0}(0)=0,$ $\displaystyle c_{4}(\pi)$
$\displaystyle=\alpha_{0}f_{0}(\pi)=\alpha_{0},$

while $c_{\ell}(0)=c_{\ell}(\pi)=0$ in any other case. The Figure 7 and the
Table 4, show the data.

$n$ | $256$ | $512$ | $1024$ | $2048$ | $4096$
---|---|---|---|---|---
$\varepsilon_{n,1}^{\textrm{MNA}}$ | $6.0007\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $3.0088\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $1.5065\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $7.5377\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $3.7702\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$
$\varepsilon_{n,1}^{\textrm{NAS}}$ | $6.0007\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $3.0088\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $1.5065\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-3}$ | $7.5377\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$ | $3.7702\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-4}$
$(n+1)\varepsilon_{n,1}^{\textrm{NAS}}$ | $1.5422\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5435\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5442\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5445\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5446\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$
$\varepsilon_{n,2}^{\textrm{MNA}}$ | $3.9869\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $1.0043\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $2.5287\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $6.3851\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $1.6273\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$
$\varepsilon_{n,2}^{\textrm{NAS}}$ | $1.5208\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-5}$ | $3.7944\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-6}$ | $9.4766\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $2.3679\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-7}$ | $5.9184\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$
$(n+1)^{2}\varepsilon_{n,2}^{\textrm{NAS}}$ | $1.0045\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $9.9858\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $9.9563\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $9.9416\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$ | $9.9342\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-1}$
$\varepsilon_{n,3}^{\textrm{MNA}}$ | $9.2045\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $1.0710\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $3.6629\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.8241\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $9.1977\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$
$\varepsilon_{n,3}^{\textrm{NAS}}$ | $8.9731\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $1.1313\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-8}$ | $1.4203\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $1.7792\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $2.2264\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$
$(n+1)^{3}\varepsilon_{n,3}^{\textrm{NAS}}$ | $1.5231\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5274\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5295\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5306\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$ | $1.5311\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{\,0}$
$\varepsilon_{n,4}^{\textrm{MNA}}$ | $3.4333\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $2.1389\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $5.3710\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$ | $2.5589\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$ | $1.2650\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$
$\varepsilon_{n,4}^{\textrm{NAS}}$ | $4.3281\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-9}$ | $2.7008\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-10}$ | $1.8110\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-11}$ | $2.3324\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-12}$ | $2.9853\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{-13}$
$(n+1)^{4}\varepsilon_{n,4}^{\textrm{NAS}}$ | $1.8881\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $1.8705\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $1.9990\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $4.1112\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$ | $8.4112\raisebox{0.56905pt}{\mbox{\tiny$\times$}}10^{1}$
Table 4: The maximum errors $\varepsilon_{n,k}^{\textrm{MNA}}$,
$\varepsilon_{n,k}^{\textrm{NAS}}$, and maximum normalized errors
$(n+1)^{k}\varepsilon_{n,k}^{\textrm{NAS}}$ for the levels $k=1,2,3,4$ and
different matrix sizes $n$, corresponding to the symbol (4.3) with
$\alpha_{0}=3$ and $\alpha_{1}=2$. We used a grid of size $n_{1}=100$.

$k=2$
$k=3$
$k=4$

Figure 7: The base-10 logarithm for the individual errors
$\varepsilon_{j,n,k}^{\operatorname{SL}}$ (thick gray),
$\varepsilon_{j,n,k}^{\textrm{MNA}}$ (green), and
$\varepsilon_{j,n,k}^{\textrm{NAS}}$ (blue) for the matrix order dependent
symbol $F_{n}$ in (4.3), with $\alpha_{0}=3$, $\alpha_{1}=2$, a matrix size
$n=4096$, a grid size $n_{1}=100$, and different levels $k$.

##  5 Conclusions

Under appropriate technical assumptions, the simple-loop theory allows to
deduce various types of asymptotic expansions for the eigenvalues of Toeplitz
matrices $T_{n}(f)$ generated by a function $f$. Independently and under the
milder hypothesis that $f$ is even and monotonic over $[0,\pi]$, matrix-less
algorithms have been developed for the fast eigenvalue computation of large
Toeplitz matrices. These procedures work with a linear complexity in the
matrix order $n$ and behind the high efficiency of such algorithms there are
the expansions predicted by the simple-loop theory, combined with the
extrapolation idea.

In this note we have focused our attention on a change of variable, followed
by the asymptotic expansion of the new variable that is

$\lambda_{j}(T_{n}(f))\equiv f(s_{j,n}),\ \ \ \ \
s_{j,n}=\theta_{j,n}+\sum_{k=1}^{\lfloor\alpha\rfloor}r_{k}(\theta_{j,n})h^{k}+E_{j,n,\alpha},$

and then we adapted the matrix-less to the considered new setting.

Numerical experiments have shown in a clear way a higher precision (till
machine precision) and the same linear computation cost, when compared with
the matrix-less procedures already presented in the relevant literature. More
specifically, among the advantages, we concisely mention the following:

  1. a)

when the coefficients of the simple-loop function are analytically known, the
algorithm computes them perfectly;

  2. b)

while the proposed algorithm is better or at worst comparable to the previous
ones for the computation of the inner eigenvalues, it is extremely better for
the computation of the extreme eigenvalues, which are essential for
determining important quantities, like the conditioning in the positive
definite case.

As next steps the following questions remain to be investigated:

  * •

a fine error analysis for having a theoretical explanation of the reason why
the new expansion leads to a much smaller errors, when compared with the
numerical results in [12, 13, 15, 14];

  * •

applications to the block cases (see [2, 3] for the theory in the block case)
and related applications [12, 13] to differential problems;

  * •

taking inspiration from [1], extension of the technique to preconditioned case
$X_{n}=T_{n}^{-1}(g)T_{n}(l)$ in which $g$ is positive over $(0,\pi)$ and is
not identically constant, making use of the ergodic theorems given in [21],
where Theorem 2.2 is extended to the preconditioned case in a very general
setting.

## References

  * [1] Ahmad, F., Al-Aidarous, E.S., Alrehaili, D.A., Ekström, S.E., Furci, I., and Serra-Capizzano, S. Are the eigenvalues of preconditioned banded symmetric Toeplitz matrices known in almost closed form?  Numer. Algorithms 78, 3 (2018), 867–893.
  * [2] Barbarino, G., Garoni, C., and Serra-Capizzano, S. Block generalized locally Toeplitz sequences: Theory and applications in the multidimensional case.  Electron. Trans. Numer. Anal. 53 (2020), 113–216.
  * [3] Barbarino, G., Garoni, C., and Serra-Capizzano, S. Block generalized locally Toeplitz sequences: Theory and applications in the unidimensional case.  Electron. Trans. Numer. Anal. 53 (2020), 28–112.
  * [4] Barrera, M., and Grudsky, S.M. Asymptotics of eigenvalues for pentadiagonal symmetric Toeplitz matrices.  Oper. Theory: Adv. Appl. 259 (2017), 51–77.
  * [5] Bogoya, M., Böttcher, A., Grudsky, S.M., and Maximenko, E.A. Eigenvalues of Hermitian Toeplitz matrices with smooth simple-loop symbols.  Oper. Theory: Adv. Appl. 422 (2015), 1308–1334.
  * [6] Bogoya, M., Böttcher, A., Grudsky, S.M., and Maximenko, E.A. Eigenvectors of Hermitian Toeplitz matrices with smooth simple-loop symbols.  Linear Algebra Appl. 493 (2016), 606–637.
  * [7] Bogoya, M., Grudsky, S.M., and Maximenko, E.A. Eigenvalues of Hermitian Toeplitz matrices generated by simple-loop symbols with relaxed smoothness.  Oper. Theory: Adv. Appl. 259 (2017), 179–212.
  * [8] Bogoya, M., and Serra-Capizzano, S. Eigenvalue superposition expansion for toeplitz matrix-sequences, generated by linear combinations of matrix-order dependent symbols, and applications to fast eigenvalue computations.  arXiv:2112.11794 (2022).
  * [9] Böttcher, A., Bogoya, M., Grudsky, S.M., and Maksimenko, r. Asymptotics of the eigenvalues and eigenvectors of Toeplitz matrices.  Mat. Sb. 208, 11 (2017), 4–28.
  * [10] Böttcher, A., and Silbermann, B. Introduction to large truncated Toeplitz matrices.  Universitext. Springer-Verlag, New York, 1999.
  * [11] Davis, P.J. Interpolation and approximation.  Dover, New York, 1975.
  * [12] Ekström, S.E., Furci, I., Garoni, C., Manni, C., Serra-Capizzano, S., and Speleers, H. Are the eigenvalues of the B-spline isogeometric analysis approximation of $-\delta u=\lambda u$ known in almost closed form?  Numer. Linear Algebra Appl. 25, 5 (2018), e2198, 34 pp.
  * [13] Ekström, S.E., Furci, I., and Serra-Capizzano, S. Exact formulae and matrix-less eigensolvers for block banded Toeplitz-like matrices.  BIT Numerical Mathematics 58, 4 (2018), 937–968.
  * [14] Ekström, S.E., and Garoni, C. A matrix-less and parallel interpolation-extrapolation algorithm for computing the eigenvalues of preconditioned banded symmetric Toeplitz matrices.  Numer. Algor. 80 (2019), 819–848.
  * [15] Ekström, S.E., Garoni, C., and Serra-Capizzano, S. Are the eigenvalues of banded symmetric Toeplitz matrices known in almost closed form?  Exper. Math. 27, 4 (2018), 478–487.
  * [16] Garoni, C., and Serra-Capizzano, S. Generalized Locally Toeplitz sequences: Theory and applications. Vol. I.  Springer, Cham, 2017.
  * [17] Garoni, C., and Serra-Capizzano, S. Generalized Locally Toeplitz sequences: Theory and applications. Vol. II.  Springer, Cham, 2018.
  * [18] Grenander, U., and Szegő, G. Toeplitz forms and their applications, second ed.  California Monographs in Mathematical Sciences. Chelsea Publishing Co., New York, 1984.
  * [19] Kac, M., Murdock, W.L., and Szegő, G. On the eigenvalues of certain Hermitian forms.  J. Rational Mech. Anal. 2 (1953), 767–800.
  * [20] Serra-Capizzano, S. The extension of the concept of the generating function to a class of preconditioned Toeplitz matrices.  Linear Algebra Appl. 267 (1997), 139–161.
  * [21] Serra-Capizzano, S. An ergodic theorem for classes of preconditioned matrices.  Linear Algebra Appl. 282, 1-3 (1998), 161–183.
  * [22] Stoer, J., and Bulirsch, R. Introduction to numerical analysis, third ed.  Springer, 2010.
  * [23] Trench, W.F. Asymptotic distribution of the spectra of a class of generalized Kac–Murdock–Szegő matrices.  Linear Algebra Appl. 294 (1999), 181–192.
  * [24] Trench, W.F. Spectral decomposition of Kac–Murdock–Szegő matrices.  The selected works of William F. Trench. http://works.bepress.com/william_trench/133, 2010.
  * [25] Tyrtyshnikov, E.E., and Zamarashkin, N.L. Spectra of multilevel Toeplitz matrices: advanced theory via simple matrix relationships.  Linear Algebra Appl. 270 (1998), 15–27.