part. Without the support from the picture of the pseudospectrum, the finite
section method does not provide information regarding the boundary of the
essential numerical range of $T$ \- there is a misleading circle of
eigenvalues of $P_{m}T|_{P_{m}\mathcal{H}}$ which do not occur along the
boundary of the essential spectrum but are simply given by the diagonal
entries $\\{t_{5},t_{6},...,t_{m}\\}$.

  * Remark 6.3

The previous examples demonstrated that, in general, the finite section method
is not always suitable for computing spectra. Rather then working with square
sections of the infinite matrix $T$, one should work with uneven sections
$P_{n}TP_{m}$, where the parameters $n$ and $m$ are allowed to vary
independently. Indeed, the algorithms presented in [28, 42] use this method.
In effect, we need to know how large $n$ should be to retain enough
information of the operator $TP_{m}$. This type of idea is also used
implicitly in the IQR algorithm (see Section 4).

###  6.2 Numerical examples I: normal operators

  * Example 6.4 (Convergence of the IQR algorithm)

We begin with two simple examples that demonstrate the linear (or exponential)
convergence proven in Theorem 3.9 and Corollary 3.12 (and its
generalisations). Consider first the one-dimensional discrete Schrödinger
operator given by

$T_{1}=\left(\begin{matrix}v_{1}&1&&&\\\ 1&v_{2}&1&&\\\ &1&v_{3}&1&\\\
&&1&v_{4}&\ddots\\\ &&&\ddots&\ddots\\\ \end{matrix}\right),$

where $v_{j}=5\sin(j)^{2}/\sqrt{j}$ if $j\leq 10$ and $v_{j}=0$ otherwise. As
a compact (in fact finite rank) perturbation of the free Laplacian,
$\sigma(T_{1})$ consists of the interval $[-2,2]$ together with isolated
eigenvalues of finite multiplicity which can be computed [77]. The second
operator, $T_{2}$, consists of taking the operator

$T_{0}=\left(\begin{matrix}2&0&0&0\\\ 0&\frac{3i}{2}&0&0\\\
0&0&-\frac{5}{4}&0\\\ 0&0&0&-\frac{9i}{8}\end{matrix}\right)\bigoplus U_{1},$

where $U_{k}$ denotes the bilateral shift $e_{j}\rightarrow e_{j+k}$, writing
this as an operator on $l^{2}(\mathbb{N})$ and then mixing the spaces via a
random unitary transformation on the span of the first $9$ basis vectors. This
ensures $T_{2}$ is not written in block form but has known eigenvalues. We
have plotted the difference in norm between the first $j\times j$ block of
each $Q_{n}^{*}T_{l}Q_{n}$ and the diagonal operator formed via the largest
$j$ eigenvalues for $j=1,2,3$ and $4$ in Fig. 2. The plot clearly shows the
exponential convergence.

Figure 2: Exponential convergence to the diagonal blocks for $T_{1}$ and
$T_{2}$.

  * Example 6.5 (Convergence to extremal parts of the spectrum)

To see why we may need some condition on $\sigma(T)$ for convergence of the
IQR algorithm to the extreme parts of the spectrum, we consider Laurent and
Toeplitz operators with symbol given by a trigonometric polynomial

$a(t)=\sum_{j=-k}^{j=k}a_{j}t^{j}.$

Given such a symbol, we define Laurent and Toeplitz operators

$L(a)=\left(\begin{array}[]{ccc|cccc}\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\\
\cdots&a_{0}&a_{-1}&a_{-2}&a_{-3}&a_{-4}&\cdots\\\
\cdots&a_{1}&a_{0}&a_{-1}&a_{-2}&a_{-3}&\cdots\\\
\hline\cr\cdots&a_{2}&a_{1}&a_{0}&a_{-1}&a_{-2}&\cdots\\\
\cdots&a_{3}&a_{2}&a_{1}&a_{0}&a_{-1}&\cdots\\\
\cdots&a_{4}&a_{3}&a_{2}&a_{1}&a_{0}&\cdots\\\
\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\end{array}\right),\quad
T(a)=\left(\begin{array}[]{cccc}a_{0}&a_{-1}&a_{-2}&\cdots\\\
a_{1}&a_{0}&a_{-1}&\cdots\\\ a_{2}&a_{1}&a_{0}&\cdots\\\
\cdots&\cdots&\cdots&\cdots\end{array}\right),$

acting on $l^{2}(\mathbb{Z})$ and $l^{2}(\mathbb{N})$ respectively. Note that
$L(a)$ is always normal whereas $T(a)$ need not be (see for example [18]). A
simple example already mentioned is $a(t)=t$ which gives rise to the bilateral
and unilateral shifts $L(a)=U_{1}$ and $T(a)=S$. In this case, both of these
operators are invariant under iterations of the IQR algorithm and hence their
finite sections $P_{m}Q_{n}^{*}TQ_{n}|_{P_{m}\mathcal{H}}$ always have
spectrum $\\{0\\}$. In the case of $L(a)$ this is an example of spectral
pollution, whereas in the case of $T(a)$ this does not capture the extremal
parts of the spectrum. Regarding pure finite section, the following beautiful
result is known:

######  Theorem 6.6 (Schmidt-Spitzer [66]).

If $a$ is a trigonometric polynomial then we have the following convergence in
the Hausdorff metric:

$\lim_{m\rightarrow\infty}\sigma(P_{m}L(a)|_{P_{m}\mathcal{H}})=\lim_{m\rightarrow\infty}\sigma(P_{m}T(a)|_{P_{m}\mathcal{H}})=\bigcap_{r\in(0,\infty)}\sigma(T(a_{r}))=:\Upsilon(a),$

where $a_{r}(t)=a(rt)$. Furthermore, this limit set is a connected finite
union of analytic arcs, each pair of which has at most endpoints in common.

It is straightforward to construct examples where it appears that both
$\lim_{n\rightarrow\infty}P_{m}Q_{n}^{*}T(a)Q_{n}|_{P_{m}\mathcal{H}}$ and
$\lim_{n\rightarrow\infty}P_{m}Q_{n}^{*}L(a)Q_{n}|_{P_{m}\mathcal{H}}$ exist
and are either the extreme parts of $\sigma(L(a))$ or of $\Upsilon(a)$. For
example consider the symbols

$a(t)=\frac{t^{3}+t^{-1}}{2},\quad\tilde{a}(t)=t+it^{-2}.$

Fig. 3 shows the outputs of the IQR algorithm and plain finite section for the
corresponding Laurent and Toeplitz operators for $m=50$ and $n=1$ and $n=300$.
In the case of $a$, it appears that both limit sets are the extremal parts of
$\sigma(L(a))$ (together with $0$ if $m$ is not a multiple of $4$). Whereas in
the case of $\tilde{a}$ it appears that
$\lim_{n\rightarrow\infty}P_{m}Q_{n}^{*}T(a)Q_{n}|_{P_{m}\mathcal{H}}$ is the
extremal parts of $\Upsilon(a)$ and
$\lim_{n\rightarrow\infty}P_{m}Q_{n}^{*}L(a)Q_{n}|_{P_{m}\mathcal{H}}$ is the
extremal parts of $\sigma(L(a))$ (again together with a finite collection of
points depending on the value of $m$ modulo $3$). Curiously, in both cases we
observed convergence in the strong operator topology to block diagonal
operators (up to unitary equivalence in each sublock), whose blocks have
spectra corresponding to the limiting sets (hence the dependence on remainder
of $m$ modulo $2$ or $3$). However, in contrast to convergence to points in
the discrete spectrum, convergence to these operators was only algebraic. This
is shown in Fig. 4 where we have plotted the Hausdorff distance between the
limiting set and the eigenvalues of the first diagonal block. We also shifted
the operators ($+1.1I$ for $a$ and $-1.5iI$ for $\tilde{a}$) so that the
extremal points correspond to exactly one point. In this scenario and for all
operators (Laurent or Toeplitz) the IQR algorithm converges strongly to a
diagonal operator whose diagonal entries are the corresponding extremal point
of $\sigma(L(a))$. This convergence is also shown in Fig. 4 and we observed a
slower rate of convergence than before. This is possibly due to points from
the other tips of the petals of $\sigma(L(a))$ converging as we increase $n$.
It would be interesting to see if some form of Theorem 6.6 holds for the IQR
algorithm (now taking $n\rightarrow\infty$). Given the examples presented
here, such a statement would likely be quite complicated. However, we
conjecture that if a normal operator has exactly one extreme point of its
essential spectrum (and finitely many eigenvalues of magnitude greater than
$r_{\mathrm{ess}}$) then this extreme point will be recovered in the limit
$n\rightarrow\infty$ for large enough $m$.

