On the infinite-dimensional QR algorithm

Spectral computations of infinite-dimensional operators are notoriously difficult, yet ubiquitous in the sciences. Indeed, despite more than half a century of research, it is still unknown which classes of operators allow for the computation of spectra and eigenvectors with convergence rates and error control. Recent progress in classifying the difficulty of spectral problems into complexity hierarchies has revealed that the most difficult spectral problems are so hard that one needs three limits in the computation, and no convergence rates nor error control is possible. This begs the question: which classes of operators allow for computations with convergence rates and error control? In this paper, we address this basic question, and the algorithm used is an infinite-dimensional version of the QR algorithm. Indeed, we generalise the QR algorithm to infinite-dimensional operators. We prove that not only is the algorithm executable on a finite machine, but one can also recover the extremal parts of the spectrum and corresponding eigenvectors, with convergence rates and error control. This allows for new classification results in the hierarchy of computational problems that existing algorithms have not been able to capture. The algorithm and convergence theorems are demonstrated on a wealth of examples with comparisons to standard approaches (that are notorious for providing false solutions). We also find that in some cases the IQR algorithm performs better than predicted by theory and make conjectures for future study.

In the paper "On the Solvability Complexity Index, the n-pseudospectrum and approximations of spectra of operators" [38] the Solvability Complexity Index (SCI) was introduced. The SCI provides a classification hierarchy [8,9,38] of spectral problems according to their computational difficulty. The SCI of a class of spectral problems is the least number of limits needed in order to compute the spectrum of operators in this class. From a classical numerical analysis point of view, such a concept may seem foreign. Indeed, the traditional sentiment is that one should have an algorithm, n , such that for an operator T ∈ B(l 2 (N)), preferably with some form of error control of the convergence. As this philosophy forms the basics of numerical analysis, it naturally permeates the classical literature on the computational spectral problem. However, as is shown in [8,9,38], an algorithm satisfying (1.1) is impossible even for the class of self-adjoint operators. Indeed, in the general case, the best possible alternative is an algorithm depending on three indices n 1 , n 2 , n 3 such that In fact, any algorithm with fewer than three limits will fail on the general class of operators. Moreover, no error control nor convergence rate on any of the limits are possible, since any such error control would reduce the number of limits needed. However, for the self-adjoint and normal cases, two limits suffice in order to recover the spectrum. This phenomenon implies that the only way to characterise the computational spectral problem is through a hierarchy classifying the difficulty of computing spectra of different subclasses of operators. This is the motivation behind the SCI hierarchy, which also covers general numerical analysis problems. Indeed, the SCI hierarchy is closely related to Smale's question on the existence of purely iterative generally convergent algorithm for polynomial zero finding [69]. As demonstrated by McMullen [49,50] and Doyle and McMullen [30], this is a case where several limits are needed in the computation, and their results become special cases of classification in the SCI hierarchy [8,9]. Informally, the SCI hierarchy is characterised as follows (see the "Appendix 1" for a more detailed summary describing the SCI hierarchy). 0 : The set of problems that can be computed in finite time, the SCI = 0. 1 : The set of problems that can be computed using one limit, the SCI = 1, however, one has error control and one knows an error bound that tends to zero as the algorithm progresses. 2 : The set of problems that can be computed using one limit, the SCI = 1, but error control may not be possible. m+1 : For m ∈ N, the set of problems that can be computed by using m limits, the SCI ≤ m.
The class 1 is of course a highly desired class, however, most spectral problems are much higher in the hierarchy. For example, we have the following known classifications [8,9,38].
(i) The general spectral problem is in 4 \ 3 .
(iii) The compact spectral problem is in 2 \ 1 .
Here, the notation\indicates the standard "setminus". Note that the SCI hierarchy can be refined. We will not consider the full generalisation in the higher part of the hierarchy in this paper, but recall the class 1 [24]. This class is defined as follows. 1 : We have 1 ⊂ 1 ⊂ 2 and 1 is the set of problems that can be computed by passing to one limit. Error control may not be possible, however, there exists an algorithm for these problems that converges and for which its output is included in the spectrum (up to an arbitrarily small accuracy parameter ).
In the context of computing σ (T ), a 1 classification means the existence of an algorithm n such that n (T ) ⊂ σ (T ) + B 2 −n (0) and n converges to σ (T ) in the Hausdorff metric. The 1 class is very important as it allows for algorithms that never make a mistake. In particular, one is always sure that the output is sound but we do not know if we have everything yet. The simplest infinite-dimensional spectral problem is that of computing the spectrum of an infinite diagonal matrix and, as is easy to see, we have the following.
(iv) The problem of computing spectra of infinite diagonal matrices is in 1 \ 1 .
Hence, the computational spectral problem becomes an infinite classification theory in order to characterise the above hierarchy. In order to do so, there will, necessarily, have to be many different types of algorithms. Indeed, characterising the hierarchy will yield a myriad of different approaches, as different structures on the various classes of operators will require specific algorithms. The key contribution of this paper is to investigate the convergence properties of the infinite-dimensional QR (IQR) algorithm, its implementation properties, and how this algorithm provides classification results in the SCI hierarchy.

Main contribution and novelty of the paper
The main contributions of the paper can be summarised as follows: New convergence results, algorithmic results (the IQR algorithm can be implemented), classification results in the SCI hierarchy and numerical examples.
(1) Convergence results We provide new convergence theorems for the IQR algorithm with convergence rates and error control. The results include eigenvalues, eigenvectors and invariant subspaces. (2) Algorithmic implementation We prove that for infinite matrices with finitely many non-zero entries in each column, it is possible to implement the IQR algorithm exactly (on a finite machine) as if one had an infinite computer at one's disposal. This can be extended to implementing the IQR algorithm with error control for general invertible operators. (3) SCI hierarchy classifications As a result of (1) and (2), we provide new classification results for the SCI hierarchy. In particular, the convergence properties of the IQR algorithm capture key structures that allow for sharp 1 classification of the problem of computing extremal points in the spectrum. Moreover, we establish sharp 1 classification of the problem of computing spectra of subclasses of compact operators. (4) Numerical examples Finally, we demonstrate the IQR algorithm and the proven convergence results on a variety of difficult problems in practical computation, illustrating how the IQR algorithm is much more than a theoretical concept. Moreover, the examples demonstrate that the IQR algorithm performs much better than predicted by our theory, working on much larger classes of operators.

Connection to previous work
Our results connect to many different approaches in the vast literature on spectral computation in infinite dimensions. The infinite-dimensional computational spectral problem is very different from the finite-dimensional computational eigenvalue problem, and even though the IQR algorithm is inspired by the finite-dimensional version, this paper solely focuses on the infinite-dimensional problem. Thus, the paper is aimed at the analysis and numerical analysis audience focusing on infinite-dimensional problems rather than the finite-dimensional numerical linear algebra discipline.

Finite sections
The IQR algorithm provides an alternative to the standard finite section method in several cases where it fails. Whereas the finite section method would extract a finite section from the infinite matrix and then apply, for example, the finite-dimensional QR algorithm, the IQR algorithm first performs the infinite QR iterations and then extracts a finite section. In general, these two processes do not commute. The finite section method (or any derivative of it) cannot work in general because of the general classification results in the SCI hierarchy mentioned in Sect. 1. Typically, it may provide false solutions. However, in the cases where it converges, it provides invaluable 2 classifications in the SCI hierarchy. The finite section method has often been viewed in connection with Toeplitz theory and the reader may want to consult the work by Böttcher [14,15], Böttcher and Silberman [18], Böttcher et al. [16], Brunner et al. [22], Hagen et al. [35], Lindner [44], Marletta [46] and Marletta and Scheichl [47]. From the operator algebra point of view, the work of Arveson [5][6][7] has been influential as well as the work of Brown [21].
Infinite-dimensional Toda flow Deift et al. [28] provided the first results on the IQR algorithm in connection with Toda flows with infinitely many variables. Their results are purely functional analytic and do not take implementation and computability issues into account. However, these results provide the fundamentals of the IQR algorithm. In [36] these results were expanded with a convergence result for eigenvectors corresponding to eigenvalues outside the essential numerical range for normal operators. Yet, this paper did not consider convergence rates, actual numerical calculation nor any classification results.
Infinite-dimensional QL algorithm Olver, Townsend and Webb have provided a practical framework for infinite-dimensional linear algebra and foundational results on computations with infinite data structures [53][54][55][56]73]. This includes efficient codes as well as theoretical results. The infinite-dimensional QL (IQL) algorithm is an important part of this program. The IQL algorithm is rather different from the IQR algorithm, although they are similar in spirit. In particular, both the implementation and the convergence results are somewhat contrasting.
Infinite-dimensional spectral computation: The results in this paper follow in the long tradition of infinite-dimensional spectral computations. This field contains a vast literature that spans more than half a century, and the references that we have cited in the first paragraph of Sect. 1 represent a small sample. However, we would like to highlight the recent work by Bögli et al. [13] who were able to computationally confirm, with absolute certainty, a conjecture on a certain oscillatory behaviour of higher auto-ionizing resonances of atoms. Note that problems that are classified as 1 and 1 in the SCI hierarchy may allow for computer assisted proofs.