Figure 3: Top: Output of IQR and finite section on $T(a)$ and $L(a)$ for
$m=50$ and $n=1$ (left), $n=300$ (right). Bottom: Same but for the symbol
$\tilde{a}$. In both cases for a given symbol $b$, $\sigma(L(b))$ is given by
$\\{b(z):z\in\mathbb{T}\\}$ (shown) and $\sigma(T(b))$ is given by
$\sigma(L(b))\cup\\{z\in\mathbb{C}\backslash
b(\mathbb{T}):\mathrm{wind}(b,z)\neq 0\\}$.

Figure 4: Left: Algebraic convergence to block diagonal operators. Right:
Algebraic convergence to diagonal operators. In both cases we have plotted the
difference in eigenvalues of the first block as we increase $n$.

  * Example 6.7 (IQR and avoiding spectral pollution)

In this example we consider whether the IQR algorithm may be used as a tool to
avoid spectral pollution. Sometimes when considering
$\sigma(P_{m}T|_{P_{m}\mathcal{H}})$, spectral pollution can be detected by
changing $m$ (edge states which correspond to spectral pollution are often
unstable but this is not always the case). In general,
$\sigma(P_{m}Q_{n}^{*}TQ_{n}|_{P_{m}\mathcal{H}})$ can be considered as a
generalised version of finite section with a finite number ($n$) of IQR
iterates being performed on the infinite dimensional operator before
truncation. If $Q_{n}$ is unitary, then this simply changes the basis before
truncation and such a change may reduce (or change) spectral pollution
allowing it to be detected. Here we consider

$T_{3}=\left(\begin{matrix}0&3&&&\\\ 3&0&1&&\\\ &1&0&3&\\\ &&3&0&\ddots\\\
&&&\ddots&\ddots\\\ \end{matrix}\right).$

The spectrum of $T_{3}$ is $[-4,-2]\cup[2,4]$. However, if $m$ is odd then
$0\in\sigma(P_{m}T_{3}|_{P_{m}\mathcal{H}})$. We shifted the operator by
considering $T_{3}+0.2I$ (and then shifted back for the spectrum). Fig. 5
shows the Hausdorff distance between
$\sigma(P_{m}Q_{n}^{*}(T_{3}+0.2I)Q_{n}|_{P_{m}\mathcal{H}})-0.2I$ and
$\sigma(T_{3})$ as $n$ varies for different $m$. The spikes in the distance
correspond to eigenvalues leaving the interval $[-4,-2]$ and crossing to
$[2,4]$ (also shown in Fig. 5). The increase in distance as $m$ decreases (for
large $n$) is due less of the interval $[-4,-2]$ being approximated. It
appears that the IQR algorithm can be an effective tool at detecting spectral
pollution - certainly a mixture of varying $m$ and $n$ will be more effective
than just varying $m$.

Figure 5: Left:
$d_{H}(\sigma(P_{m}Q_{n}^{*}(T_{3}+0.2I)Q_{n}|_{P_{m}\mathcal{H}})-0.2I,\sigma(T_{3}))$
as a function of $n$ for different $m$. Right:
$\sigma(P_{m}Q_{n}^{*}(T_{3}+0.2I)Q_{n}|_{P_{m}\mathcal{H}})-0.2I$ as a
function of $n$ for $m=201$. Note the crossing of eigenvalues across the
spectral gap.

Another example of this is given by the operator $L(a)$ considered previously.
For fixed $m$ we found that

$\lim_{n\rightarrow\infty}\sup_{z\in\sigma(P_{m}Q_{n}^{*}L(a)Q_{n}|_{P_{m}\mathcal{H}})}\mathrm{dist}(z,\sigma(L(a)))=0.$

However, for finite section, spectral pollution occurs for all large $m$

$\lim_{m\rightarrow\infty}\sup_{z\in\sigma(P_{m}L(a)|_{P_{m}\mathcal{H}})}\mathrm{dist}(z,\sigma(L(a)))>0$

and the IQR algorithm can only recover the extreme parts of the spectrum

$\lim_{n\rightarrow\infty}d_{H}(\sigma(P_{m}Q_{n}^{*}L(a)Q_{n}|_{P_{m}\mathcal{H}}),\sigma(L(a)))>0.$

Despite this, we found that for small fixed $n>0$ it appears that

$\lim_{m\rightarrow\infty}d_{H}(\sigma(P_{m}Q_{n}^{*}L(a)Q_{n}|_{P_{m}\mathcal{H}}),\sigma(L(a)))=0.$

This is shown in Fig. 6 with similar results for $L(\tilde{a})$.

Figure 6: Left:
$d_{H}(\sigma(P_{m}Q_{n}^{*}L(a)Q_{n}|_{P_{m}\mathcal{H}}),\sigma(L(a)))$ as a
function of $n$ for different $m$.
$d_{H}(\sigma(P_{m}Q_{n}^{*}L(\tilde{a})Q_{n}|_{P_{m}\mathcal{H}}),\sigma(L(\tilde{a})))$
as a function of $n$ for different $m$.

###  6.3 Numerical examples II: non-normal operators

Although Theorem 3.9 considers normal operators, Theorems 3.13 and 3.15
suggest the IQR algorithm may also be useful for non-normal operators. Indeed,
the results presented here demonstrate that in practice the IQR algorithm can
work very well for non-normal problems. If an infinite matrix $T$ has $m$
isolated eigenvalues $\\{\lambda_{1},\ldots,\lambda_{m}\\}$ (repeated
according to multiplicity) outside $r_{\mathrm{ess}}(T)$ (the essential
spectral radius), then Theorems 3.13 and 3.15 suggest that the eigenvalues
will appear on the diagonal of $P_{m}Q_{n}^{*}TQ_{n}|_{P_{m}\mathcal{H}}$ as
$n\rightarrow\infty$, i.e.

$\sigma(P_{m}Q_{n}^{*}TQ_{n}|_{P_{m}\mathcal{H}})\longrightarrow\\{\lambda_{1},\ldots,\lambda_{m}\\},\qquad\text{as
}n\rightarrow\infty.$

We will verify this numerically in the next examples. However, we will see
that not only do we get convergence to the eigenvalues, but often we also pick
up parts of the boundary of the essential spectrum (this was the case when
considering $T(a)$ but appeared not to be the case for $T(\tilde{a})$). This
phenomenon is not accounted for in the previous exposition where normality was
crucial for proving Theorem 3.9.

  * Example 6.8 (Recovering the extremal part of the spectrum)

Let us return to the infinite matrices $A$ in (6.1) and $T$ in (6.3) from
Section 6.1. We have run the IQR algorithm with $n=1000$ and $n=300$ for $A$
and $T$ respectively, shown in Fig. 7. We see that if one takes a finite
section after running the IQR algorithm, then part of the boundary of the
essential spectrum also appears, along with the discrete spectrum
$\sigma_{d}(A)$. Note that the part of the boundary that is captured is the
extreme part (points with largest modulus). It seems that after running the
IQR algorithm, the spectral information from the largest isolated eigenvalues
and the largest approximate point spectrum is “squeezed up” to the upper and
leftmost portions of the matrix. This is not completely counter-intuitive
given (2.5) and is what normally happens in finite dimensions. For both
examples, we found that the IQR iterates converges to an upper triangular
matrix (analogous to the finite dimensional case) in agreement with Theorems
3.13 and 3.15. The convergence of the upper $1\times 1$ block for $A$
(corresponding to the dominant eigenvalue) and $4\times 4$ non-diagonal block
for $T$ are shown in Fig. 8 where we have plotted the difference in norm.

Figure 7: Left: Output of the IQR algorithm
$\sigma(P_{m}Q_{n}^{*}AQ_{n}|_{P_{m}})$ for $m=300$ and $n=1000$. Right:
Output of the IQR algorithm $\sigma(P_{m}Q_{n}^{*}TQ_{n}|_{P_{m}})$ for
$m=100$ and $n=300$.

We also discuss another drawback of the algorithm, $\Gamma_{m}$, to compute
the pseudospectrum by perturbing the operator $T$. Let $\widetilde{T}$ be the
operator obtained from $T$ if we set $\widetilde{T}_{5,4}=5\times 10^{7}$
(note that this gets rid of the block form). The computation of $\Gamma_{m}$
involves squaring the operator and hence leads to matrices of norm of order
$10^{15}$, making it impossible to compute $\sigma_{\epsilon}(\widetilde{T})$
using double precision for $\epsilon\lesssim 10^{-1}$. However, the IQR
algorithm shares the pleasant feature of finite section in allowing a wider
range of magnitudes of the matrix entries of the operator. The output for
$n=300$, $m=100$ is shown in Fig. 8 as well as convergence of the upper
$2\times 2$ block (corresponding to the dominant eigenvalues). Note in this
case we can only compute this upper block to an accuracy of about $10^{-8}$ in
double precision due to the large perturbed entry. However, this is still much
better than pseudospectral techniques. All the errors in Fig. 8 were obtained
via comparison with converged matrices computed using quadruple precision.

Figure 8: Left: Output of the IQR algorithm
$\sigma(P_{m}Q_{n}^{*}\widetilde{T}Q_{n}|_{P_{m}})$ for $m=100$ and $n=300$.
The reference circle is the boundary of the essential spectrum. Right:
Convergence of upper diagonal blocks for operators $A$, $T$ and
$\widetilde{T}$.

  * Example 6.9 ($PT$-symmetry in quantum mechanics)

Finally, we consider a so called $PT$ symmetric operator (non-normal),
demonstrating the same phenomena. A Hamiltonian $H=p^{2}/2+V(x)$ is said to be
$PT$ symmetric if it commutes with the action of the operator $PT$ where $P$
is the parity operator
$\hat{x}\rightarrow{}-\hat{x},\hat{p}\rightarrow{}-\hat{p}$ and $T$ the time
operator $\hat{p}\rightarrow{}-\hat{p},i\rightarrow{}-i$. Further distinction
can be made between exact (unbroken) $PT$ symmetry when $H$ shares common
“eigenfunctions” with $PT$ and broken $PT$ symmetry when they possess
different eigenfunctions. Many $PT$ Hamiltonians possess the remarkable
property that their spectra are real for small enough $\mathrm{Im}(V)$ but
that the spectrum becomes complex above a certain threshold [11]. This phase
transition from exact to broken $PT$ phase is known as symmetry breaking.
There has been a lot of interest in recent years, both theoretically and
experimentally, in non-Hermitian $PT$ symmetric Hamiltonians [64, 49].

We consider an operator on $l^{2}(\mathbb{Z})$ of the form

$(H_{1}x)_{n}=x_{n-1}+x_{n+1}+V_{n}x_{n}.$ (6.4)

This commutes with (the discrete version of) $PT$ precisely when the potential
has even real part and odd imaginary part. We tested the IQR algorithm on the
potential

$V_{n}=\begin{cases}\cos(n)+i\gamma\sin(n),&\mod(n,2)=0\\\ 0,&\mod(n,2)=1\\\
\end{cases},$ (6.5)

and found similar results for other potentials. Fig. 9 shows the same
qualitative behaviour as the last example for $\gamma=1,2$ at $m=500,n=3000$.
We shifted by $2.2$ and $2.15$ for $\gamma=1,2$ respectively. For comparison
we have shown converged resolvent norms. We found that spectral pollution with
no IQR iterates was consistent as we varied $m$. However, for a fixed $m$,
increasing the number of iterates ($n\rightarrow\infty$) caused
$\sigma(P_{m}Q_{n}^{*}H_{1}Q_{n}|_{P_{m}\mathcal{H}})$ to approach the
extremal part of the spectrum.

Figure 9: The figures show finite sections
$\sigma(P_{m}H_{1}|_{P_{m}\mathcal{H}})$ (magenta) and (shifted)
$\sigma(P_{m}Q_{n}^{*}H_{1}Q_{n}|_{P_{m}\mathcal{H}})$ IQR iterates (blue)
along with converged resolvent norm contours for $\gamma=1$ (left) and
$\gamma=2$ (right). Both figures are for $m=500,n=3000$ and show the
convergence to the extremal parts of the spectrum.

###  6.4 Numerical examples III: random non-Hermitian operators and boundary
conditions

In this final section, we explore examples where the
$P_{m}Q_{n}^{*}TQ_{n}P_{m}$ naturally give rise to periodic boundary
conditions (this was already seen for some examples of Laurent operators in
Section 6.2). Both examples discussed here are physically motivated random
tridiagonal operators on the lattice $\mathbb{Z}$. One of the key applications
of studying such random operators can be found in condensed matter physics.
The discrete models below have been used to study conductivity of disordered
media, flux lines in superconductors and asymmetric hopping particles. Many
such operators are also the discretisation of certain stochastic differential
equations. As we will demonstrate, the IQR method can be a powerful way of
avoiding spectral pollution caused by unnatural “open” boundary conditions in
forming the finite section $P_{m}TP_{m}$. In both of these examples, periodic
boundary conditions are natural and we find that taking finite sections after
iterating the IQR algorithm captures periodic boundary conditions.

  * Example 6.10 (Hopping sign model in sparse neural networks)

The first example is a non-normal operator with random sub and super-
diagonals, first studied by Feinberg and Zee [35, 44, 23]. The usual “Hopping
Sign Model” is defined via

$(H_{2}x)_{n}=x_{n-1}+b_{n}x_{n+1},$

with $b_{n}\in\\{\pm 1\\}$ (say independent Bernoulli with parameter $p=1/2$).
This describes a particle “‘hopping” on $\mathbb{Z}$ and can be mapped into a
(complex-valued) random walk. We will consider a slightly different operator
described by

$(H_{3}x)_{n}=s_{n-1}^{-}\exp({-g})x_{n-1}+s_{n}^{+}\exp({g})x_{n+1},$ (6.6)

and appearing in [1] in the context of sparse neural networks. We shall assume
that $g$ is real and non-negative and that $s_{j}^{\pm}$ are i.d.d. random
variables with Bernoulli distribution $p$. In other words

$\mathbb{P}(s_{j}^{\pm}=1)=1-\mathbb{P}(s_{j}^{\pm}=-1)=p.$

We will only consider $g=1/10$ and $p=1/2$, but will vary $p$ in an effort to
compute the spectrum of $H_{3}$ which only depends on the support of the
distribution of the $s^{\pm{}}_{j}$’s. It is easy to prove that the spectrum
(and pseudospectrum) of $H_{3}$ is almost surely constant and that there is no
inessential spectrum. Furthermore, one can show that $\sigma(H_{3})$ is
contained in the annulus $\\{z\in\mathbb{C}:2\sinh(g)\leq\left|z\right|\leq
2\cosh(g)\\}$.

Finite section calculations associated with this operator have some
interesting properties and are extensively studied in [1]. If one projects
using the standard basis of $l^{2}(\mathbb{Z})$ then one obtains matrices of
the form

$M_{n}^{1}=\left(\begin{matrix}0&s_{-n+1}^{-}\exp({-g})&&\\\
s_{-n+1}^{+}\exp({g})&0&\ddots&\\\ &\ddots&\ddots&s_{n-1}^{-}\exp({-g})\\\
&&s_{n-1}^{+}\exp({g})&0\\\ \end{matrix}\right).$

If we use open boundary conditions (i.e. we simply project onto the space
spanned by $\\{e_{-n},...,e_{n}\\}$) then one can “gauge” away $g$ by a
similarity transformation, leading to

$M^{\prime}_{n}=\left(\begin{matrix}0&s_{-n+1}^{-}&&\\\
s_{-n+1}^{+}&0&\ddots&\\\ &\ddots&\ddots&s_{n-1}^{-}\\\ &&s_{n-1}^{+}&0\\\
\end{matrix}\right).$

On the other hand, the use of periodic boundary conditions leads to the matrix