Background and notation
Here we briefly recall some definitions used in the paper. We will consider the canonical separable Hilbert space H = l 2 (N) (the set of square summable sequences). Moreover, we write B(H) for the set of bounded operators on H. For orthogonal projections E, F, we will write E ≤ F if the range of E is a subspace of the range of F. We denote the canonical orthonormal basis of H by {e j } j∈N , and if ξ ∈ H we write ξ( j) = ξ, e j .
Note that T ∈ B(H) is uniquely determined by its matrix elements t i j = T e j , e i . Hence we will use the words bounded operator and infinite matrix interchangeably. Given a sequence of operators {T n }, we will use the notation to mean convergence in the strong and weak operator topology respectively. The spectrum of T ∈ B(H) will be denoted by σ (T ), and σ d (T ) denotes the set of isolated eigenvalues with finite multiplicity (the discrete spectrum).
In connection with the spectrum, we need to recall some definitions which will appear in the statement of our theorems. We recall that, for T ∈ B(H), the essential spectrum 1 and the essential spectral radius are given by Moreover, the numerical range and the essential numerical range of T are defined by In addition, we need the Hausdorff metric as defined by the following. Let S, T ⊂ C, be compact. Then their Hausdorff distance is where d(λ, T ) = inf ρ∈T |ρ − λ|. We also recall a generalisation of the spectrum, known as the pseudospectrum. Indeed, for > 0 define the -pseudospectrum as where we interpret S −1 as +∞ if S does not have a bounded inverse. This is easier to compute than the spectrum, converges in the Hausdorff metric to the spectrum as ↓ 0 and gives an indication of the instability of the spectrum of T . We shall use it as a comparison for the IQR algorithm and as a means to detect spectral pollution for finite section methods.
Finally, we need a notion of convergence of subspaces. We follow the notation in [41]. Let M ⊂ B and N ⊂ B be two non-trivial closed subspaces of a Banach space B. The distance between them is defined by Given subspaces M and {M k } such thatδ(M k , M) → 0 as k → ∞, we will use the notation M k →M. If we replace B with a Hilbert space H, we can express δ andδ conveniently in terms of projections and operator norms. In particular, if E and F are the orthogonal projections onto subspaces M ⊂ H and N ⊂ H respectively, then This allows us to extend the definition to allow the trivial subspace {0} and gives rise to a metric on the set of all closed subspaces of H (first introduced by Krein and Krasnoselski in [42]). We also define the (maximal) subspace angle, φ(M, N ) ∈ [0, π/2], between M and N by Finally, we will use two further well-known properties in the Hilbert space setting. First, if M and N are both finite l-dimensional subspaces, then 5) which shows that to prove convergence of finite-dimensional subspaces, it is enough to prove δ-convergence. Second, suppose we have where the N (k) j need not be orthogonal. Then a simple application of Hölder's inequality yields For further properties (including other notions of distances between subspaces) and a discussion on two projections theory, we refer the reader to the excellent article of Böttcher and Spitkovsky [19].

Organisation of the paper
The paper is organised as follows. In Sect. 2 we define the IQR algorithm (simple codes are also provided in the appendix). Section 3 contains and proves our main theorems including convergence rates. The outcome is more elaborate than the finitedimensional case, as the infinite-dimensional setting includes more intricate instances.
Our key practical result is that, despite being an algorithm dealing with infinite amount of information, it can be implemented on any standard computer and this is discussed in Sect. 4. The fact that the IQR algorithm can be computed allows for its use in order to provide new classification in the SCI hierarchy as discussed in Sect. 5. In particular, we demonstrate 1 classification for the extremal part of the spectrum and dominant invariant subspaces, as well as 1 results for spectra of certain classes of compact operators. Note that the general spectral problem for compact operators is not in 1 . The IQR algorithm and convergence theorems are demonstrated on a large collection of examples from the sciences on difficult computational spectral problems in Sect. 6, with comparisons to the finite section method. The IQR algorithm is also found to perform better than theory predicts and we conjecture conditions on the operator for this to be the case. Finally, we conclude with a discussion of the opportunities and limits of the IQR algorithm in Sect. 7.

The infinite-dimensional QR algorithm (IQR)
The IQR algorithm has existed as a pure mathematical concept for more than thirty years and it first appeared in the paper "Toda Flows with Infinitely Many Variables" [28] in 1985. However, the analysis in [28] covers only self-adjoint infinite matrices with real entries, and since the analysis is done from a pure mathematical perspective, the question regarding the actual numerical algorithm is left out. We will extend the analysis to more general operators and answer the crucial question: can one actually implement the IQR algorithm? The answer is affirmative, and we also prove convergence theorems, generalising the well-known finite-dimensional case.

The QR decomposition
The QR decomposition is the core of the QR algorithm. If T ∈ C n×n , one may apply the Gram-Schmidt procedure to the columns of T and store these columns in a matrix Q. This gives us the QR decomposition where Q is a unitary matrix and R upper triangular. It is no surprise that a QR decomposition should exist in the infinite-dimensional case, however, we need more than just the existence. A key ingredient in the QR algorithm are Householder transformations, used for computational reasons (they are backwards stable). It is crucial that we can adopt these tools in the infinite-dimensional setting. Our goal is to extend the construction of the QR decomposition, via Householder transformations, to infi-nite matrices and to find a way so that one can implement the procedure on a finite machine. To do this, we need to introduce the concept of Householder reflections in the infinite-dimensional setting.

Definition 2.1 A Householder reflection is an operator S ∈ B(H) of the form
whereξ denotes the associated functional in H * given by x → x, ξ . In the case where H = H 1 ⊕ H 2 and I i is the identity on H i then will be called a Householder transformation.
A straightforward calculation shows that S * = S −1 = S and thus also U * = U −1 = U . An important property of the operator S is that if {e j } is an orthonormal basis for H and η ∈ H, then one can choose ξ ∈ H such that In other words, one can introduce zeros in the column below the diagonal entry. Indeed, if η 1 = η, e 1 = 0 one may choose ξ = η ± η ζ, where ζ = η 1 /|η 1 |e 1 and if η 1 = 0 choose ξ = η ± η e 1 . The following theorem gives the existence of a QR decomposition, even in the case where the operator is not invertible.
This is known as the QR algorithm and is completely analogous to the finitedimensional case. Note also that we have T n =Q * n TQ n . In the finite-dimensional case and under favourable conditions,Q * n TQ n converges to a diagonal operator and the columns ofQ n converge to the corresponding eigenvectors as n → ∞ (see Theorem 3.1 below). We will see that the IQR algorithm behaves similarly for the extreme parts of the spectrum.

Remark 2.4
Note that since the Householder transformations used in the proof of Theorem 2.2 are unique up to a ± sign, we will with some abuse of language refer to the QR decomposition constructed as the QR decomposition. In general for an invertible operator, the IQR algorithm is uniquely defined up to phase-see Sect. 4.2. This will not be a problem for our theorems or numerical examples.
The following observation will be useful in the later developments. From the construction in (2.3) and (2.4) we get An easy induction gives us that T m =Q mRm . (2.5) Note thatR m must be upper triangular with respect to From this it follows immediately that

Convergence theorems
In finite dimensions we have the following well-known theorem: λ j e j ⊗ e j , as m → ∞.
In this section we will address the convergence of the IQR algorithm for normal operators under similar assumptions and prove an analogue of Theorem 3.1 in infinite dimensions (Theorem 3.9). As well as this, and for more general operators T that are not necessarily normal, we address block convergence (Theorem 3.13), relevant when the eigenvalues do not have distinct moduli, and convergence to (dominant) invariant subspaces (Theorem 3.15).

Preliminary definitions and results
To state and prove our theorems we need some preliminary results. The reader only interested in the results themselves is referred to Sect. 3.2. If T is a normal operator, we will use χ S (T ) to denote the indicator function of the set S defined via the functional calculus. Without loss of generality, we deal with the Hilbert space H = l 2 (N) and the canonical orthonormal basis {e j } j∈N . Our first set of results concerns the convergence of spanning sets under power iterations and is analogous to the finite-dimensional case. The following proposition can be found in [36] and together with Lemma 3.6 below, these are the only results we will use from [36].