$M_{n}^{2}=\left(\begin{matrix}0&s_{-n+1}^{-}\exp({-g})&&s_{n}^{+}\exp({g})\\\
s_{-n+1}^{+}\exp({g})&0&\ddots&\\\ &\ddots&\ddots&s_{n-1}^{-}\exp({-g})\\\
s_{n}^{-}\exp({-g})&&s_{n-1}^{+}\exp({g})&0\\\ \end{matrix}\right),$

which does not suffer from this setback.

Figure 10: Top: Output of finite section over a random sample of $200$
matrices of size $200$ (left) and the estimates using pseudospectral
techniques (right). Bottom: The output of IQR over $200$ samples computing
$\sigma(P_{m}Q_{n}^{*}H_{3}Q_{n}|_{P_{m}\mathcal{H}})$ for $m=200$ and $n=50$
(left), $n=2000$ (right). Note that a few iterates seems to agree with
periodic boundary conditions and then increasing the number of iterates leads
to convergence to the extremal parts of the essential spectrum.

In [1] this phenomena was studied via localisation of the eigenvalues of
$M_{n}^{2}$, in particular using the Lyapunov exponent $\kappa(z)$ which is
equal to the inverse of the localisation length. An eigenfunction $\psi$ with
eigenvalue $z$ localised around $x_{0}$ behaves approximately as

$\left|\psi(x)\right|\sim\exp({-\kappa(z)\left|x-x_{0}\right|}).$

If one defines recursively

$y_{n+1}(z)=\exp(g)\frac{\psi_{n+2}}{\psi_{n+1}}=-(s^{-}_{n-1}/s_{n}^{+})/y_{n}(z)+z/s_{n}^{+}$

then (in the limit of large system sizes)

$\kappa(z;g)=\lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_{j=-N}^{N}\big{(}\log\left|y_{j}(z)\right|-g\big{)}.$

This is known as the transfer matrix approach. For fixed $z$, as we increase
$g$, $\kappa(z;g)$ becomes negative. The heuristic is that a hole opens up in
the spectrum corresponding to a negative Lyapunov exponent. Eigenvalues of
$M_{n}^{2}$ inside the hole are swept up and become delocalised moving to the
rim of the hole, whereas those outside remain largely undisturbed. Eigenvalues
of $M_{n}^{1}$ inside the negative $\kappa$ zone correspond to edge states due
to the finite system size approximation.

Fig. 10 shows the output of a sample of $200$ finite sections with open
boundary conditions and matrix size $200$. We have also shown the annular
region that bounds the spectrum, as well as the contour $\kappa=0$. In order
to calculate $\kappa$, we calculated the above sum on a grid with large $N$ to
ensure convergence. The colour bar corresponds to the inverse participation
ratio (log scale) of normalised eigenfunctions defined by

$1/P\equiv\frac{\sum_{j}\left|\psi_{i}\right|^{4}}{\sum_{j}\left|\psi_{i}\right|^{2}}.$

Note that this has a maximum value of $1$ (localised) and a minimum value of
$1/N$ (delocalised), $N$ being the size of the matrix. Open boundary
conditions produces spectral pollution in the hole with localised
eigenfunctions and the contour $\kappa=0$ corresponds to the delocalised
region. In order to compare to the spectrum of the infinite operator on
$l^{2}(\mathbb{Z})$ we have plotted $\sigma_{\epsilon}(H_{3})$, for
$\epsilon=10^{-2}$, calculated using matrix sizes of order $10^{5}$. We note
that the spectrum is independent of $p\in(0,1)$ so we have also shown the
union of these estimates over $p=\\{k/100\\}_{k=1}^{99}$. Although the
algorithm used to compute the pseudospectrum is guaranteed to converge to
$\sigma_{\epsilon}(H_{3})$, there are regions in the complex plane where this
convergence is very slow. Taking unions over $p$ is simply a way to speed up
this convergence. We found upon taking $\epsilon$ smaller that the spectrum
appeared to have a fractal like nature. It also appears that the hole in the
spectrum corresponds to the boundary of two ellipses. It is easy to prove that
the ellipse

$E_{1}=\\{\exp({g+i\theta})+\exp({-g-i\theta}):\theta\in[0,2\pi)\\}$

is contained in $\sigma(H_{3})$ and that the spectrum (and pseudospectrum) of
$H_{3}$ has fourfold rotational symmetry. Denoting the rotation of $E_{1}$ by
$\pi/4$ as $E_{2}$ we have shown $E_{1}\cup E_{2}$ in the figure.

Fig. 10 also shows the effect of IQR iterations over random samples of size
$200$ for $m=200$ and $n=50$ and $2000$. Remarkably, as we increase $n$, a few
iterations is enough to capture periodic boundary conditions and sweep away
the localised edge states. We have also shown the inverse participation ratio
which, although now is defined with respect to a new basis, still gives an
indication of how “diagonal” the matrix
$P_{m}Q_{n}^{*}H_{3}Q_{n}|_{P_{m}\mathcal{H}}$ is. If we increase $n$ further,
the output approaches the edge of the spectrum with eigenvectors becoming more
localised (in the new basis). We found exactly the same phenomena to occur if
we shifted the operator $H_{3}$ with convergence to the corresponding extremal
part of the essential spectrum.

  * Example 6.11 (NSA Anderson model in superconductors)

Finally, we consider a non-normal operator with no inessential spectrum where
the IQR algorithm does not seem to converge to the boundary of the essential
spectrum, but rather to a curve associated with periodic boundary conditions
in the large system size limit.

Over the past twenty years there has been considerable interest in non self-
adjoint random operators, sparked by Hatano and Nelson studying a non self-
adjoint Anderson model in the context of vortex pinning in type-II
superconductors [43]. Their model showed that an imaginary gauge field in a
disordered one-dimensional lattice can induce a delocalisation transition. The
operator in $\mathcal{B}(l^{2}(\mathbb{Z}))$ can be written as

$(H_{4}x)_{n}=\exp({-g})x_{n-1}+\exp({g})x_{n+1}+V_{n}x_{n}$ (6.7)

where $g>0$ and $V$ is a random potential. This operator also has applications
in population biology [56] and the self-adjoint version of this model is
widely studied for the phenomenon of Anderson localisation (absence of
diffusion of waves) [2, 12]. In the non self-adjoint case, complex values of
the spectrum indicate delocalisation. Note that we now have randomness on the
diagonal with fixed coupling coefficients $\exp({\pm g})$.

Standard finite section produces real eigenvalues since the matrix
$P_{m}H_{4}|_{P_{m}\mathcal{H}}$ is similar to a real symmetric matrix.
However, truncating the operator and adopting periodic boundary conditions
gives rise to the famous “bubble and wings”. If $V=0$ then the spectrum is an
ellipse $E=\\{\exp({g+i\theta})+\exp({-g-i\theta}):\theta\in[0,2\pi)\\}$, but
as we increase the randomness wings appear on the real axis. For a study of
this phenomenon and the described phase transition we refer the reader to
[35]. Goldsheid and Khoruzhenko have studied the convergence of the spectral
measure in the periodic case as $N\rightarrow\infty$ in [37], $N$ being the
number of sites. In general, the support of these measures as
$N\rightarrow\infty$ can be very different to the spectrum of the operator on
$l^{2}(\mathbb{Z})$ given by (6.7), highlighting the difficulty in computing
the spectrum.

We consider the case $g=1/2$ with $V_{n}$ i.i.d. Bernoulli random variables
taking values in $\\{\pm 1\\}$ with equal probability $p=1/2$. Again, there is
no inessential spectrum and the spectrum/pseudospectrum is constant almost
surely, depending only on the support of the distribution of the $V_{n}$. The
following inclusion is also known, which gives boundaries to the spectrum:

$\sigma(H_{4})\subset(\overline{\mathrm{conv}}(E)+[-1,1])\cap(E+B_{1}),$

where $\overline{\mathrm{conv}}(E)$ its closed convex hull of $E$ and $B_{1}$
denotes the closed unit disk. The choice of $g$ ensures the spectrum has a
hole in it. One may calculate the Lyapunov exponent, either by the transfer
matrix approach or by calculating a potential related to the density of
states. The limiting distribution of the eigenvalues of finite section with
periodic boundary conditions is given by the complex curve

$\\{z\in\mathbb{C}\backslash\mathbb{R}:\kappa(z)=0\\}\cup\\{x\in\mathrm{supp}(dN):\kappa(x+i0)>0\\}.$