Proposition 3.2 Suppose that T ∈ B(H) is normal, is invertible and that
consists of finitely many isolated eigenvalues of T with |λ 1 | > |λ 2 | > · · · > |λ N |. Suppose further that sup{|z| : z ∈ } < |λ N |. Let l ∈ N and suppose that {ξ i } l i=1 are linearly independent vectors in H such that {χ ω (T )ξ i } l i=1 are also linearly independent. Then are linearly independent and there exists an ldimensional subspace B ⊂ ranχ ω (T ) such that where D is an (l − 1)-dimensional subspace, then where ξ ∈ ranχ ω (T ) is an eigenvector of T .
In order to extend this proposition to describe rates of convergence and prove our main theorems, we need to describe the space B in more detail. This is done inductively as follows. The first step is to choose ν 1,1 ∈ {λ i } N i=1 of maximum modulus such that We then let ξ 1,1 be a linear multiple of ξ 1 such that χ ν 1 We seek to add the space spanned by the vector ξ m+1 whilst preserving these properties.
First we deal with (2).
After l steps the above process terminates giving a new basis The subspace B can then be described as

Definition 3.3
With respect to the above construction we define the following: Since the Gram-Schmidt process is defined uniquely up to phases we see that are linearly independent then We can now prove the following refinement of Proposition 3.2: Since the spaces are ldimensional, it follows from (1.5) that we have the convergence ratê Proof Consider the subspaces By construction, we have for any suchξ i in the above sum that This gives T kξ i = ν k j χ ν j (T )ξ j,i + T k χ θ j (T )ξ i . Now, by the assumption on σ (T ), we have Thus, since Here we have used Hölder's inequality together with the fact that . The right-hand side gives an upper bound for δ(E j , E k j ). Analogous rates of convergence hold for the other subspaces and from (1.6) we have since the spaces E j are orthogonal.
To apply Propositions 3.2 and 3.4 to prove the main result Theorem 3.9, we need to take care of the case that some of the e j may have χ ω (T )e j = 0. Definition 3.5 Suppose that (A1) and (A2) hold and let K ∈ N ∪ {∞} be minimal with the property that dim(span{χ ω (T )e j } K j=1 ) = M. Define Define also the corresponding subset {ê j } M j=1 ⊂ {e j } K j=1 such that {ê j } M j=1 = ω \˜ ω and such that upon writingê j = e p j , the p j are increasing.
Note that we have the following decomposition of T into is an orthonormal set of eigenvectors of T . The following simple lemma extends Lemma 39 in [36] to infinite M but the proof is verbatim so omitted.

Lemma 3.6 If e m ∈
∪˜ ω , then The following theorem is the key step of the proof of Theorem 3.9 and concerns convergence to the eigenvectors of T . ∪˜ ω and μ is maximal with p μ < m (recall thatê j = e p j ), then we have In the case that m < p 1 , we interpret this as χ ω (T )q k,m = 0 which holds from Lemma 3.6. (b) For any j < M + 1, and henceδ Here, as in Lemma 3.6, q k, j = Q k e j andq k, j = Q kê j . Finally, if M is finite then we must have span{q j } M j=1 = ranχ ω (T ).
We will provide an inductive proof of Theorem 3.7 which requires the following for the inductive step of part (a).
Proof First note that from (2.6), invertibility of T and the fact that {χ ω (T )ê j } μ j=1 are linearly independent, it must hold that {χ ω (T )q k, j } μ j=1 are linearly independent also. Then by using the assumptions stated and the fact that χ ω (T )q j =q j we have Also, we have that s(m) ≤ μ and Lemma 3.6 implies Using the fact that χ ω (T )q k,m ≤ 1 and the definition of δ (along with the fact that span{q j } μ j=1 is finite-dimensional), it follows that there exists some v k = μ j=1 β j,kq j ∈ span{q j } μ j=1 with v k ≤ 1 and We also have from assumption (b) that since q k,m is orthogonal toq k, j . This together with (3.7) gives that β j,k ≤ C(μ) Using (3.7) again then gives the result. Note that we have used orthonormality of {q j } μ j=1 which will be proven as part of the induction.
Proof of Theorem 3. 7 We begin with the initial step of the induction for (b) and (c). Note that (a) trivially holds by construction with A(m) = 0 for any m < p 1 where e p 1 =ê 1 and this provides the initial step for (a). By Propositions 3.2 and 3.4, there exists a unit eigenvectorq 1 ∈ ranχ ω (T ) such that Thus, it follows that are orthonormal (recall that Q k is unitary) and hence by (3.9) there exists some coefficients α k,i with If e m ∈ ∪˜ ω , where m < p 1 then by Lemma 3.6 q k,m ,q 1 = 0. It follows that we must have Hence we can take B(1) = 1 and C(1) = 1 in (b) and (c) respectively which completes the initial step.
For the induction step we will argue simultaneously for (a), (b) and (c) using induction on μ. Suppose that (a) holds for m < p μ with e p μ =ê μ together with (b) and (c) for j ≤ μ and some μ < M. Let e p μ+1 =ê μ+1 then we can use Lemma 3.8 to extend (a) to all m < p μ+1 and this provides the step for (a). For (b), we note that Propositions 3.2 and 3.4 imply that (3.11) where ξ is a unit eigenvector of T . We may also assume without loss of generality that ξ is orthogonal toq j for j = 1, . . . , μ. As before, since span{T kê and hence by invertibility of T If e m ∈ ∪˜ ω , where m < p μ+1 then as shown above we have Taking the inner product of ξ −η k with q k,m and using (3.13) together with the orthonormality of the q k, j s, it follows that α k,m ≤ A(m) To finish the inductive step, we defineq μ+1 = ξ . Recall that ξ is orthogonal to anyq l with l ≤ μ. Hence it follows that {q i } μ+1 i=1 are orthonormal and we can take are orthonormal and (1.6) imply we can take

Main results
Our first result generalises Theorem 3.1 to infinite dimensions and relies on Theorem 3.7 (which concerns convergence to eigenvectors).

) is a collection of orthonormal eigenvectors of T and if
where P M denotes the orthogonal projection onto span{ê j } M j=1 . Furthermore, if we define then r < 1 and for any fixed x ∈ span{ê j } M j=1 we have the following rate of convergence (3.14) If M is finite then we can write (after possibly re-ordering)

15)
and in part (ii) we have the rate of convergence If {χ ω (T )e l } M l=1 are linearly independent, then we can takeê j = e j . Remark 3.10 What Theorem 3.9 essentially says is that if we take the n-th iteration of the IQR algorithm and truncate to an m × m matrix (i.e. P m Q * n T Q n P m ) then, as n grows, the eigenvalues of this matrix will converge to the extremal parts of the spectrum of T . In particular, the theorem suggests that the IQR algorithm can locate the extremal parts of the spectrum.

Proof of Theorem 3.9
To prove (i), since a closed ball in B(l 2 (N)) is weakly sequentially compact, it follows that any subsequence of {Q * n T Q n } n∈N must have a weakly convergent subsequence {Q * n k T Q n k } k∈N . In particular, there exists a W ∈ B(l 2 (N)) such that Let P M denote the projection onto span{ê j } M j=1 . Note that part (i) of the theorem will follow if we can show that Tq j ,q j ê j ⊗ê j , (3.17) and We will indeed show this, and we start by observing that, due to the weak convergence and the standard functional calculus, we have that We then have the following Thus, by (3.18), (3.20), (3.21) and Theorem 3.7 we get (3.17) and also that P ⊥ M W P M = 0. Also, by (3.19), (3.20), (3.21) and Theorem 3.7 we get that P M W P ⊥ M = 0. Note that in all of these cases, Theorem 3.7 implies that the rate of convergence is such that the difference between Wê j , e i , W e i ,ê j and their limiting values is O(r n k ) (however, not necessarily uniformly over the indices). Now suppose that T has finitely many non-zero entries in each column. This can be described by a function f : N → N nondecreasing with f (n) ≥ n such that T e j , e i = 0 when i > f ( j) as in Definition 4.1. Proposition 4.2 shows that this is preserved under the iteration in the IQR algorithm, i.e. Q * n k T Q n k also has this property. So let x ∈ l 2 (N) and > 0. Choose y of finite support such that x − y ≤ . It is then clear that Q * n k T Q n k y − W y → 0 as n k → ∞ (since we only require convergence in finitely many entries). Hence lim sup Since > 0 and x were arbitrary, we have Q * To prove (ii), suppose that x ∈ span{ê j } M j=1 , then x can be written as with at most finitely many x j non-zero. We have thatδ(span{Q nê j }, span{q j }) = O(r n ) and hence there exists some a n, j of unit modulus such that Q nê j − a n, jq j = O(r n ). Since Q n is unitary, we then have where we have used the fact that T is bounded in the last line. We therefore have convergence on span{ê j } M j=1 , and, since the operators are uniformly bounded, we must have convergence on span{ê j } M j=1 which implies that For the last parts, suppose that M is finite. Theorem 3.7 then implies (3.15) after a possible re-ordering. The rate of convergence in (3.14) also implies that Recall also from the proof of Theorem 3.7 that are linearly independent then˜ ω = ∅, and therefore {ê j } M j=1 = {e j } M j=1 , which yields that the projection P M in (3.17) is the projection onto span{e j } M j=1 . Theorems 3.9 and 3.7 also give us convergence to the eigenvectors. With the use of (possibly countably many) shifts and rotations, the above theorem allows us to find all eigenvalues, their multiplicities and eigenspaces outside the convex hull of the essential spectrum, i.e. outside the essential numerical range.