The output of the IQR algorithm for $m=30$ and $n=15$ and $n=300$ over $200$
random samples are shown in Fig. 11. Note that if we took $n=0$, the spectrum
would be real in stark contrast to Fig. 11. Taking a small number of IQR
iterates approximates the bubble and wings with a few remaining real
eigenvalues. However, upon increasing $n$, the output does not seem to
converge to the extremal parts of the spectrum, but seems to remain stuck on
the limit curve with the operator
$P_{m}Q^{*}_{n}H_{4}Q_{n}|_{P_{m}\mathcal{H}}$. Shifting by $+4iI$ caused the
output to recover the top part of the limit curve.

Figure 11: The output of IQR over $200$ samples computing
$\sigma(P_{m}Q_{n}^{*}H_{4}Q_{n}|_{P_{m}\mathcal{H}})$ for $m=30$ and $n=15$
(left), $n=300$ (right). Note that we appear to recover the periodic limit
curve and increasing the number of iterates converges to the extremal parts.
Applying shifts allowed us to recover the extremal parts of the limit curves.

  * Remark 6.12

For any operator $T$ that has $Q_{n}$ unitary, the essential spectrum and
spectrum of $Q_{n}^{*}TQ_{n}$ is equal to that of $T$. As the above two
examples suggest, taking a small value of $n$ could be used as a method of
testing eigenvalues of finite section methods that correspond to finite system
size effects, such as open boundary conditions. This could be used in quasi-
periodic systems or systems with very few symmetries, where there is no
obvious choice of appropriate boundary conditions. However, detecting isolated
eigenvalues of finite multiplicity within the convex hull of the essential
spectrum still remains a challenge.

##  7 Concluding remarks and open problems

This paper discussed the generalisation of the famous QR algorithm to infinite
dimensions. It was shown that for a large class of operators, encompassing
many in scientific applications, the iterates of the IQR algorithm can be
computed efficiently on a computer. For matrices with finitely many entries in
each column, the computation collapses to a finite one. In general, for an
invertible operator we can compute the iterates to any given accuracy in
finite time. Furthermore, it was proven that for normal operators, the
algorithm converges to the discrete spectrum outside the convex hull of the
essential spectrum, with the rate of convergence generalising the well known
result in finite dimensions. These were extended to more general invariant
subspaces and non-normal operators in Theorems 3.13 and 3.15. Unfortunately
the IQR algorithm cannot in general be sped up with the use of shift
strategies, which considerably speed up the finite dimensional algorithm [61].
This is due to two reasons. The first is simply that there is no final column
of an infinite matrix, hence the usual link with inverse iteration cannot be
made. Second, it is also possible for part of the spectrum to be lost in the
limit (see Example 3.2) and below the essential spectra radius there is no
guarantee of convergence.

Despite these inherent drawbacks of the infinite dimensional setting, we
showed how the IQR algorithm can be used to gain new classification results
and convergent algorithms in the SCI hierarchy. In particular, we showed how
to compute eigenvalues and eigenspaces outside the essential spectrum with
error control for normal operators. This was extended to dominant invariant
subspaces for general (possibly non-normal) operators as well as the spectrum
of a large class of operators that includes compact normal operators with
eigenvalues of distinct magnitude. These results present the first such
algorithms that tackle these problems with error control.

Finally, we demonstrated that the IQR algorithm can be implemented both in
theory and in practice. We demonstrated the convergence theorems in Section 3
as well as some examples (normal and non-normal) where the extremal parts of
the essential spectrum also appear to be recovered. Based on this, we
conjecture that there may be a large class of operators for which the IQR
algorithm converges to the extreme parts of the essential spectrum.555This is
false in general as is easily seen by considering the shift operator. In
particular, we conjecture that this holds for normal operators if the set of
extremal points of the essential spectrum has size one. However, an example
was given where convergence to the essential spectrum was only algebraic
$O(n^{-\alpha})$ as $n\rightarrow\infty$ as opposed to the linear convergence
rate $O(r^{n})$ to the discrete spectrum/eigenvalues. It was also demonstrated
that the truncations of the IQR algorithm have spectra agreeing with periodic
boundary conditions for a range of operators in the class of “pseudoergodic”
NSA random operators. In some cases, the algorithm performed much better than
standard finite section methods. We should stress that, as in the case of the
finite section method, for fixed $n$ and $m\rightarrow\infty$, the output
$\sigma(P_{m}Q_{n}^{*}TQ_{n}|_{P_{m}})$ will in general still suffer from the
spectral pollution phenomenon and in some cases not recover the full spectrum.
This was apparent in the numerical examples and is likely to hold true for
many operators even when taking a mixture of double limits
$m,n\rightarrow\infty$. However, examples of Laurent operators were given
where it appears $\sigma(P_{m}Q_{n}^{*}TQ_{n}|_{P_{m}})$ converges to the
spectrum as $m\rightarrow\infty$ for fixed $n>0$ but not $n=0$. We hope that
the algorithm’s potential use in sifting out spectral pollution/complying with
appropriate boundary conditions via a canonical unitary transformation can
also be exploited.

Based on our findings, we end with a list of open problems for further study
on the theoretical properties of the IQR algorithm:

  * •

Which conditions are needed on a possibly non-normal operator in order for the
IQR algorithm to pick up the extreme points of the essential spectrum?

  * •

Is the convergence rate to non-isolated points of the spectrum algebraic?

  * •

For operators which do not have a trivial QR decomposition, is there a way of
choosing $n=n(m)$ such that $\sigma(P_{m}Q_{n(m)}^{*}TQ_{n(m)}|_{P_{{m}}})$
converges to the spectrum as $m\rightarrow\infty$? If not, then for which
classes of operators does such a choice exist?

  * •

Is there a link between the IQR algorithm and the finite section method with
periodic boundary conditions for the class of pseudoergodic operators?

  * •

Are there other cases where the IQR algorithm alleviates the need to provide
natural boundary conditions when applying the finite section method?

  * •

Extending the IQR algorithm to unbounded operators. Can the IQR algorithm also
be extended to a continuous version for differential operators?

## Acknowledgments

MJC acknowledges support from the UK Engineering and Physical Sciences
Research Council (EPSRC) grant EP/L016516/1. ACH acknowledges support from a
Royal Society University Research Fellowship as well as EPSRC grant
EP/L003457/1.We would also like to thank the referees whose comments and
suggestions led to the improvement of the manuscript.

##  Appendix A Appendix

###  A.1 Example codes

Here we show example code for the IQR algorithm in the case that the matrix
has $k$ subdiagonals. The code can easily be adapted for the more general case
considered in Section 4.1.

######  Algorithm A.1.

    % The Infinite_QR(A,n,k,m) takes a section P_{nk+m}AP_{nk+m}
    % of an infinite matrix A with k subdiagonals, performs n iterations
    % of the infinite dimensional QR algorithm and returns
    % J = P_mQ_nAQ*_nP_m.

    function J = Infinite_QR(A,n,k,m)
    d = size(A,2);
     for j=1:n
       A = Inf_QR(A,d-j*k,k);   % The output in each loop is actually
     end                        % U_(d-j*k)...U_1A_(j-1)U_1...U_(d-j*k)
    J = A(1:m,1:m);             % if A_j is the j-th term in the QR iteration.

######  Algorithm A.2.

    % Inf_QR(A,n,k) takes a matrix A with k subdiagonals and performs
    % multiplication by n Householder transformation from the left and
    % right, i.e. B = U_n...U_1AU_1...U_n.

    function B = Inf_QR(A,n,k)
    B = A; d = size(A,1);
     for j = 1:n
        u = House(A(j:j+k,j));
        A(j:j+k,j:d) = A(j:j+k,j:d) - 2*u*(u’*A(j:j+k,j:d));
        B(j:j+k,1:d) = B(j:j+k,1:d) - 2*u*(u’*B(j:j+k,1:d));
        B(1:d,j:j+k) = B(1:d,j:j+k) - 2*(B(1:d,j:j+k)*u)*u’;
     end