Example 3.11
It is possible in the case of infinite M that theq j do not form an orthonormal basis of ranχ ω (T ) and we can even lose part of ω in the convergence of P M Q * n T Q n P M to a diagonal operator. This is to be contrasted to the finite-dimensional case. For example, suppose that with respect to an initial orthonormal basis {v j } j∈N , T is given by the diagonal matrix Diag(1/2, 1, 1, . . .). Now define f j = v 1 + (1/ j)v j+1 and apply Gram-Schmidt to the sequence { f j } j∈N to generate orthonormal vectors {e j } j∈N . It is easy to see that any v j can be approximated to arbitrary accuracy using finite linear combinations of e j and hence {e j } j∈N is an orthonormal basis of our Hilbert space. We also have that the χ 1 (T )( f j ) = (1/ j)v j+1 are linearly independent and hence so are χ 1 (T )(e j ). It follows that the IQR iterates converge in the strong operator topology to the identity operator. However, we could equally take ω = {1, 1/2} in Theorem 3.9. Hence we have the curious case that span{q j } j∈N ⊂ span{v j } j>1 and we lose the eigenvalue 1/2.
The following corollary is entirely analogous to the finite-dimensional case.
are linearly independent. In the notation of Theorem 3.9, let ρ = sup{|z| : z ∈ }. For j < N define r j = max{|λ k+1 /λ k | : k ≤ j} and for j = N define r N = max{|λ k+1 /λ k |, |λ N /ρ| : k ≤ j}. We then have the following rates of convergence to the diagonal operator for i, j ≤ M: Proof The result follows from Theorem 3.9 applied successively to ω 1 , ω 2 , . . . , ω N where ω j = {λ k : k ≤ j}. In general, analogous results follows from Theorem 3.9 when M is infinite and with other linear independence conditions on χ ω (T )e i with ω ⊂ ω but the statements become less succinct.
In the finite-dimensional case and the case of distinct eigenvalues of the same magnitude, the QR algorithm applied to a normal matrix will 'converge' to a block diagonal matrix (without necessarily converging in each block). This can be extended to infinite dimensions by inductively using the following theorem which also extends to non-normal operators.

if T has only finitely many non-zero entries in each column then we can replace W O T convergence by S O T convergence.
Remark 3.14 Theorem 3.13 essentially says that the IQR algorithm can compute the invariant subspace ran(P) of such an operator if there is enough separation between T restricted to ran(P) and ran(I − P). In other words, provided the existence of a dominant invariant subspace.

Proof of Theorem 3.13
The main ideas of the proof of Theorem 3.13 have already been presented so we sketch the proof. We first define the vectors {ê j } M j=1 in a similar way to Definition 3.5 inductively byê j = e p j where Let r = β/α < 1. We will prove inductively that for some constants C 1 (μ) and C 2 ( j). Suppose that this has been done. Part (i) of Theorem 3.13 now follows since span{P Q nê j } μ j=1 ⊂ ran(P). We then argue as in the proof of Theorem 3.9 to gain Then by studying the inner products T Q n k e j , Q n k e i using the invariance of ran(P), ran(I − P) under T and from (b), part (ii) of Theorem 3.13 easily follows (note that (a) implies that (I − P)Q nê j ≤ C 1 ( j)r n ). The final part of the theorem then follows from the same arguments in the proof of Theorem 3.9. Hence we only need to prove (a) and (b). We first claim that P commutes with T which is invertible and hence both of these spaces have dimension μ by the construction of theê j . It follows that (3.22) implieŝ To show (3.22), let x n 1 , . . . , x n μ be an orthonormal basis for span{PT nê j } μ j=1 and let ξ = μ j=1 α j x n j have norm at most 1. Now, we may choose coefficients β j,n such that T n μ j=1 β j,n x n j = ξ since T | ran(P) is invertible when viewed as an operator acting on ran(P). By the assumptions on T we must have that We may change basis from Then clearly by Hölder's inequality proving (3.22) and hence (3.23).
Note that the proof of Lemma 3.6 carries over (replacing the projection χ ω (T ) by P) to prove that where we have used (2.6) to reach the second line and the fact that span{PT nê j } μ j=1 ⊂ span{PT n e j } p μ j=1 to reach the third line. Again, both spaces have dimension μ so we have With these arguments out of the way (these are the analogue of Proposition 3.4) we can now form our inductive argument, similar to the proof of Theorem 3.7. Suppose first that (a) holds for μ (allowing μ = 0 for the initial step) and let j ∈ have Then we must have a n,i + v n , Q nêi = P Q n e j , Q nêi = Q n e j , P Q nêi .
Using (a) again, along with the fact that Q n e j is orthogonal to {Q nêi } μ i=1 , we must have a n,i ≤ 2C 1 (μ)r n . It follows that we can take have unit norm and assume that p μ+1 < ∞ (else there is nothing to prove since then μ = M). Then there exists b n, j and w n such that ξ = p μ+1 j=1 b n, j Q n e j + w n and w n ≤ C 5 (μ + 1)r n . Now let j ∈ with j < p μ+1 then we must have ξ, P Q n e j = ξ, Q n e j = b n, j + w n , Q n e j .
We have proven (b) for such j and hence we have b n, j ≤ C 2 ( j) + C 5 (μ + 1) r n . It follows that we can take where the square root factor appears since the relevant spaces are μ-dimensional. This completes the inductive step (the initial step is identical) and hence the proof of the theorem.
Theorem 3.13 can be made sharper (under a slightly stricter assumption on the linear independence of {e j } M j=1 ) with the following theorem which includes the case that ran(I − P) is not necessarily invariant. N)) be an invertible operator (not necessarily normal) and suppose that there exists an orthogonal projection P of finite rank M such that the range of P is invariant under T . Suppose also that there exists α > β > 0 such that Under these conditions, there exists a canonical M-dimensional T * −invariant subspace S and we letP denote the orthogonal projection onto S (in the special case that ran(I − P) is also T -invariant such as in Theorems 3.9 and 3.13, then S = ran(P)). Suppose also that {Pe j } M j=1 are linearly independent. Let {Q n } n∈N and {R n } n∈N be Q-and R-sequences of T with respect to {e i }. Then

Furthermore, if T has only finitely many non-zero entries in each column then we can replace W O T convergence by S O T convergence.
Remark 3.16 Theorem 3.15 says that the IQR algorithm can be used to approximate dominant invariant subspaces. In particular, we shall use the bound (3.26) to build a 1 algorithm in Sect. 5. Note in the normal case that Theorem 3.9 is more precise, both in giving convergence of individual vectors to eigenvectors and in the less restrictive assumptions on spanning sets and M. In the normal case (and that of Theorem 3.13) we also have that the limit operator has a block diagonal form.

Proof of Theorem 3.15
In this section we will prove Theorem 3.15. The proof technique is different from those used above, and hence we have given it a separate section. Throughout, we will denote the ratio β/α by r . Note that since M is finite, the bound α implies that j=1 is a basis for ran(P). Then as matrices with respect to the original basis we can write where T 11 ∈ C M×M and T 12 has M rows. Our assumptions imply that T −1

11
≤ 1/α and T 22 ≤ β. The next lemma shows that we can change the basis further to eliminate the sub-block T 12 . This is needed to apply a power iteration type argument.

It is then straightforward to check
The canonical T * −invariant subspace alluded to in Theorem 3.15 is then simply S = span{Y e j } M j=1 . The space is canonical since it is easily seen that it is unchanged if we use a different basis for ran(P 1 ) and ran(P 2 ) in the definition of Q. Now let P 0 = e 1 e 2 . . . e M ∈ B(C M , l 2 (N)) denote the matrix whose columns are the first M basis elements {e j } M j=1 . Since the {R i } are upper triangular, it is easy to see that T n P 0 = Q n R n P 0 = Q n P 0 P * 0 R n P 0 .
We will denote the (invertible) matrix P * 0 R n P 0 ∈ C M×M by Z n . Now define then we have the relation Unwinding the definitions, this implies the matrix identities Lemma 3. 18 The following identity holdŝ Proof Note that span{Q n e j } M j=1 = ran(Q n P 0 ) and ran(P) = ran(P 1 ). Since P 1 P * 1 and Q n P 0 P * 0 Q * n are orthogonal projections, it follows that δ(span{Q n e j } M j=1 , ran(P)) = Q n P 0 P * 0 Q * n − P 1 P * But we have that P * 0 Q * n P 2 = V 2 n and hence we are done if we can show P * 0 Q * n P 2 = (I − P 0 ) * Q * n P 1 . Consider the unitary matrix Now let x ∈ C M be of unit norm, then U 11 x 2 + U 21 x 2 = 1. It follows that where σ 0 denotes the smallest singular value. Applying the same argument to U * we see that U 12 2 = 1 − σ 0 (U 11 ) 2 = U 21 2 , completing the proof. .

(3.31)
Proof First note that since {Pe j } M j=1 are linearly independent, we must have φ(span{e j } M j=1 , S) < π/2 and hence the bound in (3.31) is finite. l 2 (N)). By considering W * W = I ∈ C M×M , we see that the columns of W are orthonormal. In fact, expanding Y we have and hence the columns of W are a basis for the subspace S. Arguing as in the proof of Lemma 3.18, we have that This implies that W * P 0 is invertible with We also have the identity Since (I + A A * ) −1/2 has norm at most 1, we see that (V 1 0 − AV 2 0 ) is invertible and (3.31) holds.

Proof of Theorem 3.15
Using Lemma 3.19 and the matrix identities (3.28) and (3.29), we can write Using (3.30) and (3.31), this implieŝ δ(span{Q n e j } M j=1 , ran(P)) ≤ . (3.32) It is clear by summing a geometric series that . Substituting this into (3.32) proves part (i) of the theorem.
Next we argue that if i > M then P Q n e i → 0 as n → ∞. We have that P Q n e i = M j=1 α j,n Q n e j + v n with v n → 0 by part (i). Note that we then have α j,n = P Q n e i , Q n e j + j,n = Q n e i , P Q n e j + j,n with { j,n } null. But again by (i) we have that P Q n e j approaches span{Q n e k } M k=1 which is orthogonal to Q n e i and hence {α j,n } is null. The proof of part (ii) now follows the same argument as in the proof of part (i) of Theorem 3.9 and of the final part of Theorem 3.13. The key property being that if j ≤ M and i > M then Q * n T Q n e j , e i → 0 due to the invariance of ran(P) under T . Note that it does not necessarily follow (as is easily seen by considering upper triangular T ) that Q * n T Q n e i , e j → 0 for such i, j.

The IQR algorithm can be computed
The previous section gives a theoretical justification for why the IQR algorithm may work, but we are faced with the possibly unpleasant problem of how to compute with infinite data structures on a computer. Fortunately, there is a way to overcome such a problem. The key is to impose some structural requirements on the infinite matrix. This is the class of infinite matrices with a finite number of non-zero elements in each column (and not necessarily in each row) which is captured by the function f . It is for this class that the computation of the IQR algorithm is feasible on a finite machine. For this class of operators one can actually compute (without any approximation or any extra discretisation) the matrix elements of the n-th iteration of the IQR algorithm as if it was done on an infinite computer (meaning the computation collapses to a finite one). The following result of independent interest is needed in the proof and generalises the well-known fact in finite dimensions that the QR algorithm preserves bandwidth (see [57] for a good discussion of the tridiagonal case). T ∈ B(l 2 (N)) and let T n be the n-th element in the IQR iteration, such that T n = Q * n · · · Q * 1 T Q 1 · · · Q n , where (4.1)

Proposition 4.2 Let
Using the fact that f is increasing, it follows that each U 1 m has quasi-banded subdiagonals with respect to f , as does the product U 1 1 · · · U 1 m . It follows that Q 1 must have quasi-banded subdiagonals with respect to f and hence so does T 1 = R 1 Q 1 since R 1 is upper triangular. (4.2)

Remark 4.4
What Theorem 4.3 says is that to compute the finite section of size m of the n-th iteration of the IQR algorithm (i.e. P m T n P m ), one only needs information from the finite section of size f n (m) (i.e. P f n (m) T P f n (m) ) since the relevant Householder reflections can be computed from this information. In other words, the IQR algorithm can be computed.

Proof of Theorem 4.3 By induction it is enough to prove that
To see why this is true, note that by the assumption that T has quasi-banded subdiagonals with respect to f , Proposition 4.2 shows that T n has quasi-banded subdiagonals with respect to f for all n ∈ N. Thus, it follows from the construction in the proof of Theorem 2.2 that each U j l is of the form where I l, j,1 denotes the identity on P l−1 H, I l, j,2 denotes the identity on span{e k : l ≤ k ≤ f (l)}, I l, j,3 denotes the identity on P ⊥ f (l) H and ξ l, j ∈ span{e k : l ≤ k ≤ f (l)}. Since P m is compact, it then follows that  only affects finitely many columns or rows of A if multiplied either on the left or the right. In computer science, it is often referred to as "Lazy evaluation" when one computes with infinite data structures, but defers the use of the information until needed. A simple implementation is shown in the appendix for the case that the matrix has k subdiagonals (i.e. we have f (n) = n + k).
The next question is how restrictive is the assumption in Definition 4.1? In particular, suppose that T ∈ B(H) and that ξ ∈ H is a cyclic vector for T (i.e. span{ξ, T ξ, T 2 ξ, . . .} is dense in H). Then by applying the Gram-Schmidt procedure to {ξ, T ξ, T 2 ξ, . . .} we obtain an orthonormal basis {η 1 , η 2 , η 3 , . . .} for H such that the matrix representation of T with respect to {η 1 , η 2 , η 3 , . . .} is upper Hessenberg, and thus the matrix representation has only one subdiagonal. The question is therefore about the existence of a cyclic vector. Note that if T does not have invariant subspaces then every vector ξ ∈ B(H) is a cyclic vector. Now what happens if ξ is not cyclic for T ? We may still form {η 1 , η 2 , η 3 , . . .} as above, however, H 1 = span{η 1 , η 2 , η 3 , . . .} is now an invariant subspace for T and H 1 = H. We may still form a matrix representation of T with respect to {η 1 , η 2 , η 3 , . . .}, but this will now be a matrix representation of T | H 1 . Obviously, we can have that σ (T | H 1 ) σ (T ).
However, the following example shows that the class of matrices for which we can compute the IQR algorithm covers a wide number of applications. In particular, it includes all finite interaction Hamiltonians on graphs. Such operators play a prominent role in solid state physics [48,51] describing propagation of waves and spin waves as well as encompassing Jacobi operators studied in many physical models and integrable lattices [71].

Example 4.6
Consider a connected, undirected graph G, such that each vertex degree is finite and the set of vertices V (G) is countably infinite. Consider the set of all bounded operators A on l 2 (V (G)) ∼ = l 2 (N) such that the set S(v) := {w ∈ V : w, Av = 0} is finite for any v ∈ V . Suppose our enumeration of the vertices obeys the following pattern. e 1 's neighbours (including itself) are S 1 = {e 1 , e 2 , . . . , e q 1 } for some finite q 1 . The set of neighbours of these vertices is S 2 = {e 1 , . . . , e q 2 } for some finite q 2 where we continue the enumeration of S 1 and this process continues inductively enumerating S m . If we know S(v) for all v ∈ V then we can find an f : We simply choose f (n) = q r n where r n is minimal such that ∪ j≤n S(e j ) ⊂ S r n .

Invertible operators
More generally, given an invertible operator T with information on how its columns decay at infinity, we can compute finite sections of the IQR iterates with error control. For computing spectral properties, we can assume, by shifting T → T + λI then translating by −λ back, that the operator we are interested in is invertible, hence the invertibility criterion is not that restrictive. Throughout, we will use the following lemma which says that for invertible operators, the QR decomposition is essentially unique.