######  Algorithm A.3.

    % House(x) takes a vector x and creates a unit vector u
    % such that (I - 2u*u’)x = ce_1 where c is some complex
    % number (depending on x) and e_1 = [1,0...].

    function u = House(x)
    v = x;
    if v(1) == 0
       v(1) = v(1) + norm(v);             %This is the classical way
    else                                  %of creating Householder reflections
       v(1) = x(1) + sign(x(1))*norm(x);  %as in finite dimensions.
    end
    u = v/norm(v);

###  A.2 Recalling the basics of the SCI hierarchy

The corner stone in the SCI hierarchy is the definition of a computational
problem, a general algorithm and towers of algorithms. The basic objects in a
computational problem are as follows:

  * (i)

$\Omega$ is some set, called the _domain_.

  * (ii)

$\Lambda$ is a set of complex valued functions on $\Omega$ called the
_evaluation_ set.

  * (iii)

$\mathcal{M}$ is a metric space with metric $d_{\mathcal{M}}$.

  * (iv)

$\Xi:\Omega\to\mathcal{M}$ is called the _problem function_.

The set $\Omega$ is the set of objects that give rise to our computational
problems. The problem function $\Xi:\Omega\to\mathcal{M}$ is what we are
interested in computing. Moreover, the set $\Lambda$ is the collection of
functions that provide us with the information we are allowed to read. This
leads to the following definition.

######  Definition A.4 (Computational Problem).

Given a primary set $\Omega$, an evaluation set $\Lambda$, a metric space
$\mathcal{M}$ and a problem function $\Xi:\Omega\to\mathcal{M}$ we call the
collection $\\{\Xi,\Omega,\mathcal{M},\Lambda\\}$ a _computational problem_.

For instance, when computing the spectrum of bounded operators on
$l^{2}(\mathbb{N})$, we let $\Omega$ be a subset of
$\mathcal{B}(l^{2}(\mathbb{N}))$ (for example the set of self-adjoint
operators or compact operators), $(\mathcal{M},d)$ be the set of all non-empty
compact subsets of $\mathbb{C}$ provided with the Hausdorff metric $d=d_{H}$
in (1.2). The evaluation functions in $\Lambda$ consist of the family of all
functions $f_{i,j}:A\mapsto\langle Ae_{j},e_{i}\rangle$, $i,j\in\mathbb{N}$,
which provide the entries of the matrix representation of $A$ with respect to
the canonical basis $\\{e_{i}\\}_{i\in\mathbb{N}}$. Finally,
$\Xi:A\mapsto\sigma(A)$.

The goal is to find algorithms which approximate the function $\Xi$. More
generally, the main pillar of our framework is the concept of a tower of
algorithms, which is needed to describe problems that need several limits in
the computation. However, first one needs the definition of a general
algorithm.

######  Definition A.5 (General Algorithm).

Given a computational problem $\\{\Xi,\Omega,\mathcal{M},\Lambda\\}$, a
_general algorithm_ is a mapping $\Gamma:\Omega\to\mathcal{M}$ such that for
each $A\in\Omega$

  * (i)

there exists a finite subset of evaluations
$\Lambda_{\Gamma}(A)\subset\Lambda$,

  * (ii)

the action of $\,\Gamma$ on $A$ only depends on
$\\{A_{f}\\}_{f\in\Lambda_{\Gamma}(A)}$ where $A_{f}:=f(A),$

  * (iii)

for every $B\in\Omega$ such that $B_{f}=A_{f}$ for every
$f\in\Lambda_{\Gamma}(A)$, it holds that
$\Lambda_{\Gamma}(B)=\Lambda_{\Gamma}(A)$.

Note that the definition of a general algorithm is more general than the
definition of a Turing machine or a Blum-Shub-Smale (BSS) machine. A general
algorithm has no restrictions on the operations allowed. The only restriction
is that it can only take a finite amount of information, though it is allowed
to _adaptively_ choose the finite amount of information it reads depending on
the input. Condition (iii) assures that the algorithm reads the information in
a consistent way. Note that the purpose of such a general definition is to get
strong lower bounds. In particular, the more general the definition is, the
stronger a proven lower bound will be.

With a definition of a general algorithm we can define the concept of towers
of algorithms. However, before we define that, we will discuss the cases for
which we may have a set valued function.

  * Remark A.6 (Set valued functions)

Occasionally we will consider a function $\Xi$ such that for $T\in\Omega$ we
have that $\Xi(T)\subset\mathcal{M}$. In this case we will still require that
a general algorithm produces a single valued out put i.e
$\Gamma(T)\in\mathcal{M}$ for $T\in\Omega$. However, we replace the metric in
order to define convergence. In particular, $\Gamma_{n}(T)\rightarrow\Xi(T)$,
as $n\rightarrow\infty$ means

$\inf_{y\in\Xi(T)}d_{\mathcal{M}}(\Gamma_{n}(T),y)\rightarrow 0.$

######  Definition A.7 (Tower of Algorithms).

Given a computational problem $\\{\Xi,\Omega,\mathcal{M},\Lambda\\}$, a _tower
of algorithms of height $k$ for $\\{\Xi,\Omega,\mathcal{M},\Lambda\\}$_ is a
family of sequences of functions

$\Gamma_{n_{k}}:\Omega\rightarrow\mathcal{M},\
\Gamma_{n_{k},n_{k-1}}:\Omega\rightarrow\mathcal{M},\dots,\
\Gamma_{n_{k},\ldots,n_{1}}:\Omega\rightarrow\mathcal{M},$

where $n_{k},\ldots,n_{1}\in\mathbb{N}$ and the functions
$\Gamma_{n_{k},\ldots,n_{1}}$ at the “lowest level” of the tower are general
algorithms in the sense of Definition A.5. Moreover, for every $A\in\Omega$,

$\Xi(A)=\lim_{n_{k}\rightarrow\infty}\Gamma_{n_{k}}(A),\quad\Gamma_{n_{k},\ldots,n_{j+1}}(A)=\lim_{n_{j}\rightarrow\infty}\Gamma_{n_{k},\ldots,n_{j}}(A)\quad
j=k-1,\dots,1.$

In addition to a general tower of algorithms (defined above), we will focus on
radical towers. The definition of a general algorithm allows for strong lower
bounds, however, to produce upper bounds we must add structure to the
algorithm and towers of algorithms. A radical tower allows for arithmetic
operations, comparisons and radicals.

######  Definition A.8 (Radical Towers).

Given a computational problem $\\{\Xi,\Omega,\mathcal{M},\Lambda\\}$, a
_Radical Tower of Algorithms_ of height $k$ for
$\\{\Xi,\Omega,\mathcal{M},\Lambda\\}$ is a tower of algorithms where the
lowest level functions

$\Gamma=\Gamma_{n_{k},\ldots,n_{1}}:\Omega\rightarrow\mathcal{M}$

satisfy the following: For each $A\in\Omega$ the action of $\,\Gamma$ on $A$
consists of only finitely many arithmetic operations, comparisons and radicals
($\sqrt{\cdot}$) of positive numbers on
$\\{A_{f}\\}_{f\in\Lambda_{\Gamma}(A)}$, where $A_{f}=f(A)$.

In other words one may say that for the finitely many steps of the computation
of the lowest functions
$\Gamma=\Gamma_{n_{k},\ldots,n_{1}}:\Omega\rightarrow\mathcal{M}$ only the
four arithmetic operations $+,-,\cdot,/$ within the smallest (algebraic) field
which is generated by the input $\\{A_{f}\\}_{f\in\Lambda_{\Gamma}(A)}$ are
allowed. In addition we allow to extract radicals of positive real numbers. We
implicitly assume that any complex number can be decomposed into a real and an
imaginary part, and moreover we can determine whether $a=b$ or $a>b$ for all
real numbers $a,b$ which can occur during the computations. Given the
definitions above we can now define the key concept, namely, the Solvability
Complexity Index:

######  Definition A.9 (Solvability Complexity Index).

A computational problem $\\{\Xi,\Omega,\mathcal{M},\Lambda\\}$ is said to have
_Solvability Complexity Index
$\mathrm{SCI}(\Xi,\Omega,\mathcal{M},\Lambda)_{\alpha}=k$_, with respect to a
tower of algorithms of type $\alpha$, if $k$ is the smallest integer for which
there exists a tower of algorithms of type $\alpha$ of height $k$. If no such
tower exists then
$\mathrm{SCI}(\Xi,\Omega,\mathcal{M},\Lambda)_{\alpha}=\infty.$ If there
exists a tower $\\{\Gamma_{n}\\}_{n\in\mathbb{N}}$ of type $\alpha$ and height
one such that $\Xi=\Gamma_{n_{1}}$ for some $n_{1}<\infty$, then we define
$\mathrm{SCI}(\Xi,\Omega,\mathcal{M},\Lambda)_{\alpha}=0$. We may sometimes
write $\mathrm{SCI}(\Xi,\Omega)_{\alpha}$ to simplify notation when
$\mathcal{M}$ and $\Lambda$ are obvious.