Lemma 4.7 Let T be an invertible operator (viewed as a matrix acting on l 2 (N)), then there exists a unique decomposition T = Q R with Q unitary and R invertible,
upper triangular such that R ii ∈ R >0 . Furthermore, any other "QR" decomposition T = Q R has a diagonal matrix D = Diag(t 1 , t 2 , . . .) such that |t i | = 1 and Q = Q D. In other words, the QR decomposition is unique up to phase choices.
Proof Consider the QR decomposition already discussed in this paper, T = Q R . T is invertible, and hence Q is a surjective isometry so is unitary. Hence R = Q * T is invertible. Being upper triangular, it follows that R ii = 0 for all i. Choose t i ∈ T such that t i R ii ∈ R >0 and set D = Diag(t 1 , t 2 , . . .). Letting Q = Q D * and R = D R we clearly have the decomposition as claimed. Now suppose that T = Q R then we can write Q = Q R R −1 . It follows that R R −1 is a unitary upper triangular matrix and hence must be of the form D = Diag(t 1 , t 2 , . . .) with |t i | = 1.
Another way to see this result is to note that the columns of Q are obtained by applying the Gram-Schmidt procedure to the columns of T . The restriction that R ii ∈ R >0 can also be incorporated into Theorem 4.3. Theorem 4.3 (in this subcase of invertibility) is then a consequence of the fact that if T has quasi-banded subdiagonals with respect to f then P m T n P m = P m (P f n (m) T P f n (m) ) n P m and the relations (2.5)-we can apply Gram-Schmidt (or a more stable modified version) to the columns of P f n (m) T P f n (m) and truncate the resulting matrix.
Assume that given T ∈ B(l 2 (N)) invertible (not necessarily with quasi-banded subdiagonals), we can evaluate an increasing family of increasing functions g j : N → N such that defining the matrix T ( j) with columns {P g j (n) T e n } we have that T ( j) is invertible and It is easy to see that such a sequence of functions must exist since any S with S − T < T −1 −1 is invertible. Given this information, without loss of generality by increasing the g j s pointwise if necessary, applying Hölder's inequality and taking subsequences, we may assume that T ( j) − T ≤ 1/ j. In other words, given a sequence of functions satisfying (4.5) we can evaluate a sequence of functions with this stronger condition. The following says that given such a sequence of functions, we can compute the truncations P m T n P m to a given precision.

Proof of Theorem 4.8
First consider the error when applying Theorem 4.3 to T ( j) with g j for any fixed j. We will show that we can compute an error bound which converges to zero as j → ∞ and from this the theorem easily follows by successively computing the bound and halting when this bound is less than .
We have T − T ( j) ≤ 1/ j and hence, by writing whereC = (C + 1) n . The columns ofQ n andQ ( j),n are simply the columns of the matrices T n and (T ( j) ) n after the application of Gram-Schmidt. Let the first m columns of T n and (T ( j) ) n be denoted by {t k } m k=1 and {t j k } m k=1 respectively and let {q k } m k=1 and {q j k } m k=1 be the vectors obtained after applying Gram-Schmidt to these sequences of vectors. We then have For a vector v of unit norm, let P ⊥v denote the orthogonal projection onto the space of vectors perpendicular to v. Note that for two such vectors v, w, we have then q k are just the normalised version of v k and likewiseq j k are just the normalised version ofṽ j k . Suppose that for μ < k we have q μ −q j μ ≤ δ for some δ > 0. Then applying the above products of projections we have In the last line we have used the fact that if the operators {A l } m l=1 and {B l } m l=1 have norm bounded by 1, then Applying the same argument as in the inequalities (4.6) we see that We must have q k −q j k ≤ δ m ( j) for 1 ≤ k ≤ m where we have now shown the j dependence as an argument.
It follows that (Q n −Q ( j),n )P m ≤ √ mδ m ( j) and hence that P m T n P m − P m T ( j),n P m ≤ P m (Q n −Q ( j),n ) * TQ n P m So we need only show that δ m ( j) → 0 as j → ∞. Note that as j → ∞, the columns of (T ( j) ) n converge to that of T n . It follows thatt j k converge to t k andq j 1 converges to q 1 . An easy inductive argument using (4.7) and (4.8) shows that the vectorsq j k converge to q k and ṽ j k are bounded below. The convergence δ m ( j) → 0 now follows.

SCI classification theorems
In this section we will apply the above results to prove three new classification theorems in the SCI hierarchy. First, assume that T ∈ B(l 2 (N)) is an invertible normal operator with σ (T ) = ω ∪ , where ω ∩ = ∅, ω = {λ i } N i=1 , and the λ i 's are isolated eigenvalues with multiplicity m i satisfying |λ 1 | > · · · > |λ N |. As usual, we also assume that sup{|θ | : θ ∈ } < |λ N | and set M := m 1 + · · · + m N ∈ N ∪ {∞}. (5.1) In this section we will assume for simplicity that all the m i except possibly m N are finite. To be able to obtain the classification results we need two key assumptions.
(I) Column decay We assume a much weaker condition than bandedness of the infinite matrix. Indeed, we suppose a known decay of the elements in the columns of T that is described through a family of increasing functions {g j } j∈N . In particular, g j : N → N is such that defining the infinite matrix T ( j) with columns {P g j (n) T e n } n∈N we have that T ( j) is invertible and (II) Distance to span of eigenvectors In order to obtain error control ( 1 classification) one needs to control the hidden constant in the O(r n ) estimate in (3.16). This is done as follows, where {Q n } n∈N is a Q-sequence of T with respect to {e j } j∈N . Given finite k < M + 1 with m 1 + · · · + m N −1 < k, we will assume that if l < N then {χ {λ 1 ,...,λ l } (T )e j } m 1 +···+m l j=1 are linearly independent. We also assume that {χ {λ 1 ,...,λ N } (T )e j } k j=1 are linearly independent. This simply ensures that the IQR algorithm converges with the expected ordering (largest eigenvalue in the first diagonal entry then in descending order). It follows from Theorems 3.9 and 3.7, that there exist eigenspaces E 1 , . . . , E N (with the last space depending on k and the vectors {e j }) corresponding to the eigenvalues λ 1 , . . . , λ N such that We then define the initial supremum subspace angle by where φ, defined by (1.4), denotes the subspace angle. Our assumptions and the proofs in Sect. 3 show that (T , {e j } k j=1 ) < π/2 and hence the key quantity tan (T , {e j } k j=1 ) is finite.

Remark 5.1
The quantity tan (T , {e j } k j=1 ) can be viewed as a measure of how far {e j } k j=1 is from {q j } k j=1 , the k eigenvectors of T corresponding to the first k eigenvalues (including multiplicity and preserving order). Hence it gives an estimate of how good the initial approximation {e j } k j=1 to {q j } k j=1 is. Indeed, we know from (3.16) that the convergence rate is O(r n ), and the hidden constant C depends exactly on this behaviour. In particular, if e j = q j for j ≤ k then C = 0.
Define also We can now define the class of operators k t,L for the classification theorem.
As discussed in Remark A.6 in an Appendix, where we review the SCI hierarchy, when we speak of convergence of M n (T ) to 1 (T ), we define, with a slight abuse of notation, Having established the basic definition we can now present the classification theorem.

Remark 5.4
Note that this means we converge to the k largest magnitude eigenvalues in order with error control, and not just arbitrary points of the spectrum. This is in contrast to most 1 classifications in the SCI hierarchy where the best we can hope for is to bound dist(z, σ (T )) for z ∈ C.

Proof of Theorem 5.3
Let T ∈ k t,L then by the definition of k t,L , we may takeê j = e j for j = 1, . . . , k in the arguments in Sect. 3.1. The first step is to bound Z (T , {e j } k j=1 ) in terms of (T , {e j } k j=1 ). Let {ẽ j } k j=1 denote the basis described in Sect. 3.1. In our case: Let δ j = ẽ j then we must have that if m 1 + · · · , m l−1 < j ≤ m 1 + · · · + m l then Where the first line holds since the nearest point toẽ j in is simply χ λ l (T )ẽ j and the E i are defined as above and in (3.1). Rearranging, this implies that .
Hence it follows that In particular, Theorem 3.7 and its proof now implies that where {q j } k j=1 are orthonormal eigenvectors of T and Q m is a Q−sequence of T . In particular, {B( j)} k j=1 can be computed in finitely many arithmetic operations from the induction proof of Theorem 3.7. It follows that there exists z j,m ∈ C of unit modulus such that defining β = max{B(1), . . . , Note that we do not need to assume knowledge of N for this bound (trivially N ≤ k). Using that Q m is an isometry, this implies that where Tq j = λ a j . Note that we must have {λ a j } m 1 +···+m l j=m 1 +···+m l−1 +1 = λ l and {λ a j } k j=m 1 +···+m N −1 +1 = λ N by 3. in the definition of k t,L .
Given any > 0, choose m large enough so that 2Lβt m ≤ and βt m ≤ . The fact that T ≤ L and (5.2) hold implies that we can compute Q * m T Q m e j , e j to accuracy using finitely many arithmetical and square root operations using Theorem 4.8. Call these approximationsλ 1 ,λ 2 , . . . ,λ k . Furthermore, the proof of Theorem 4.8 also makes clear that we can compute Q m e j ∈ l 2 (N) to accuracy using finitely many arithmetical and square root operations (the approximations have finite support). Call these approximationsq 1 ,q 2 , . . . ,q k . Then set (λ 1 ,λ 2 , . . . ,λ k ) × (q 1 ,q 2 , . . . ,q k ).
Next, suppose that we have a continuous increasing function g : R ≥0 → R ≥0 diverging at ∞ such that g(0) = 0 and g(x) ≤ x. Let g IQR be the set of all operators T acting on l 2 (N) (i.e. we fix the representation w.r.t. the canonical basis) for which the IQR algorithm converges in the weak operator topology to a diagonal matrix with the same spectrum as T and such that Note that by Theorem 3.9 this includes all normal compact operators, T , such that {z ∈ σ (T ) : |z| = s} has size at most 1 for all s > 0 (where we can take g(x) = x). 2 We will allow evaluations of g in our algorithms and also assume that we are given functions that satisfy (5.2) and have an upper bound for T . We consider computing 2 (T ) = σ (T ) in the space of compact non-empty subsets of C with the Hausdorff metric.

It is also clear that n (T ) → σ (T ) in the Hausdorff metric.
The final result considers dominant invariant subspaces discussed in Theorem 3.15. Let M ∈ N, t ∈ (0, 1) and L > 0. We let˜ M t,L denote the class of operators such that the assumptions of Theorem 3.15 hold (same M) and such that: We also assume that we are given functions that satisfy (5.2) and consider computing the dominant invariant subspace 3 (T ) = ran(P) in the space of M-dimensional subspaces of l 2 (N) equipped with the metricδ.

Examples and numerical simulations
The aim of the is section is threefold: 1. To demonstrate the convergence and implementation results of Sects. 3-5 on practical examples. 2. To demonstrate that, as well as the proven results, the IQR algorithm performs better than theoretically expected in many cases. In particular, we conjecture that for normal operators whose essential spectrum has exactly one extremal point, the IQR algorithm will also converge to this point. We also demonstrate cases where this seems to hold even if there are multiple extreme points of the essential spectrum and even in non-normal cases. 3. To compare the IQR algorithm to the finite section method and show that in some cases it considerably outperforms it. In general, one can view σ (P m Q * n T Q n | P m H ) as a generalised version of the finite section method, now with two parameters (m and n) that can be varied with n controlling the number of IQR iterates. In some cases we find this avoids spectral pollution whilst still converging to the entire spectrum.
Before embarking with some numerical examples, two remarks are in order. First, extra care has been taken in the case of non-self-adjoint operators whose finite truncations can be non-normal, and hence the computation of their spectra can be numerically unstable. Unless stated otherwise, all calculations were performed in double precision (in MATLAB) and have been checked against extended precision [32] to ensure that none of the results are due to numerical artefacts. Second, when dealing with operators acting on l 2 (Z) we use N as an index set by listing the canonical basis as e 0 , e 1 , e −1 , e 2 , e −2 , . . ., allowing us to apply the IQR algorithm on l 2 (N). Of course different indexing is possible and in general, this would lead to different implementations of the IQR algorithm, 3 but we stick with this ordering throughout.

The finite section method
We first briefly say a few words on the finite section method, the standard means to discretise infinite matrices, since comparisons will be made later. If {P m } m∈N is a sequence of finite-rank projections such that P m+1 ≥ P m and P m → I strongly, where I is the identity, then the idea is to replace T by the finite square matrix P m T | P m H (typically, one takes P m to be the orthogonal projection onto span{e 1 , . . . , e m }). Thus, to find σ (T ), we instead compute σ (P m T | P m H ). However, there can be significant issues when using the finite section method. In general, there is no guarantee that the computed spectra σ (P m T | P m H ) need converge to σ (T ).
For example, consider the shift operator Se j = e j+1 on l 2 (N). If P m projects onto span{e 1 , . . . , e m }, we would get that σ ((P m S| P m H ) = {0} for all m, whereas σ (S) is the closed unit disc. We can also have that σ (P m T | P m H ) σ (T ). For example, let where a j = 5 cos( j)/4 + 2i sin( j). To gain an accurate picture of the spectrum, note that A is banded and hence we can compute approximates to the pseudospectrum [38]. In order to approximate the spectrum in the best possible way we must take as small as possible. Unfortunately, there is a restriction to how small can be depending on mach (machine precision) of the software used. To illustrate this, observe that the approximates are given by (a discrete version of) Thus, ignoring the additional error in computing the smallest eigenvalues of the squared operator and assuming A to have matrix entries of order 1, computing m (A) will have the same challenges as if one squares a real number and then takes its square root. In particular, due to the floating point arithmetic used in the software and (6.2), we must ensure that and this puts a serious restriction on our computation, particularly for the non-normal case where the distance d H (σ (T ), σ (T )) may be large (though we always have σ (T ) ⊂ σ (T )). However, it is possible to detect spectral pollution outside of σ (T ) if we can approximate it well. The phenomenon of "spectral pollution" occurs for A: namely, the computed spectrum σ (P m A| P m H ) contains elements that have nothing to do with σ (A). This is  Fig. 1, an example with spectral pollution z / ∈ σ 1/10 (A), where the same phenomenon occurs for larger m. The spectral pollution phenomenon is well-known. As the following theorem suggests, such pollution can be arbitrarily bad. Despite this result, the finite section can perform quite well. This is the case for self-adjoint operators [6,21,36] and it is also well suited for the computation of pseudospectra of Toeplitz operators [14,18]. Moreover, in general, we have the following (recall that W e (T ) is the convex hull of the essential spectrum for T normal): However, if we want to use the finite section method and rely on Theorem 6.2, we must know W e (T ), and that may be unpleasant to compute. Alternatively, we could hope that σ ess (T ) is close to W e (T ). For example, if T is hypo-normal (T * T − T T * ≥ 0) then conv(σ ess (T )) = W e (T ), where conv(σ ess (T )) denotes the convex hull of σ ess (T ). But what if we have a "very non-normal" operator?
Another problem we may encounter when using the finite section method is that even though σ d (T ) may be recovered, one may get a very misleading picture of the rest of the spectrum. Such problems are illustrated in the following simple example. Let where t j = 1 + 0.5(sin( j) + i cos( j)) for j ≥ 4. This operator decomposes into an upper 4 × 4 block and an operator acting on the perpendicular subspace. It is also possible to compute the spectrum analytically (it consists of a disc of radius 1 centred at 1 together with two isolated eigenvalues). Again, we can compute the pseudospectrum of T (Fig. 1) to reveal that whilst the eigenvalues produced by the finite section method are correct, they do not capture the entire spectrum. It is straightforward to adapt this example (e.g. by changing basis) to have the same phenomena without an obvious decomposition of the operator into a finite part and triangular 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 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 T P m , where the parameters n and m are allowed to vary independently. Indeed, the algorithms presented in [24,38] use this method. In effect, we need to know how large n should be to retain enough information of the operator T P m . This type of idea is also used implicitly in the IQR algorithm (see Sect. 4).

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 where v j = 5 sin( j) 2 / √ j if j ≤ 10 and v j = 0 otherwise. As a compact (in fact finite rank) perturbation of the free Laplacian, σ (T 1 ) consists of the interval [−2, 2] together with isolated eigenvalues of finite multiplicity which can be computed [73]. The second operator, T 2 , consists of taking the operator where U k denotes the bilateral shift e j → e j+k , writing this as an operator on l 2 (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 × 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. Example 6.5 (Convergence to extremal parts of the spectrum) To see why we may need some condition on σ (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 Given such a symbol, we define Laurent and Toeplitz operators · · · · · · · · · · · · · · · · · · · · · · · · a 0 a −1 a −2 a −3 a −4 · · · · · · a 1 a 0 a −1 a −2 a −3 · · · · · · a 2 a 1 a 0 a −1 a −2 · · · · · · a 3 a 2 a 1 a 0 a −1 · · · · · · a 4 a 3 a 2 a 1 a 0 · · · · · · · · · · · · · · · · · · · · · · · · a 0 a −1 a −2 · · · a 1 a 0 a −1 · · · a 2 a 1 a 0 · · · · · · · · · · · · · · · ⎞ ⎟ ⎟ ⎠ , acting on l 2 (Z) and l 2 (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 T Q n | P m 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 and Spitzer [62]) If a is a trigonometric polynomial then we have the following convergence in the Hausdorff metric: 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→∞ P m Q * n T (a)Q n | P m H and lim n→∞ P m Q * n L(a)Q n | P m H exist and are either the extreme parts of σ (L(a)) or of ϒ(a). For example, consider the symbols Figure 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 σ (L(a)) (together with 0 if m is not a multiple of 4). Whereas in the case of a it appears that lim n→∞ P m Q * n T (ã)Q n | P m H is the extremal parts of ϒ(ã) and lim n→∞ P m Q * n L(ã)Q n | P m H is the extremal parts of σ (L(ã)) (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 the 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.5i I for 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 σ (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 σ (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 → ∞). 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 ess ) then this extreme point will be recovered in the limit n → ∞ for large enough m.
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 σ (P m T | P m H ), spectral pollution can be detected by changing m (edge states which correspond to spectral pollution are often unstable, but this is not 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 always the case). In general, σ (P m Q * n T Q n | P m 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 The spectrum of T 3 is [− 4, − 2]∪ [2,4]. However, if m is odd then 0 ∈ σ (P m T 3 | P m H ). We shifted the operator by considering T 3 + 0.2I (and then shifted back for the spectrum). Figure 5 shows the Hausdorff distance between σ (P m Q * n (T 3 +0.2I )Q n | P m H )− 0.2I and σ (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.
Another example of this is given by the operator L(a) considered previously. For fixed m we found that Despite this, we found that for small fixed n > 0 it appears that This is shown in Fig. 6 with similar results for L(ã).

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 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 (ã)). 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 Sect. 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 σ 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) Fig. 8 where we have plotted the difference in norm. We also discuss another drawback of the algorithm, m , to compute the pseudospectrum by perturbing the operator T . Let T be the operator obtained from T if we set T 5,4 = 5 × 10 7 (note that this gets rid of the block form). The computation of m involves squaring the operator and hence leads to matrices of norm of order 10 15 , making it impossible to compute σ ( T ) using double precision for 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 × 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. Example 6.9 (P T -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 operatorx → −x,p → −p and T the time operatorp → −p, i → − 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 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 [45,60].
We consider an operator on l 2 (Z) of the form 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 = cos(n) + iγ sin(n), mod (n, 2) = 0 0, mod (n, 2) = 1 , (6.5) and found similar results for other potentials. Figure 9 shows the same qualitative behaviour as the last example for γ = 1, 2 at m = 500, n = 3000. We shifted by 2.2 and 2.15 for γ = 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 → ∞) caused σ (P m Q * n H 1 Q n | P m H ) to approach the extremal part of the spectrum.

Numerical examples III: random non-Hermitian operators and boundary conditions
In this final section, we explore examples where the P m Q * n T Q n P m naturally give rise to periodic boundary conditions (this was already seen for some examples of Laurent operators in Sect. 6.2). Both examples discussed here are physically motivated random tridiagonal operators on the lattice 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 T P 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 superdiagonals, first studied by Feinberg and Zee [23,31,40]. The usual "Hopping Sign Model" is defined via with b n ∈ {± 1} (say independent Bernoulli with parameter p = 1/2). This describes a particle "'hopping" on Z and can be mapped into a (complex-valued) random walk. We will consider a slightly different operator described by 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 are i.d.d. random variables with Bernoulli distribution p. In other words 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 ± 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 σ (H 3 ) is contained in the annulus {z ∈ C : 2 sinh(g) ≤ |z| ≤ 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 (Z) then one obtains matrices of the form 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 . Note that after a few iterates, the output 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 On the other hand, the use of periodic boundary conditions leads to the matrix which does not suffer from this setback.
In [1] this phenomenon was studied via localisation of the eigenvalues of M 2 n , in particular using the Lyapunov exponent κ(z) which is equal to the inverse of the localisation length. An eigenfunction ψ with eigenvalue z localised around x 0 behaves approximately as If one defines recursively y n+1 (z) = exp(g) ψ n+2 ψ n+1 = −(s − n−1 /s + n )/y n (z) + z/s + n then (in the limit of large system sizes) This is known as the transfer matrix approach. For fixed z, as we increase g, κ(z; g) becomes negative. The heuristic is that a hole opens up in the spectrum corresponding to a negative Lyapunov exponent. Eigenvalues of M 2 n 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 1 n inside the negative κ zone correspond to edge states due to the finite system size approximation. Figure 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 κ = 0. In order to calculate κ, 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 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 produce spectral pollution in the hole with localised eigenfunctions and the contour κ = 0 corresponds to the delocalised region. In order to compare to the spectrum of the infinite operator on l 2 (Z) we have plotted σ (H 3 ), for = 10 −2 , calculated using matrix sizes of order 10 5 . We note that the spectrum is independent of p ∈ (0, 1) so we have also shown the union of these estimates over p = {k/100} 99 k=1 . Although the algorithm used to compute the pseudospectrum is guaranteed to converge to σ (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. Upon taking smaller, we found 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 is contained in σ (H 3 ) and that the spectrum (and pseudospectrum) of H 3 has fourfold rotational symmetry. Denoting the rotation of E 1 by π/4 as E 2 we have shown E 1 ∪ E 2 in the figure. 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 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 nonnormal 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 [39]. Their model showed that an imaginary gauge field in a disordered one-dimensional lattice can induce a delocalisation transition. The operator in B(l 2 (Z)) can be written as where g > 0 and V is a random potential. This operator also has applications in population biology [52] 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(±g). Standard finite section produces real eigenvalues since the matrix P m H 4 | P m 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θ) + exp(−g − iθ) : θ ∈ [0, 2π)}, 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 [31]. Goldsheid and Khoruzhenko have studied the convergence of the spectral measure in the periodic case as N → ∞ in [33], N being the number of sites. In general, the support of these measures as N → ∞ can be very different from the spectrum of the operator on l 2 (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 {±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 bounds the spectrum:  Fig. 11 The output of IQR over 200 samples computing σ (P m Q * n H 4 Q n | P m 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 where 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 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 H . Shifting by +4i I caused the output to recover the top part of the limit curve.

Remark 6.12
For any operator T that has Q n unitary, the essential spectrum and spectrum of Q * n T Q 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.
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 [57]. 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.11), and below the essential spectral 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 Sect. 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. 4 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 −α ) as n → ∞ 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 → ∞, the output σ (P m Q * n T Q 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 for many operators even when taking a mixture of double limits m, n → ∞. However, examples of Laurent operators were given where it appears σ (P m Q * n T Q n | P m ) converges to the spectrum as m → ∞ 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 σ (P m Q * n(m) T Q n(m) | P m ) converges to the spectrum as m → ∞? 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?

A.2 Recalling the basics of the SCI hierarchy
The cornerstone 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: The set is the set of objects that give rise to our computational problems. The problem function : → M is what we are interested in computing. Moreover, the set is the collection of functions that provide us with the information we are allowed to read. This leads to the following definition.  The goal is to find algorithms which approximate the function . 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. 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 such that for T ∈ we have that (T ) ⊂ M. In this case, we will still require that a general algorithm produces a single valued out put i.e (T ) ∈ M for T ∈ . However, we replace the metric in order to define convergence. In particular, n (T ) → (T ), as n → ∞ means 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. In other words one may say that for the finitely many steps of the computation of the lowest functions = n k ,...,n 1 : → M only the four arithmetic operations +, −, ·, / within the smallest (algebraic) field which is generated by the input {A f } f ∈ (A) are allowed. In addition, we allow the extraction of 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 { , , M, } is said to have Solvability Complexity Index SCI( , , M, ) α = k, with respect to a tower of algorithms of type α, if k is the smallest integer for which there exists a tower of algorithms of type α of height k. If no such tower exists then SCI( , , M, ) α = ∞. If there exists a tower { n } n∈N of type α and height one such that = n 1 for some n 1 < ∞, then we define SCI( , , M, ) α = 0. We may sometimes write SCI( , ) α to simplify notation when M and are obvious.
The definition of the SCI immediately induces the SCI hierarchy: Definition A.10 (The Solvability Complexity Index Hierarchy) Consider a collection C of computational problems and let T be the collection of all towers of algorithms of type α for the computational problems in C. Define Remark A.11 (The k notation) Note that in this paper we only consider radical towers and hence the superscript α will be omitted throughout. Thus we will always write k .
Finally, we recall the definition of α 1 .