The definition of the SCI immediately induces the SCI hierarchy:

######  Definition A.10 (The Solvability Complexity Index Hierarchy).

Consider a collection $\mathcal{C}$ of computational problems and let
$\mathcal{T}$ be the collection of all towers of algorithms of type $\alpha$
for the computational problems in $\mathcal{C}$. Define

$\begin{split}\Delta^{\alpha}_{0}&:=\\{\\{\Xi,\Omega\\}\in\mathcal{C}\ |\
\mathrm{SCI}(\Xi,\Omega)_{\alpha}=0\\}\\\
\Delta^{\alpha}_{m+1}&:=\\{\\{\Xi,\Omega\\}\in\mathcal{C}\ |\
\mathrm{SCI}(\Xi,\Omega)_{\alpha}\leq m\\},\qquad\quad
m\in\mathbb{N},\end{split}$

as well as

$\Delta^{\alpha}_{1}:=\\{\\{\Xi,\Omega\\}\in\mathcal{C}\ |\ \exists\
\\{\Gamma_{n}\\}_{n\in\mathbb{N}}\in\mathcal{T}\text{ s.t. }\forall A\
d(\Gamma_{n}(A),\Xi(A))\leq 2^{-n}\\}.$

  * Remark A.11 (The $\Delta_{k}$ notation)

Note that in this paper we only consider radical towers and hence the
superscript $\alpha$ will be omitted throughout. Thus we will always write
$\Delta_{k}$.

Finally, we recall the definition of $\Sigma_{1}^{\alpha}$.

$\begin{split}\Sigma^{\alpha}_{1}&=\\{\\{\Xi,\Omega\\}\in\Delta^{\alpha}_{2}\
|\ \exists\ \\{\Gamma_{n}\\}_{n\in\mathbb{N}}\in\mathcal{T}\text{ s.t. }\
\Gamma_{n}(A)\subset\mathcal{N}_{2^{-n}}(\Xi(A))\ \text{and}\
\Gamma_{n}(A)\rightarrow\Xi(A)\ \forall A\in\Omega\\}\\\ \end{split}$

where $\mathcal{N}_{\delta}(\omega)$ denotes the $\delta$-neighbourhood of
$\omega\subset\mathcal{M}.$

## References

  * [1] A. Amir, N. Hatano, and D. R. Nelson.  Non-Hermitian localization in biological networks.  Physical Review E, 93(4):042310, 2016.
  * [2] P. W. Anderson.  Absence of diffusion in certain random lattices.  Physical review, 109(5):1492, 1958.
  * [3] N. Aronszajn.  Approximation methods for eigenvalues of completely continuous symmetric operators.  In Proceedings of the Symposium on Spectral Theory and Differential Problems, pages 179–202, 1951.
  * [4] W. Arveson.  Improper filtrations for $C^{*}$-algebras: spectra of unilateral tridiagonal operators.  Acta Sci. Math. (Szeged), 57(1-4):11–24, 1993.
  * [5] W. Arveson.  Noncommutative spheres and numerical quantum mechanics.  In Operator algebras, mathematical physics, and low-dimensional topology (Istanbul, 1991), volume 5 of Res. Notes Math., pages 1–10. A K Peters, Wellesley, MA, 1993.
  * [6] W. Arveson.  $C^{*}$-algebras and numerical linear algebra.  J. Funct. Anal., 122(2):333–360, 1994.
  * [7] W. Arveson.  The role of $C^{\ast}$-algebras in infinite-dimensional numerical linear algebra.  In $C^{\ast}$-algebras: 1943–1993 (San Antonio, TX, 1993), volume 167 of Contemp. Math., pages 114–129. Amer. Math. Soc., Providence, RI, 1994.
  * [8] J. Ben-Artzi, M. J. Colbrook, A. C. Hansen, O. Nevanlinna, and M. Seidel.  Computing Spectra–On the Solvability Complexity Index Hierarchy and Towers of Algorithms.  arXiv preprint arXiv:1508.03280, 2015.
  * [9] J. Ben-Artzi, A. C. Hansen, O. Nevanlinna, and M. Seidel.  New barriers in complexity theory: On the solvability complexity index and the towers of algorithms.  Compt. Rend. Math., 353(10):931 – 936, 2015.
  * [10] C. M. Bender.  Making sense of non-Hermitian Hamiltonians.  Reports on Progress in Physics, 70(6):947, 2007.
  * [11] C. M. Bender and S. Boettcher.  Real spectra in non-Hermitian Hamiltonians having PT symmetry.  Physical Review Letters, 80(24):5243, 1998.
  * [12] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, and A. Aspect.  Direct observation of Anderson localization of matter waves in a controlled disorder.  Nature, 453(7197):891–894, 2008.
  * [13] S. Bögli, B. M. Brown, M. Marletta, C. Tretter, and M. Wagenhofer.  Guaranteed resonance enclosures and exclosures for atoms and molecules.  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 470(2171), 2014.
  * [14] A. Böttcher.  Pseudospectra and singular values of large convolution operators.  J. Integral Equations Appl., 6(3):267–301, 1994.
  * [15] A. Böttcher.  Infinite matrices and projection methods.  In Lectures on operator theory and its applications (Waterloo, ON, 1994), volume 3 of Fields Inst. Monogr., pages 1–72. Amer. Math. Soc., Providence, RI, 1996.
  * [16] A. Böttcher, H. Brunner, A. Iserles, and S. P. Nørsett.  On the singular values and eigenvalues of the Fox-Li and related operators.  New York J. Math., 16:539–561, 2010.
  * [17] A. Böttcher, A. V. Chithra, and M. N. N. Namboodiri.  Approximation of approximation numbers by truncation.  Integral Equations Operator Theory, 39(4):387–395, 2001.
  * [18] A. Böttcher and B. Silbermann.  Introduction to large truncated Toeplitz matrices.  Universitext. Springer-Verlag, New York, 1999.
  * [19] A. Böttcher and I. M. Spitkovsky.  A gentle guide to the basics of two projections theory.  Linear Algebra and its Applications, 432(6):1412–1459, 2010.
  * [20] L. Boulton.  Projection methods for discrete Schrödinger operators.  Proc. London Math. Soc. (3), 88(2):526–544, 2004.
  * [21] N. Brown.  Quasi-diagonality and the finite section method.  Math. Comp., 76(257):339–360, 2007.
  * [22] H. Brunner, A. Iserles, and S. P. Nørsett.  The computation of the spectra of highly oscillatory Fredholm integral operators.  J. Integral Equations Appl., 23(4):467–519, 2011.
  * [23] S. Chandler-Wilde, R. Chonchaiya, and M. Lindner.  Eigenvalue problem meets Sierpinski triangle: computing the spectrum of a non–self–adjoint random operator.  Operators and Matrices, 5(4):633–648, 2011.
  * [24] M. J. Colbrook.  Computing spectral measures and spectral types.  arXiv:1908.06721v2, 2019.
  * [25] M. J. Colbrook.  The foundations of spectral computations via the solvability complexity index hierarchy: Part II.  arXiv:1908.09598, 2019.
  * [26] M. J. Colbrook.  Pseudoergodic operators and periodic boundary conditions.  Mathematics of Computation, 89(322):737–766, 2020.
  * [27] M. J. Colbrook and A. C. Hansen.  The foundations of spectral computations via the solvability complexity index hierarchy: Part I.  arXiv:1908.09592, 2019.
  * [28] M. J. Colbrook, B. Roman, and A. C. Hansen.  How to compute spectra with error control.  Physical Review Letters, 122(25):250201, 2019.
  * [29] E. B. Davies.  Spectral enclosures and complex resonances for general self-adjoint operators.  LMS J. Comput. Math., 1:42–74, 1998.
  * [30] E. B. Davies.  Linear operators and their spectra, volume 106.  Cambridge University Press, 2007.
  * [31] C. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, et al.  Hofstadter’s butterfly and the fractal quantum Hall effect in moire superlattices.  Nature, 497(7451):598–602, 2013.
  * [32] P. Deift, L. Li, and C. Tomei.  Toda flows with infinitely many variables.  Journal of functional analysis, 64(3):358–402, 1985.
  * [33] T. Digernes, V. S. Varadarajan, and S. Varadhan.  Finite approximations to quantum systems.  Reviews in Mathematical Physics, 6(04):621–648, 1994.
  * [34] P. Doyle and C. McMullen.  Solving the quintic by iteration.  Acta Math., 163(3-4):151–180, 1989.
  * [35] J. Feinberg and A. Zee.  Non-Hermitian localization and delocalization.  Physical Review E, 59(6):6433, 1999.
  * [36] M. C. T. for MATLAB 4.5.3.12856.  Advanpix LLC., Yokohama, Japan.
  * [37] I. Y. Goldsheid and B. A. Khoruzhenko.  Distribution of eigenvalues in non-Hermitian Anderson models.  Physical review letters, 80(13):2897, 1998.
  * [38] R. M. Gray et al.  Toeplitz and circulant matrices: A review.  Foundations and Trends in Communications and Information Theory, 2(3):155–239, 2006.
  * [39] R. Hagen, S. Roch, and B. Silbermann.  $C^{*}$-algebras and numerical analysis, volume 236 of Monographs and Textbooks in Pure and Applied Mathematics.  Marcel Dekker Inc., New York, 2001.
  * [40] A. C. Hansen.  On the approximation of spectra of linear operators on Hilbert spaces.  J. Funct. Anal., 254(8):2092–2126, 2008.
  * [41] A. C. Hansen.  Infinite-dimensional numerical linear algebra: theory and applications.  Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 466(2124):3539–3559, 2010.
  * [42] A. C. Hansen.  On the solvability complexity index, the n-pseudospectrum and approximations of spectra of operators.  J. Amer. Math. Soc., 24(1):81–124, 2011.
  * [43] N. Hatano and D. R. Nelson.  Localization transitions in non-Hermitian quantum mechanics.  Physical review letters, 77(3):570, 1996.
  * [44] D. E. Holz, H. Orland, and A. Zee.  On the remarkable spectrum of a non-Hermitian random matrix model.  Journal of Physics A: Mathematical and General, 36(12):3385, 2003\.
  * [45] T. Kato.  Perturbation theory for linear operators, volume 132.  Springer Science & Business Media, 2013.
  * [46] M. Krein and M. Krasnoselski.  Fundamental theorems concerning the extension of Hermitian operators and some of their applications to the theory of orthogonal polynomials and the moment problem.  Uspekhi Mat. Nauk., 2:60–106, 1947.
  * [47] M. Levitin and E. Shargorodsky.  Spectral pollution and second-order relative spectra for self-adjoint operators.  IMA Journal of Numerical Analysis, 24(3):393–416, 2004.
  * [48] M. Lindner.  Infinite matrices and their finite sections.  Frontiers in Mathematics. Birkhäuser Verlag, Basel, 2006.  An introduction to the limit operator method.
  * [49] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani.  Beam dynamics in PT symmetric optical lattices.  Physical Review Letters, 100(10):103904, 2008.
  * [50] M. Marletta.  Neumann-Dirichlet maps and analysis of spectral pollution for non-self-adjoint elliptic PDEs with real essential spectrum.  IMA J. Numer. Anal., 30(4):917–939, 2010.
  * [51] M. Marletta and R. Scheichl.  Eigenvalues in spectral gaps of differential operators.  J. Spectr. Theory, 2(3):293–320, 2012.
  * [52] D. C. Mattis.  The few-body problem on a lattice.  Reviews of Modern Physics, 58(2):361, 1986.
  * [53] C. McMullen.  Families of rational maps and iterative root-finding algorithms.  Ann. of Math. (2), 125(3):467–493, 1987.
  * [54] C. McMullen.  Braiding of the attractor and the failure of iterative algorithms.  Invent. Math., 91(2):259–272, 1988.
  * [55] A. Mogilner.  Hamiltonians in solid state physics as multiparticle discrete Schrödinger operators.  Advances in Societ Mathematics, 5:139–194, 1991.
  * [56] D. R. Nelson and N. M. Shnerb.  Non-Hermitian localization and population biology.  Physical Review E, 58(2):1383, 1998.
  * [57] S. Olver.  ApproxFun.jl v0.8.  github (online) https://github.com/JuliaApproximation/ApproxFun.jl, 2018.
  * [58] S. Olver and A. Townsend.  A Fast and Well-Conditioned Spectral Method.  SIAM Rev., 55(3):462–489, 2013.
  * [59] S. Olver and A. Townsend.  A Practical Framework for Infinite-dimensional Linear Algebra.  In Proceedings of the 1st First Workshop for High Performance Technical Computing in Dynamic Languages, HPTCDL ’14, pages 57–62, Piscataway, NJ, USA, 2014. IEEE Press.
  * [60] S. Olver and M. Webb.  SpectralMeasures.jl.  github (online) https://github.com/JuliaApproximation/SpectralMeasures.jl, 2018.
  * [61] B. N. Parlett.  The symmetric eigenvalue problem, volume 20.  siam, 1998.
  * [62] A. Pokrzywa.  Method of orthogonal projections and approximation of the spectrum of a bounded operator.  Studia Math., 65(1):21–29, 1979.
  * [63] L. Ponomarenko, R. Gorbachev, G. Yu, D. Elias, R. Jalil, A. Patel, A. Mishchenko, A. Mayorov, C. Woods, J. Wallbank, et al.  Cloning of Dirac fermions in graphene superlattices.  Nature, 497(7451):594–597, 2013.
  * [64] A. Regensburger, C. Bersch, M.-A. Miri, G. Onishchukov, D. N. Christodoulides, and U. Peschel.  Parity-time synthetic photonic lattices.  Nature, 488(7410):167–171, 2012.
  * [65] R. Riddell.  Spectral concentration for self-adjoint operators.  Pacific Journal of Mathematics, 23(2):377–401, 1967.
  * [66] P. Schmidt and F. Spitzer.  The Toeplitz matrices of an arbitrary Laurent polynomial.  Mathematica Scandinavica, 8(1):15–38, 1960.
  * [67] M. Seidel.  On $(N,\epsilon)$-pseudospectra of operators on Banach spaces.  J. Funct. Anal., 262(11):4916–4927, 2012.
  * [68] M. Seidel and B. Silbermann.  Finite sections of band-dominated operators – norms, condition numbers and pseudospectra.  In Operator theory, pseudo-differential equations, and mathematical physics, volume 228 of Oper. Theory Adv. Appl., pages 375–390. Birkhauser/Springer Basel AG, Basel, 2013.
  * [69] E. Shargorodsky.  Geometry of higher order relative spectra and projection methods.  J. Operator Theory, 44(1):43–62, 2000.
  * [70] E. Shargorodsky.  On the limit behaviour of second order relative spectra of self-adjoint operators.  J. Spectr. Theory, 3:535–552, 2013.
  * [71] P. Shivakumar, K. Sivakumar, and Y. Zhang.  Infinite matrices and their recent applications, 2016.
  * [72] B. Simon.  The classical moment problem as a self-adjoint finite difference operator.  Advances in Mathematics, 137(1):82–203, 1998.
  * [73] S. Smale.  The fundamental theorem of algebra and complexity theory.  Bull. Amer. Math. Soc. (N.S.), 4(1):1–36, 1981.
  * [74] A. Szabo and N. S. Ostlund.  Modern quantum chemistry: introduction to advanced electronic structure theory.  Courier Corporation, 2012.
  * [75] G. Teschl.  Jacobi operators and completely integrable nonlinear lattices.  American Mathematical Soc., 2000.
  * [76] L. N. Trefethen and M. Embree.  Spectra and pseudospectra: the behavior of nonnormal matrices and operators.  Princeton University Press, 2005.
  * [77] M. Webb and S. Olver.  Spectra of Jacobi operators via connection coefficient matrices.  arXiv preprint arXiv:1702.03095, 2017.