The foundations of spectral computations via the Solvability Complexity Index hierarchy: Part I
TTHE FOUNDATIONS OF SPECTRAL COMPUTATIONSVIA THE SOLVABILITY COMPLEXITY INDEX HIERARCHY: PART I
M. J. COLBROOK AND A. C. HANSENA
BSTRACT . The problem of computing spectra of operators is arguably one of the most investigated areas ofcomputational mathematics. However, the question of computing spectra of general infinite matrices has untilrecently been open. This recent progress and the current paper reveal that, unlike the finite-dimensional case,infinite-dimensional problems yield a highly intricate infinite classification theory determining which spectralproblems can be solved and with which type of algorithms. Classifying spectral problems and providing optimalalgorithms is uncharted territory in the foundations of computational mathematics, and this paper is the first of atwo-part series establishing the foundations of computational spectral theory through the Solvability ComplexityIndex (SCI) hierarchy and has three purposes. First, we establish answers to many longstanding open questions onthe existence of algorithms. We show that for large classes of partial differential operators on unbounded domains,spectra can be computed with error control from point sampling operator coefficients. Further results includecomputing spectra of (possibly unbounded) operators on graphs with error control, the computational spectralgap problem, computing spectral classifications, and computing discrete spectra, multiplicities and eigenspaces.Second, these classifications determine which types of problems can be used in computer-assisted proofs. Thetheory for this is virtually non-existent, and we provide some of the first results in this infinite classification theory.Third, our proofs are constructive, yielding a library of new algorithms and techniques that handle problems thatbefore were out of reach. We show several examples on contemporary problems in the physical sciences. Ourapproach is closely related to Smale’s program on the foundations of computational mathematics initiated inthe 1980s, as many spectral problems can only be computed via several limits, a phenomenon shared with thefoundations of polynomial root finding with rational maps, as proved by McMullen. C ONTENTS
1. Introduction 22. Classifications in the SCI Hierarchy 33. Main Results 64. Connection to Previous Work 145. Mathematical Preliminaries 156. Proofs of Theorems on Unbounded Operators on Graphs 187. Proofs of Theorems on Differential Operators on Unbounded Domains 268. Proofs of Theorems on Discrete Spectra 359. Proof of Theorem on the Spectral Gap and Spectral Classification 4210. Computational Examples 44References 50Appendix A. Computational Routines 54 D EPARTMENT OF A PPLIED M ATHEMATICS AND T HEORETICAL P HYSICS , U
NIVERSITY OF C AMBRIDGE
E-mail addresses : [email protected], [email protected] . Key words and phrases.
Computational spectral problem, Solvability Complexity Index hierarchy, Smale’s program on the founda-tions of computational mathematics, computer-assisted proofs2020
Mathematics Subject Classification. a r X i v : . [ m a t h . SP ] A ug FOUNDATIONS OF SPECTRAL COMPUTATIONS
1. I
NTRODUCTION
The problem of computing spectra of operators has fascinated and frustrated mathematicians for severaldecades since the 1950s resulting in a vast literature (see § Unfortunately, there is a dearth of literature on this basic problem, and so far as we have been ableto tell, there are no proven techniques ” [5]. This longstanding problem for general infinite matrices, and thequestion why there have been no known general techniques, have recently been addressed [8,73] and are dueto classification results in the newly established Solvability Complexity Index (SCI) hierarchy [8, 43, 73].The fact that algorithms were not found for the general computational spectral problem has a potentiallysurprising cause: one needs several limits in the computation. Traditional approaches have been dominatedby techniques based on one limit, and this is the reason behind Arveson’s observation. Moreover, the factthat several limits are required is a phenomenon that is shared by other areas of computational mathematics.For example, the problem of root-finding of polynomials with rational maps initiated by S. Smale [102] isalso subject to the issue of requiring several limits. This was established by C. McMullen [83, 84] and P.Doyle & C. McMullen in [53], and their results become classification results in the SCI hierarchy.The recent results in [8,40,42,43,73] establishing the SCI hierarchy reveal that the computational spectralproblem becomes an infinite classification theory. Hence, there is an infinite well of open problems, some ofwhich have been open for decades. For example, the following issue, even when neglecting the requirementof an error parameter, has been open since the early days of spectral computations in the 1950s:
For which classes of differential operators on unbounded domains do there exist algorithmsthat converge to the true spectrum, and also guarantee that the output is in the spectrum upto an arbitrary small (cid:15) > parameter (the problem is in Σ in the SCI hierarchy language)?In other words, the algorithm will never make a mistake. There is, of course, a vast literature on computing spectra of differential operators on bounded domains,but these techniques will typically yield non-convergent methods in the unbounded domain case, and evenfor bounded domains obtaining error bounds is, in general, well known to be very difficult. The main topicof this paper is to provide solutions to many of these problems, and this program has two main motivations:(I)
Classifications and new algorithms : Sharp classifications of problems in the SCI hierarchy establishthe boundaries of what computers can achieve. Constructive classifications, which we always providein this paper, provide algorithms that realise these boundaries. Moreover, with new classifications, suchalgorithms will solve problems in the sciences that before were not possible. We provide several examplesin this paper.(II)
Computer-assisted proofs : Computer-assisted proofs, where computers are used to solve numericalproblems rigorously, have become essential in modern mathematics. What may be surprising is that unde-cidable or non-computable problems can be used in computer-assisted proofs. Indeed, the recent proof ofKepler’s conjecture (Hilbert’s 18th problem) [69, 70], led by T. Hales, on optimal packings of -spheres,relies on such undecidable problems. Moreover, the Dirac–Schwinger conjecture on the asymptotic be-haviour of ground states of certain Schr¨odinger operators was proven in a series of papers by C. Feffermanand L. Seco [55–63] using computer assistance. Fascinatingly, this proof also relies on computing non-computable problems. This may seem like a paradox, but can indeed be explained by the SCI hierarchy.In particular, it is the Σ A class described below that is crucial. In fact, Hales, Fefferman and Seco implic-itly prove Σ A classifications in the SCI hierarchy in their papers. Our classifications of spectral problemsprovide new results regarding which spectral problems can be used in computer-assisted proofs.A summary of the main results of this paper is given in Table 1. The main theorems are contained in § OUNDATIONS OF SPECTRAL COMPUTATIONS 3 (i)
Computing spectra of differential operators . Linked to computational PDE theory, there is a rich litera-ture on computing spectra of differential operators on bounded domains (see [17,18,34–36,39,90,91,114]for a small sample). However, it is, in general, unknown how to compute spectra of differential operatorson unbounded domains. We provide a sharp solution to this problem for large classes of differential opera-tors, meaning that we achieve the boundary of what computers can achieve on these problems. We provideconvergent algorithms that are guaranteed to produce output that is in the spectrum up to an arbitrarilysmall error chosen by the user. As such, these algorithms can be used in computer-assisted proofs.(ii)
Computing spectra of unbounded operators on graphs . In the discrete setting, operators on l ( N ) or,more generally, graphs or lattices are ubiquitous in mathematics and physics. We establish sharp classi-fications of spectral problems for such operators, and in many cases, we establish convergent algorithmswith guaranteed error control on the output. Hence these algorithms may also be used in computer-assistedproofs. We also consider the decision problem of determining whether the spectrum (or pseudospectrum)intersects a given compact set.(iii) The spectral gap problem . The spectral gap problem has a long tradition and is linked to many impor-tant conjectures and problems such as the Haldane conjecture [68] or the Yang–Mills mass gap problemin quantum field theory [20]. The problem consists of determining whether there is a gap between thelowest element in the spectrum and the next element. We show why this problem is notoriously difficult,as this problem is higher up in the SCI hierarchy, even for the simplest of operators. This means that noalgorithm can provide verifiable results on a digital computer, and hence these problems cannot be usedin computer-assisted proofs. This is extended to the problem of spectral classification at the bottom of thespectrum.(iv)
Computing discrete spectra and multiplicities . We demonstrate why this is also a difficult problem byestablishing the correct classification high up in the SCI hierarchy. However, the sharp algorithm weprovide is still practical, since its first limit is always contained in the discrete spectrum and one can obtaina distance function of each point of the output to the spectrum. We extend these results to computingmultiplicities, eigenspaces and determining if the discrete spectrum is non-empty.The rest of this paper is organised as follows. In § §
5. The main results aregiven in § §
4. Proofs are given in § §
9. Finally, somecomputational examples are given in §
10, and pseudocode is provided in Appendix A.2. C
LASSIFICATIONS IN THE
SCI H
IERARCHY
The SCI hierarchy.
The fundamental notion of the SCI hierarchy is that of a computational problem.The SCI of a class of computational problems is the smallest number of limits needed in order to computethe solution to the problem. The basic objects of a computational problem are: Ω , called the domain , Λ a setof complex-valued functions on Ω , called the evaluation set , ( M , d ) a metric space, and Ξ : Ω → M the problem function . The set Ω is the set of objects that give rise to our computational problems. The problemfunction Ξ : Ω → M is what we are interested in computing. Finally, the set Λ is the collection of functionsthat provide us with the information we are allowed to read as input to the algorithm. This leads to thefollowing definition: Definition 2.1 (Computational problem) . Given a domain Ω ; an evaluation set Λ , such that for any A , A ∈ Ω , A = A if and only if f ( A ) = f ( A ) for all f ∈ Λ ; a metric space M ; and a problem function Ξ : Ω → M , we call the collection { Ξ , Ω , M , Λ } a computational problem. The definition of a computational problem is deliberately general in order to capture any computationalproblem in the literature. However, the set-up of this paper has the following typical form: Ω will be a class FOUNDATIONS OF SPECTRAL COMPUTATIONS
Description of Problem SCI Hierarchy Classification Theorems
Computing spectrum/pseudospectrum of differentialoperators whose coefficients have bounded total varia-tion from point evaluations of coefficients. ∈ Σ A , / ∈ ∆ G (see Theorem for relaxations) 3.3Computing spectrum/pseudospectrum of differentialoperators whose coefficients are analytic from powerseries of coefficients. ∈ Σ A , / ∈ ∆ G (see Theorem for relaxations) 3.5Computing spectrum/pseudospectrum of unboundedoperators with known bounded dispersion and knownresolvent bound. ∈ Σ A , / ∈ ∆ G (same for diagonal operators) 3.8Determining if the spectrum/pseudospectrum of anoperator with known bounded dispersion intersects acompact set. ∈ Π A , / ∈ ∆ G (same for diagonal operators) 3.9Spectral gap problem. ∈ Σ A , / ∈ ∆ G (same for diagonal operators) 3.11Spectral classification problem. ∈ Π A , / ∈ ∆ G (same for diagonal operators) 3.11Computing Sp d ( A ) (and multiplicities of eigenvalues)for bounded normal operators. With bounded dispersion: ∈ Σ A , / ∈ ∆ G (same for diagonal operators)Multiplicities have generalised ∈ Π A Without bounded dispersion: ∈ Σ A , / ∈ ∆ G ∈ Σ A , / ∈ ∆ G Without bounded dispersion: ∈ Σ A , / ∈ ∆ G ABLE
1. Summary of the main results, pseudocode for the algorithms is provided inAppendix A. Bounded dispersion roughly means that we know the asymptotic off-diagonaldecay on the matrix elements in the matrix representation of the operator, see (3.9). Also,known resolvent bound means control of the growth of the resolvent near the spectrum,see (3.4) and (3.10).of operators on a separable Hilbert space H , Ξ( A ) = Sp( A ) (the spectrum or other related maps), ( M , d ) is the collection of closed subsets of C with an appropriate generalisation of the Hausdorff metric (see (3.1)and (3.2)), and Λ may be the set of complex functions that could provide the matrix elements of A ∈ Ω givensome orthonormal basis { e j } of H . In particular, Λ consists of f i,j : A (cid:55)→ (cid:104) Ae j , e i (cid:105) , i, j ∈ N , which providethe entries of the matrix representation of A with respect to the basis. Moreover, Λ could be the collection offunctions providing point samples of a potential (or coefficient) function of a Schr¨odinger (or more general)partial differential operator.In words, the SCI hierarchy [8, 73] for spectral problems can be informally described as follows, and fordecision problems, the description is similar (see § The SCI hierarchy:
Given a collection C of computational problems,(i) ∆ α = Π α = Σ α is the set of problems that can be computed in finite time, the SCI = 0 .(ii) ∆ α is the set of problems that can be computed using one limit (the SCI = 1 ) with control of theerror, i.e. ∃ a sequence of algorithms { Γ n } such that d (Γ n ( A ) , Ξ( A )) ≤ − n , ∀ A ∈ Ω .(iii) Σ α : We have ∆ α ⊂ Σ α ⊂ ∆ α and Σ α is the set of problems for which there exists a sequence ofalgorithms { Γ n } such that for every A ∈ Ω we have Γ n ( A ) → Ξ( A ) as n → ∞ . However, Γ n ( A ) is always contained in a set X n ( A ) such that d ( X n , Ξ( A )) ≤ − n . OUNDATIONS OF SPECTRAL COMPUTATIONS 5 F IGURE
1. Meaning of Σ and Π convergence for problem function Ξ computed in theHausdorff metric. The red areas represent Ξ( A ) , whereas the green areas represent theoutput of the algorithm Γ n ( A ) . Σ convergence means convergence as n → ∞ but eachoutput point in Γ n ( A ) is at most distance − n from Ξ( A ) . Similarly, in the case of Π , wehave convergence as n → ∞ but any point in Ξ( A ) is at most distance − n from Γ n ( A ) .(iv) Π α : We have ∆ α ⊂ Π α ⊂ ∆ α and Π α is the set of problems for which there exists a sequence ofalgorithms { Γ n } such that for every A ∈ Ω we have Γ n ( A ) → Ξ( A ) as n → ∞ . However, thereexists sets X n ( A ) such that Ξ( A ) ⊂ X n ( A ) and d ( X n , Γ n ( A )) ≤ − n .(v) ∆ α is the set of problems that can be computed using one limit (the SCI = 1 ) without error control,i.e. ∃ a sequence of algorithms { Γ n } such that lim n →∞ Γ n ( A ) = Ξ( A ) , ∀ A ∈ Ω .(vi) ∆ αm +1 , for m ∈ N , is the set of problems that can be computed by using m limits, (the SCI ≤ m ),i.e. ∃ a family of algorithms { Γ n m ,...,n } such that lim n m →∞ . . . lim n →∞ Γ n m ,...,n ( A ) = Ξ( A ) , ∀ A ∈ Ω . (vii) Σ αm is the set of problems that can be computed by passing to m limits, and computing the m -thlimit is a Σ α problem.(viii) Π αm is the set of problems that can be computed by passing to m limits, and computing the m -thlimit is a Π α problem.Schematically, the SCI hierarchy can be viewed in the following way.(2.1) Π α Π α Π α ∆ α ∆ α Σ α ∪ Π α ∆ α Σ α ∪ Π α ∆ α · · · Σ α Σ α Σ α == (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) (cid:40) Note that the Σ α and Π α classes become crucial in computer-assisted proofs, as they guarantee algorithmsthat will not make mistakes, see § Remark 2.2 (The model of computation α ) . The α in the superscript indicates the model of computation,which is described in §
5. For α = G , the underlying algorithm is general and can use any tools at its disposal.The reader may think of a Blum–Shub–Smale (BSS) machine [14] or a Turing machine [108] with accessto any oracle, although a general algorithm is even more powerful. However, for α = A this means thatonly arithmetic operations and comparisons are allowed. In particular, if rational inputs are considered, thealgorithm is a Turing machine, and in the case of real inputs, a BSS machine. Hence, a result of the form / ∈ ∆ Gk is stronger than / ∈ ∆ Ak . FOUNDATIONS OF SPECTRAL COMPUTATIONS
Indeed, a / ∈ ∆ Gk result is universal and holds for any model of computation. Moreover, ∈ ∆ Ak is stronger than ∈ ∆ Gk , and similarly for the Π k and Σ k classes. The main results are sharp classification results in this hierarchythat are summarised in Table 1.2.2. The SCI hierarchy and computer-assisted proofs.
Note that ∆ A is the class of problems that arecomputable according to Turing’s definition of computability [108]. In particular, there exists an algorithmsuch that for any (cid:15) > , the algorithm can produce an (cid:15) -accurate output. Most infinite-dimensional spectralproblems, unlike the finite-dimensional case, are / ∈ ∆ A . The simplest way to see this is to consider theproblem of computing spectra of infinite diagonal matrices. Since this problem is the simplest of the infinitecomputational spectral problems and does not lie in ∆ A , very few interesting infinite-dimensional spectralproblems are actually in ∆ A . This is why most of the literature on spectral computations provides algorithmsthat yield ∆ A classification results. In particular, an algorithm will converge, but error control may not bepossible.Problems that are not in ∆ A are computed daily in the sciences, simply because numerical simulationsmay be suggestive rather than providing a rock-solid truth. Moreover, the lack of error control may becompensated for by comparing with experiments. However, this is not possible in computer-assisted proofs,where rigour is the only approach accepted. It may, therefore, be surprising that there are examplesof famous conjectures that have been proven with numerical calculations of problems that are not in ∆ A ,i.e. problems that are non-computable according to Turing. A striking example is the proof of Kepler’sconjecture [69, 70], where the decision problems computed are not in ∆ A . The decision problems are ofthe form of deciding feasibility of linear programs given irrational inputs, shown in [7] to not lie in ∆ A .Similarly, the problem of obtaining the asymptotics of the ground state of the operator H dZ = d (cid:88) k =1 ( − ∆ x k − Z | x k | − ) + (cid:88) ≤ j ≤ k ≤ d | x j − x k | − , as Z → ∞ was obtained by a computer-assisted proof [55–63] by Fefferman and Seco, proving the Dirac-Schwinger conjecture, that relied on problems that were not in ∆ A . The SCI hierarchy can describe theseparadoxical phenomena.2.2.1. The Σ A and Π A classes. The key to the paradoxical phenomena lies in the Σ A and Π A classes.These classes of problems are larger than ∆ A , but can still be used in computer-assisted proofs. Indeed,if we consider computational spectral problems that are in Σ A , then there is an algorithm that will neverprovide incorrect output. The output may not include the whole spectrum, but it is always sound. Thus,conjectures about operators never having spectra in a certain area could be disproved by a computer-assistedproof. Similarly, Π A problems would always be approximated from above, and thus conjectures on thespectrum being in a certain area could be disproved by computer simulations.In both of the above examples (the proof of the Dirac-Schwinger conjecture and Kepler’s conjecture), oneimplicitly shows that the relevant computational problems in the computer-assisted proofs are in Σ A .3. M AIN R ESULTS
The main results are sharp classifications in the SCI hierarchy with corresponding algorithms that settlesome of the many open classification problems in computational spectral theory. We are concerned with thefollowing problem:
Given a computational spectral problem, where is it in the SCI hierarchy?
In addition to the spectrum, we also consider the pseudospectrum defined by Sp (cid:15) ( A ) := cl( { z ∈ C : (cid:107) ( A − zI ) − (cid:107) > /(cid:15) } ) , (cid:15) > . OUNDATIONS OF SPECTRAL COMPUTATIONS 7
In this paper, we consider the following four main problems: computing spectra of general differentialoperators, computing spectra of unbounded operators on graphs, the computational spectral gap problem,and computing discrete spectra with multiplicities. A formal definition of the SCI hierarchy is given in § ( M , d ) be the set of all non-empty compactsubsets of C provided with the Hausdorff metric d = d H :(3.1) d H ( X, Y ) = max (cid:26) sup x ∈ X inf y ∈ Y d ( x, y ) , sup y ∈ Y inf x ∈ X d ( x, y ) (cid:27) , where d ( x, y ) = | x − y | is the usual Euclidean distance. In the case of unbounded operators, we use the Attouch–Wets metric defined by(3.2) d AW ( C , C ) = ∞ (cid:88) n =1 − n min (cid:40) , sup | x |≤ n | dist( x, C ) − dist( x, C ) | (cid:41) , for C , C ∈ Cl( C ) , where Cl( C ) denotes the set of closed non-empty subsets of C .3.1. Computing spectra of differential operators on unbounded domains.
There is a rich literature onhow to compute spectra of differential operators on bounded domains that is intimately linked to compu-tational PDE theory. The computation is often done with finite element, finite difference or spectral meth-ods by discretising the operator on a suitable finite-dimensional space and then using algorithms for finite-dimensional matrix eigenvalue problems on the discretised operator [17,18,34–36,39,90,91,114]. However,it is in general unknown how to compute spectra of differential operators on unbounded domains, or wherethis problem is in the SCI hierarchy.For N ∈ N , consider the operator formally defined on L ( R d ) by(3.3) T u ( x ) = (cid:88) k ∈ Z d ≥ , | k |≤ N a k ( x ) ∂ k u ( x ) , where throughout we use multi-index notation with | k | = max {| k | , ..., | k d |} and ∂ k = ∂ k x ∂ k x ...∂ k d x d . Wewill assume that the coefficients a k ( x ) are complex-valued measurable functions on R d . Suppose also that T can be defined on an appropriate domain D ( T ) such that T is closed and has a non-empty spectrum. Ouraim is to compute the spectrum and (cid:15) -pseudospectrum from the functions a k . We consider two cases. First,the algorithm can access point samples of the functions, and second, the algorithm can access coefficients inthe series expansion of the functions in the case that the a k are analytic. Note that these are very differentcomputational problems. Remark 3.1 (The open problem of computing spectra of differential operators) . There is no existing theoryguaranteeing even a finite SCI for this problem even when each a k is a polynomial. For simple polynomials,arising for example as potentials of Schr¨odinger operators, the standard procedure is to discretise the differ-ential operator via finite differences, truncate the resulting infinite matrix and then handle the finite matrixwith standard algorithms designed for finite-dimensional problems. Such an approach would at best give a ∆ A classification, and, in general, this approach may not always converge. Despite this, we prove below thatone can actually achieve Σ A classification for a large class of operators.3.1.1. The set-up.
To make our problems well-defined, we let Ω consist of all such T such that the followingassumptions hold:(1) The set C ∞ ( R d ) of smooth, compactly supported functions forms a core of T and its adjoint T ∗ .(2) The adjoint operator T ∗ can be initially defined on C ∞ ( R d ) via T ∗ u ( x ) = (cid:88) k ∈ Z d ≥ , | k |≤ N ˜ a k ( x ) ∂ k u ( x ) , where ˜ a k ( x ) are complex-valued measurable functions on R d . FOUNDATIONS OF SPECTRAL COMPUTATIONS (3) For each of the functions a k ( x ) and ˜ a k ( x ) , there exists a positive constant A k and an integer B k such that | a k ( x ) | , | ˜ a k ( x ) | ≤ A k (cid:16) | x | B k (cid:17) , almost everywhere on R d , that is, we have at most polynomial growth.(4) We have access to a sequence { g m } of strictly increasing continuous functions g m : R ≥ → R ≥ vanishing only at and diverging at ∞ such that(3.4) g m (dist( z, Sp( T ))) ≤ (cid:107) R ( z, T ) (cid:107) − , ∀ z ∈ B m (0) . In this case we say that T has resolvent bounded by { g m } . Note that this implicitly assumes that Sp( T ) (and hence Sp (cid:15) ( T ) ) is non-empty.Hence, we consider the operator T defined as the closure of T acting on C ∞ ( R d ) . The initial domain C ∞ ( R d ) is commonly encountered in applications, and it is straightforward to adapt our methods to otherinitial domains such as Schwartz space. Remark 3.2.
In order to handle non-self-adjoint operators, we need to be able to control the resolvent asin (3.4). Without such control, the spectral problem is still not in ∆ G . If T has Sp( T ) (cid:54) = ∅ , then a simplecompactness argument implies the existence of such a sequence of continuous functions. We may not beable to control the growth of the resolvent across the whole complex plane by a single function, but forconvergence of our algorithms, it is enough to control the resolvent on compact balls. For self-adjoint (and,more generally, normal) T , we can take g m ( x ) = x . In what follows, the functions { g m } are not needed tocompute pseudospectra.3.1.2. General case with function evaluations.
In this section we consider the computation of spectra andpseudospectra of operators T ∈ Ω from evaluations of the functions a k and ˜ a k . For dimension d and r > consider the space(3.5) A r = { f ∈ M ([ − r, r ] d ) : (cid:107) f (cid:107) ∞ + TV [ − r,r ] d ( f ) < ∞} , where M ([ − r, r ] d ) denotes the set of measurable functions on the hypercube [ − r, r ] d and TV [ − r,r ] d the totalvariation norm in the sense of Hardy and Krause [85]. This space becomes a Banach algebra when equippedwith the norm (cid:107) f (cid:107) A r = (cid:107) f (cid:107) ∞ + σ TV [ − r,r ] d ( f ) with σ = 3 d + 1 [16]. We will assume that each of the(appropriate restrictions of) a k and ˜ a k lie in A r for all r > and that we are given a sequence of positivenumbers such that(3.6) (cid:107) a k (cid:107) A n , (cid:107) ˜ a k (cid:107) A n ≤ c n , c n > , n ∈ N , | k | ≤ N. This information is completely analogous to using bounded dispersion for matrix problems encountered in § Σ . Let Ω = { T ∈ Ω | such that assumptions (1) – (4) and (3.6) hold } . In this case, we let Λ contain functions that allow us to access sample of the functions { g m } m ∈ N and { a k , ˜ a k } | k |≤ N , as well as the constants { A k , B k } | k |≤ N , { c n } n ∈ N . Consider the weaker assumption on Λ that we can evaluate b n > (and not the A k , B k ’s and the c n ’s) such that(3.7) sup n ∈ N max {(cid:107) a k (cid:107) A n , (cid:107) ˜ a k (cid:107) A n : | k | ≤ N } b n < ∞ . With a slight abuse of notation, we use Ω to denote the class of problems where we have this weakerrequirement. We can now define the mappings Ξ j , Ξ j : Ω , Ω (cid:51) T (cid:55)→ Sp( T ) ∈ M AW , j = 1Sp (cid:15) ( T ) ∈ M AW , j = 2 , and state the first theorem. OUNDATIONS OF SPECTRAL COMPUTATIONS 9
Theorem 3.3 (Differential operators and point samples) . Let Ξ j , Ξ j , Ω and Ω be as above. Then for j = 1 , , we have that ∆ G (cid:54)(cid:51) { Ξ j , Ω } ∈ Σ A and Σ G ∪ Π G (cid:54)(cid:51) { Ξ j , Ω } ∈ ∆ A . Remark 3.4.
The proof also shows that even if we had included the information { A k , B k } | k |≤ N for oper-ators in Ω , we would still have { Ξ j , Ω } / ∈ Σ G ∪ Π G . Moreover, though we have chosen R d as thegeometrical domain of our operators, the proof will make clear that other domains (and correspondinglyother discretisations) can also be treated, yielding a Σ A classification. Indeed, the results can be extended toany domain where we can build a suitable basis to represent the operator, such as the half-line for radiallysymmetric Dirac operators in quantum chemistry, intervals using orthogonal polynomial series, or productsof any of the above geometries. It is also possible to extend our results to more complicated domains usingfinite elements, non-orthogonal bases and generalised pencil eigenvalue problems, but this will be the topicof future work.3.1.3. Analytic coefficients.
In this section, we assume that the functions a k and ˜ a k are analytic. In partic-ular, we assume we can evaluate { c j } j ∈ N , an enumeration (where we know the ordering) of the coefficients a mk where a k ( x ) = (cid:80) m ∈ ( Z ≥ ) d a mk x j . In this special case, we can compute the corresponding coefficients ofthe ˜ a k ( x ) using finitely many arithmetic operations on { c j } . We will assume that as well as the information { g m } , { c j } and { A k , B k } , our algorithms can read the following information. Given a k ( x ) = (cid:88) m ∈ ( Z ≥ ) d a mk x m , ˜ a k ( x ) = (cid:88) m ∈ ( Z ≥ ) d ˜ a mk x m , for each n ∈ N we know a constant d n such that(3.8) | a mk | , | ˜ a mk | ≤ d n ( n + 1) −| m | , ∀ m ∈ ( Z ≥ ) d , | k | ≤ N. It is straightforward to show that such a d n must exist using the fact that the power series converges absolutelyon the whole of R d . Let Ω = { T ∈ Ω | such that (1) – (4), the functions a k are analytic and (3.8) hold } . Moreover, in this case we let Λ contain functions that allow us to access sample of the functions { g m } m ∈ N ,the constants { A k , B k } | k |≤ N , { c n } n ∈ N , and { d n } n ∈ N . As the proof makes clear, the information d n canbe replaced by any suitable information that allows us to control the remainder term in the truncated Taylorseries uniformly on compact subsets of R d . For example, we could use Cauchy’s formula, together withbounds on the functions a k on compact subsets of C d . One could consider a weaker requirement on Λ byreplacing knowledge of A k , B k and d n by some sequence of positive numbers b n with sup n ∈ N sup m ∈ ( Z ≥ ) d max {| a mk | ( n + 1) | m | , | ˜ a mk | ( n + 1) | m | : | k | ≤ N } b n < ∞ . With a slight abuse of notation, we use Ω to denote the class of problems where we have this weakerrequirement. Moreover, let Ω p denote the class of operators in Ω such that each a k is a polynomial(where we can let b n be n ! say). We can now define the mappings Ξ j , Ξ j : Ω , Ω , Ω p (cid:51) T (cid:55)→ Sp( T ) ∈ M AW , j = 1Sp (cid:15) ( T ) ∈ M AW , j = 2 , and state the second theorem. Theorem 3.5 (Differential operators and analytic coefficients) . Let Ξ j , Ξ j , Ω , Ω and Ω p be as above.Then for j = 1 , , we have that ∆ G (cid:54)(cid:51) { Ξ j , Ω } ∈ Σ A , Σ G ∪ Π G (cid:54)(cid:51) { Ξ j , Ω } ∈ ∆ A , and Σ G ∪ Π G (cid:54)(cid:51) { Ξ j , Ω p } ∈ ∆ A . The new algorithms that yield the Σ A results are useful not only on unbounded domains but also onbounded domains. Indeed, whilst standard algorithms for computing spectra of differential operators onbounded domains often have results on qualitative rates of convergence, typically they do not have the abovefeature of error control and it can be difficult to determine which portion of the computation can be trusted.This is a well-known problem, see for example [113], which occurs even if the algorithm is convergent. Forexample, this means that such algorithms cannot be used for computer-assisted proofs. In the language ofthe SCI hierarchy, these standard algorithms provide, at best, ∆ A classifications of the problems, and notthe correct Σ A classification. In other words, the standard algorithms do not realise the boundary of whatalgorithms can actually do. Hence, we can draw the following conclusion: Computing spectra of differential operators through discretisation of the operator that yield a finite-dimensional matrix, for which one computes its eigenvalues, is typically not an optimal method. Suchmethods will typically not yield the sharp Σ A classification providing certainty about the output. How-ever, as demonstrated above, there do exist algorithms that are optimal in the way that they provideerror control and certainty about the computed output. Computing spectra of unbounded operators on graphs.
Consider a possibly unbounded operator A with domain D ( A ) ⊂ l ( N ) and non-empty spectrum. We consider the problems of computing Ξ ( A ) = Sp( A ) , Ξ ( A ) = Sp (cid:15) ( A ) . To define the computational problem, we have to define the domain Ω as well as Λ , the set of evaluationfunctions. Let C ( l ( N )) denote the set of closed, densely defined operators on l ( N ) , and consider thefollowing assumptions.(1) The subspace span { e n : n ∈ N } forms a core for A and A ∗ , where { e j } j ∈ N is the canonical basisfor l ( N ) .(2) Given any f : N → N with f ( n ) ≥ n define D f,n ( A ) := max (cid:8)(cid:13)(cid:13) ( I − P f ( n ) ) AP n (cid:13)(cid:13) , (cid:13)(cid:13) ( I − P f ( n ) ) A ∗ P n (cid:13)(cid:13)(cid:9) , (3.9) where P n is the projection onto the span of { e , . . . , e n } of the canonical basis. We say that an oper-ator has bounded dispersion with respect to f if lim n →∞ D f,n ( A ) = 0 . We will assume knowledgeof a null sequence { c n } n ∈ N ⊂ Q with D f,n ( A ) ≤ c n .(3) As in the case of § { g m } (see (3.4) and the assumptions on { g m } )such that(3.10) g m (dist( z, Sp( A ))) ≤ (cid:107) R ( z, A ) (cid:107) − , ∀ z ∈ B m (0) . In this case we say that A has resolvent bounded by { g m } . Note that this implicitly assumes that thespectrum is non-empty. Remark 3.6.
The concept of bounded dispersion in (3.9) generalises the notion of a banded matrix. More-over, given any operator with assumption (1), there exists an f such that lim n →∞ D f,n ( A ) = 0 . The theoremwe prove is for the class of operators that have lim n →∞ D f,n ( A ) = 0 given a fixed f . The function f willbe used to construct certain rectangular truncations of our operators, which is a key difference to previousmethods that typically use square truncations. OUNDATIONS OF SPECTRAL COMPUTATIONS 11
Defining Ω and Λ . Let f be as described in (2), and ˆΩ to be the class of all A ∈ C ( l ( N )) such that(1) and (2) hold and such that the spectrum is non-empty. Given a sequence as described in (3), let Ω g be theclass of all A ∈ ˆΩ such that (3.10) holds. We also let Ω D denote the operators in ˆΩ that are diagonal. Operators on graphs.
For operators on graphs, consider any connected, undirected graph G , such the setof vertices V = V ( G ) is countably infinite. We consider operators on l ( V ) that are closed, densely definedand of the form(3.11) A = (cid:88) v,w ∈ V α ( v, w ) | v (cid:105) (cid:104) w | , for some α : V × V → C . We have also used the classical Dirac notation in (3.11) and identified any v ∈ V by the element in ψ v ∈ l ( V ) such that ψ v ( v ) = 1 and ψ v ( w ) = 0 for w (cid:54) = v . When writing this, weassume that the linear span of such vectors forms a core of both A and its adjoint. We also assume that forany v ∈ V , the set of vertices w with α ( v, w ) (cid:54) = 0 or α ( w, v ) (cid:54) = 0 is finite. We then let Ω G the class of allsuch A with non-empty spectrum and Ω G g operators in Ω G of known { g m } such that (3.10) holds. We alsoassume that with respect to some given enumeration v , v , ... of V , we have access to a function S : N → N such that if m > S ( n ) then α ( v n , v m ) = α ( v m , v n ) = 0 . Remark 3.7 (Defining Λ ) . For operators on l ( N ) , Λ contains the collection of matrix value evaluationfunctions, the functions describing the dispersion and the family of the functions { g m } controlling the growthof the resolvent. For operators on l ( V ) , Λ contains the functions α , the function S and, in the case of Ω G g ,the family { g m } .We can now state our main result in this section: Theorem 3.8 (Unbounded operators on graphs) . Let Ξ be the problem function Sp( · ) and Ξ be the problemfunction Sp (cid:15) ( · ) for (cid:15) > , where these map into the metric space (Cl( C ) , d AW ) . Then ∆ G (cid:54)(cid:51) { Ξ , Ω D } ∈ Σ A , ∆ G (cid:54)(cid:51) { Ξ , Ω g } ∈ Σ A , ∆ G (cid:54)(cid:51) { Ξ , Ω G g } ∈ Σ A , and ∆ G (cid:54)(cid:51) { Ξ , Ω D } ∈ Σ A , ∆ G (cid:54)(cid:51) { Ξ , ˆΩ } ∈ Σ A , ∆ G (cid:54)(cid:51) { Ξ , Ω G } ∈ Σ A . Furthermore the routines
CompSpecUB and
PseudoSpecUB in Appendix A realise the sharp Σ A inclu-sions, and in the case of Ξ , the output is guaranteed to be inside the true pseudospectrum. The algorithm used to compute the pseudospectrum can be applied to cases where the spectrum or pseu-dospectrum are empty, and we provide a computational example of this below. Finally, we consider twodiscrete problems which also include the case when the spectrum may be empty. Let K be a non-empty andcompact set in C and denote the collection of such subsets by K ( C ) . Consider Ξ : ( A, K ) → “Is Sp( A ) ∩ K = ∅ ?” Ξ : ( A, K ) → “Is Sp (cid:15) ( A ) ∩ K = ∅ ?”More precisely, the information we consider available to the algorithms in the l ( N ) ( l ( V ( G )) ) case are thematrix elements of A (the functions α ), the dispersion function f and dispersion bounds { c n } (the finite sets S v ) and a sequence of finite sets K n ⊂ Q + i Q , with the property that d H ( K n , K ) ≤ − ( n +1) . The followingshows that these discrete problems are harder than computing the spectrum. Theorem 3.9 (Does a set intersect the spectrum/pseudospectrum?) . We have the following classifications for j = 3 , : ∆ G (cid:54)(cid:51) { Ξ j , ˆΩ × K ( C ) } ∈ Π A , ∆ G (cid:54)(cid:51) { Ξ j , Ω D × K ( C ) } ∈ Π A , ∆ G (cid:54)(cid:51) { Ξ j , Ω G × K ( C ) } ∈ Π A . The routines
TestSpec and
TestPseudoSpec in Appendix A, used for Ξ and Ξ respectively, realisethe sharp Π A classifications. Furthermore, the proof will make clear that the lower bounds also hold whenwe restrict the allowed compact sets to any fixed compact subset of R . Remark 3.10.
By considering singletons K = { z } , we can test whether a point lies in the spectrum orpseudospectrum. Even when restricting to such K , the proof shows that the classification remains the same.3.3. The spectral gap problem and classifications of the spectrum.
The spectral gap problem has a longtradition and is linked to many important conjectures and problems such as the Haldane conjecture [68] orthe Yang–Mills mass gap problem in quantum field theory [20]. The spectral gap question is fundamentalin physics, and in the seminal paper [44] it was shown that the spectral gap problem is undecidable whenconsidering the thermodynamic limit of finite-dimensional Hamiltonians.In this paper, we consider the general infinite-dimensional problem. The question can be formulated inthe following way. Let (cid:98) Ω SA be the set of all bounded below, self-adjoint operators A on l ( N ) , for which thelinear span of the canonical basis form a core of A (we do not assume A is bounded above) and such thatone of the two following cases occur:(1) The minimum of the spectrum, a , is an isolated eigenvalue with multiplicity one.(2) There is some (cid:15) > such that [ a, a + (cid:15) ] ⊂ Sp( A ) .In the former case, we say the spectrum is gapped, whereas in the latter we say it is gappless. Note thatbecause we have restricted ourselves to the class where either (1) or (2) must hold, our problem is well-defined as a decision problem. Moreover, this definition is in line with the definitions in [44] and the physicsliterature. We also let (cid:98) Ω D denote the operators in (cid:98) Ω SA that are diagonal and define(3.12) Ξ gap : (cid:98) Ω SA , (cid:98) Ω D (cid:51) A (cid:55)→ “Is the spectrum of A gapped?”The above spectral gap problem can also be extended as follows. Let (cid:101) Ω f SA denote the class of operatorsthat are bounded below, self-adjoint, for which the linear span of the canonical basis form a core, and thathave (known) bounded dispersion with respect to the function f . Let a ( A ) = inf { x : x ∈ Sp( A ) } andconsider the following four cases(1) a ( A ) lies in the discrete spectrum and has multiplicity ,(2) a ( A ) lies in the discrete spectrum and has multiplicity > ,(3) a ( A ) lies in the essential spectrum but is an isolated point of the spectrum,(4) a ( A ) is a cluster point of Sp( A ) .We will consider the classification problem Ξ class which maps (cid:101) Ω f SA (or relevant subclasses) to the discretespace { , , , } (with the natural ordering). We denote by (cid:101) Ω D the class of diagonal operators in (cid:101) Ω f SA . Theorem 3.11 (Spectral gap and classification) . Let Ξ gap be as in (3.12) and (cid:98) Ω SA , (cid:98) Ω D as above. Similarly,Let Ξ class , (cid:101) Ω f SA and (cid:101) Ω D be as above. Then ∆ G (cid:54)(cid:51) { Ξ gap , (cid:98) Ω SA } ∈ Σ A , ∆ G (cid:54)(cid:51) { Ξ gap , (cid:98) Ω D } ∈ Σ A . In particular, the routine
SpecGap in Appendix A realises the sharp Σ A inclusions. Moreover, ∆ G (cid:54)(cid:51) { Ξ class , (cid:101) Ω f SA } ∈ Π A , ∆ G (cid:54)(cid:51) { Ξ class , (cid:101) Ω D } ∈ Π A , and SpecClass in Appendix A realises the sharp Π A inclusions. Remark 3.12 (Diagonal vs. full matrix) . It is worth noting that Theorem 3.11 shows that there is no differ-ence in the classification of the spectral gap problem between the set of diagonal matrices and the collectionof full matrices.3.4.
Computing discrete spectra, multiplicities and approximate eigenvectors.
Let Ω d N denote the classof bounded normal operators on l ( N ) with (known) bounded dispersion and with non-empty discrete spec-trum, and denote by Ω d D the class of bounded diagonal self-adjoint operators in Ω d N . For a normal operator A ,there is a simple decomposition of Sp( A ) into the discrete spectrum and the essential spectrum, denoted by Sp d ( A ) and Sp ess ( A ) respectively. The discrete spectrum consists of isolated points of the spectrum that are OUNDATIONS OF SPECTRAL COMPUTATIONS 13 eigenvalues of finite multiplicity. The essential spectrum has numerous definitions in the non-normal case,but for the normal case is defined as the set of z such that A − zI is not a Fredholm operator. Define theproblem function(3.13) Ξ d : Ω d N , Ω d D (cid:51) A (cid:55)→ Sp d ( A ) . We have taken the closure and restricted to operators with non-empty discrete spectrum since we wantconvergence with respect to the Hausdorff metric. However, the algorithm we build, Γ n ,n , has the propertythat lim n →∞ Γ n ,n ( A ) ⊂ Sp d ( A ) , so this is not restrictive in practice.Let also Ω f N denote the class of bounded normal operators with (known) bounded dispersion with respectto the function f . In addition, let Ω D denote the class of bounded self-adjoint diagonal operators and considerthe following discrete problem function(3.14) Ξ d : Ω f N , Ω D (cid:51) A (cid:55)→ “Is Sp d ( A ) (cid:54) = ∅ ?” Theorem 3.13.
Let Ξ d , Ω d N and Ω d D , as well as Ξ d , Ω f N and Ω D , be as above. Then, ∆ G (cid:54)(cid:51) { Ξ d , Ω d N } ∈ Σ A , ∆ G (cid:54)(cid:51) { Ξ d , Ω d D } ∈ Σ A and ∆ G (cid:54)(cid:51) { Ξ d , Ω f N } ∈ Σ A , ∆ G (cid:54)(cid:51) { Ξ d , Ω D } ∈ Σ A . In particular, the routines
DiscreteSpec and
DiscSpecEmpty in Appendix A realise the sharp Σ A inclusions for Ξ d and Ξ d respectively. The constructed algorithm Γ n ,n (routine DiscreteSpec ) has the property that given A ∈ Ω d N and z ∈ Sp d ( A ) , the following holds. If (cid:15) > is such that Sp( A ) ∩ B (cid:15) ( z ) = { z } , then there is at most one pointin Γ n ,n ( A ) that also lies in B (cid:15) ( z ) . Furthermore, the limit lim n →∞ Γ n ,n ( A ) = Γ n ( A ) is containedin the discrete spectrum and increases to Sp d ( A ) in the Hausdorff metric. In other words, a given point of Sp d ( A ) has at most one point in Γ n ,n ( A ) approximating it.We also want to compute multiplicities. Suppose that we have z n ,n ∈ Γ n ,n ( A ) with lim n →∞ z n ,n = z n = z ∈ Sp d ( A ) (the limit independent of n ). Our tower also computes a function h n ,n ( A, · ) over the output Γ n ,n ( A ) such that lim n →∞ lim n →∞ h n ,n ( A, z n ,n ) = h ( A, z ) (where h ( A, z ) denotes the multiplicity of the eigenvalue z ) in Z ≥ with the discrete metric. The routine Multiplicity in Appendix A computes h n ,n .Finally, ApproxEigenvector in Appendix A approximates eigenvectors. For simplicity we shall stickto eigenspaces of multiplicity , but note that these ideas can be easily extended to higher multiplicities toapproximate the whole eigenspace. The question is whether, given a z n in the output Γ n ,n ( A ) of thealgorithm DiscreteSpec and an approximation(3.15) σ ( P f ( n ) ( A − z n I ) | P n H ) ≤ E ( n , z n ) , where σ denotes the smallest singular value, we can find a x n of unit norm satisfying (cid:107) ( A − z n I ) x n (cid:107) ≤ E ( n , z n )+ c n (recall that c n is the dispersion bound). The discussion in § Theorem 3.14.
Suppose A ∈ Ω f N . Let δ > and z n ∈ Γ n ,n ( A ) such that z n → z ∈ Sp d ( A ) . Supposewe also have the computed bound (3.15), then we can compute a corresponding vector x n (of finite support)satisfying (cid:107) ( A − z n I ) x n (cid:107) < (cid:107) x n (cid:107) (cid:0) E ( n , z n ) + c n + δ (cid:1) and − δ < (cid:107) x n (cid:107) < δ in finitely many arithmetic operations. What happens when we cannot bound the dispersion?
Whilst Theorem 3.13 shows that computingthe discrete spectrum requires two limits, the constructed algorithm is still useful since the constructed algo-rithm has lim n →∞ Γ n ,n ( A ) ⊂ Sp d ( A ) . This is reflected in Theorem 3.14, which shows that we can stilleffectively approximate eigenspaces with error control. But what happens if we do not know a dispersionfunction f as in (3.9)? To investigate this case we let Ω d denote the class of bounded normal operators withnon-empty discrete spectrum and Ω d the class of bounded normal operators. As the next theorem reveals,we get a jump in the SCI hierarchy. Theorem 3.15.
Let Ξ di and Ω di be as above. Then, ∆ G (cid:54)(cid:51) { Ξ d , Ω d } ∈ Σ A and ∆ G (cid:54)(cid:51) { Ξ d , Ω d } ∈ Σ A .
4. C
ONNECTION TO P REVIOUS W ORK
The SCI hierarchy:
Our paper is part of the program on the SCI hierarchy [8, 42, 43, 73], which isa direct continuation of S. Smale’s work and his program on the foundations of computational mathemat-ics [14,15,101,103]. Other results in this program related to our paper are the results by C. McMullen [83,84]and P. Doyle & C. McMullen [53] on polynomial root-finding, that are classification results in the SCI hier-archy, and the contributions by L. Blum, F. Cucker, M. Shub & S. Smale [14, 15, 45, 97]. Further examplesare the results by C. Fefferman and L. Seco [55–63], proving the Dirac-Schwinger conjecture on the asymp-totical behaviour of the ground state energy of a family of Schr¨odinger operators that implicitly prove Σ A classifications in the SCI hierarchy. This is also the case in T. Hales’ Flyspeck program [69, 70] leading tothe proof of Kepler’s conjecture (Hilbert’s 18th problem) which also implicitly proves Σ A classifications. Itshould also be noted that many other problems in the foundations of computations such as the work by S.Weinberger [110], can be viewed in the context of the SCI hierarchy. Classical results on computing spectra:
Due to the vast literature on spectral computation, we can onlycite a small subset here that are related to this paper. The ideas of using computational and algorithmicapproaches to obtain spectral information date back to leading physicists and mathematicians such as H.Goldstine [67], T. Kato [77], F. Murray [67], E. Schr¨odinger [93], J. Schwinger [94] and J. von Neumann[67]. Schwinger introduced finite-dimensional approximations to quantum systems in infinite-dimensionalspaces that allow for spectral computations. An interesting observation is that Schwinger’s ideas were alreadypresent in the work of H. Weyl [112]. The work by H. Goldstine, F. Murray and J. von Neumann [67]was one of the first to establish rigorous convergence results, and their work yields ∆ A classification forcertain self-adjoint finite-dimensional problems. In [51] T. Digernes, V. S. Varadarajan and S. R. S. Varadhanproved convergence of spectra of Schwinger’s finite-dimensional discretisation matrices for a specific class ofSchr¨odinger operators with certain types of potential, which yields a ∆ A classification in the SCI hierarchy.The finite-section method, which has been intensely studied for spectral computation, and has often beenviewed in connection with Toeplitz theory, is very similar to Schwinger’s idea of approximating in a finite-dimensional subspace. The reader may want to consult the pioneering work by A. B¨ottcher [21, 22] andA. B¨ottcher & B. Silberman [26, 27]. W. Arveson [1–5] and N. Brown [28–30] pioneered the combinationof spectral computation and the C ∗ -algebra literature (which dates back to the work of A. B¨ottcher & B.Silberman [25]), both for the general spectral computation problem as well as for Schr¨odinger operators. Seealso the work by N. Brown, K. Dykema, and D. Shlyakhtenko [31], where variants of finite section analysisare implicitly used. Arveson also considered spectral computation in terms of densities, which is related toSzeg¨o’s work [104] on finite section approximations. Similar results are also obtained by A. Laptev and Y.Safarov [78]. Typically, when applied to appropriate subclasses of operators, finite section approaches yield ∆ A classification results. There are also other approaches based on the infinite QR algorithm in connectionwith Toda flows with infinitely many variables pioneered by P. Deift, L. C. Li, and C. Tomei [50]. See alsothe work by P. Deift, J. Demmel, C. Li, and C. Tomei [49]. E. B. Davies considered second order spectra OUNDATIONS OF SPECTRAL COMPUTATIONS 15 methods [46, 48], and E. Shargorodsky [96] demonstrated how second order spectra methods [46] will neverrecover the whole spectrum.
Recent results on computing spectra:
There are many recent directions in computational spectral theorythat are related to our work.(i)
Infinite-dimensional numerical linear algebra:
S. Olver, A.Townsend and M. Webb have provideda foundational and practical framework for infinite-dimensional numerical linear algebra and foun-dational results on computations with infinite data structures [86–89, 109]. This includes efficientcodes as well as theoretical results. The infinite-dimensional QL and QR algorithms, inspired by thework of Deift et. al. [49,50] mentioned above, are important parts of this program that yield classifi-cations in the SCI hierarchy of computing extreme elements in the spectrum, see also [42,71] for theinfinite-dimensional QR algorithm. The recent work of M. Webb and S. Olver [109] on computingspectra of Jacobi operators is also formulated in the SCI hierarchy.(ii)
Finite section approaches:
In the cases where the finite section method works, it will typicallyyield ∆ A classifications in the SCI hierarchy, and occasionally ∆ A classifications, see for examplethe work by A. B¨ottcher, H. Brunner, A. Iserles & S. Nørsett [23], A. B¨ottcher, S. Grudsky & A.Iserles [24], H. Brunner, A. Iserles & S. Nørsett [32, 33], M. Marletta [81] and M. Marletta & R.Scheichl [82]. The latter papers also discuss the failure of the finite section approach for certainclasses of operators, see also [71, 72].(ii) Resonances:
We would like to mention the recent work by M. Zworski [115, 116] on computingresonances that can be viewed in terms of the SCI hierarchy. In particular, the computational ap-proach [116] is based on expressing the resonances as limits of non-self-adjoint spectral problems,and hence the SCI hierarchy is inevitable, see also [100]. The recent work of J. Ben-Artzi, M.Marletta & F. R¨osler [9] on computing resonances is also formulated in terms of the SCI hierarchy.(ii)
Computer-assisted proofs:
We have already mentioned the results by C. Fefferman and L. Seco[55–63] on computer-assisted proofs proving classification results in the SCI hierarchy. However,recent results using computer-assisted proofs in spectral theory also includes the work of M. Brown,M. Langer, M. Marletta, C. Tretter, & M. Wagenhofer [80] and S. B¨ogli, M. Brown, M. Marletta, C.Tretter & M. Wagenhofer [19].
Acknowledgements.
The authors would like to thank Percy Deift, Charlie Fefferman, Tom Hales, AriLaptev, Steve Smale, Maciej Zworski and Shmuel Weinberger for helpful discussions. MJC acknowledgessupport from the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/L016516/1.ACH acknowledges support from a Royal Society University Research Fellowship, the Leverhulme Prize2017, as well as the UK Engineering and Physical Sciences Research Council (EPSRC) grant EP/L003457/1.5. M
ATHEMATICAL P RELIMINARIES
In this section, we formally define the SCI hierarchy. We have already presented the definition of acomputational problem { Ξ , Ω , M , Λ } . The goal is to find algorithms that approximate the function Ξ . Moregenerally, the main pillar of our framework is the concept of a tower of algorithms, which is needed todescribe problems that need several limits in the computation. However, first one needs the definition of ageneral algorithm. Definition 5.1 (General Algorithm) . Given a computational problem { Ξ , Ω , M , Λ } , a general algorithm isa mapping Γ : Ω → M such that for each A ∈ Ω (i) there exists a (non-empty) finite subset of evaluations Λ Γ ( A ) ⊂ Λ ,(ii) the action of Γ on A only depends on { A f } f ∈ Λ Γ ( A ) where A f := f ( A ) , (iii) for every B ∈ Ω such that B f = A f for every f ∈ Λ Γ ( A ) , it holds that Λ Γ ( B ) = Λ Γ ( A ) . Note that the definition of a general algorithm is more general than the definition of a Turing machine[108] or a Blum–Shub–Smale (BSS) machine [14]. A general algorithm has no restrictions on the operationsallowed. 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) ensuresthat the algorithm consistently reads the information. With a definition of a general algorithm, we can definethe concept of towers of algorithms.
Definition 5.2 (Tower of Algorithms) . Given a computational problem { Ξ , Ω , M , Λ } , a tower of algorithmsof height k for { Ξ , Ω , M , Λ } is a family of sequences of functions Γ n k : Ω → M , Γ n k ,n k − : Ω → M , . . . , Γ n k ,...,n : Ω → M , where n k , . . . , n ∈ N and the functions Γ n k ,...,n at the lowest level of the tower are general algorithms inthe sense of Definition 5.1. Moreover, for every A ∈ Ω , Ξ( A ) = lim n k →∞ Γ n k ( A ) , Γ n k ,...,n j +1 ( A ) = lim n j →∞ Γ n k ,...,n j ( A ) j = k − , . . . , . In addition to a general tower of algorithms (defined above), we will focus on arithmetic towers.
Definition 5.3 (Arithmetic Tower) . Given a computational problem { Ξ , Ω , M , Λ } , where Λ is countable, wedefine the following: An arithmetic tower of algorithms of height k for { Ξ , Ω , M , Λ } is a tower of algorithmswhere the lowest functions Γ = Γ n k ,...,n : Ω → M satisfy the following: For each A ∈ Ω the mapping ( n k , . . . , n ) (cid:55)→ Γ n k ,...,n ( A ) = Γ n k ,...,n ( { A f } f ∈ Λ ) is recursive, and Γ n k ,...,n ( A ) is a finite string ofcomplex numbers that can be identified with an element in M . For arithmetic towers we let α = A Remark 5.4.
By recursive we mean the following. If f ( A ) ∈ Q (or Q + i Q ) for all f ∈ Λ , A ∈ Ω , and Λ is countable, then Γ n k ,...,n ( { A f } f ∈ Λ ) can be executed by a Turing machine [108], that takes ( n k , . . . , n ) as input, and that has an oracle tape consisting of { A f } f ∈ Λ . If f ( A ) ∈ R (or C ) for all f ∈ Λ , then Γ n k ,...,n ( { A f } f ∈ Λ ) can be executed by a BSS machine [14] that takes ( n k , . . . , n ) , as input, and that hasan oracle that can access any A f for f ∈ Λ .Given the definitions above we can now define the key concept, namely, the Solvability Complexity Index: Definition 5.5 (Solvability Complexity Index) . A computational problem { Ξ , Ω , M , Λ } is said to have Solv-ability Complexity Index
SCI(Ξ , Ω , M , Λ) α = k , with respect to a tower of algorithms of type α , if k is thesmallest 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 forsome n < ∞ , then we define SCI(Ξ , Ω , M , Λ) α = 0 . The type α may be General, or Arithmetic, denotedrespectively G and A. We may sometimes write SCI(Ξ , Ω) α to simplify notation when M and Λ are obvious. We will let
SCI(Ξ , Ω) A and SCI(Ξ , Ω) G denote the SCI with respect to an arithmetic tower and a generaltower, respectively. Note that a general tower means just a tower of algorithms as in Definition 5.2, wherethere are no restrictions on the mathematical operations. Thus, clearly SCI(Ξ , Ω) A ≥ SCI(Ξ , Ω) G . Thedefinition of the SCI immediately induces the SCI hierarchy: Definition 5.6 (The Solvability Complexity Index Hierarchy) . Consider a collection C of computationalproblems and let T be the collection of all towers of algorithms of type α for the computational problems in C . Define ∆ α := {{ Ξ , Ω } ∈ C | SCI(Ξ , Ω) α = 0 } ∆ αm +1 := {{ Ξ , Ω } ∈ C | SCI(Ξ , Ω) α ≤ m } , m ∈ N , as well as ∆ α := {{ Ξ , Ω } ∈ C | ∃ { Γ n } n ∈ N ∈ T s.t. ∀ A d (Γ n ( A ) , Ξ( A )) ≤ − n } . OUNDATIONS OF SPECTRAL COMPUTATIONS 17
When there is additional structure on the metric space, such as in the spectral case when one considersthe Attouch–Wets or the Hausdorff metric, one can extend the SCI hierarchy.
Definition 5.7 (The SCI Hierarchy (Attouch–Wets/Hausdorff metric)) . Given the set-up in Definition 5.6,and suppose in addition that ( M , d ) has the Attouch–Wets or the Hausdorff metric induced by another metricspace ( M (cid:48) , d (cid:48) ) , define, for m ∈ N , Σ α = Π α = ∆ α , Σ α = {{ Ξ , Ω } ∈ ∆ α | ∃ { Γ n } ∈ T , { X n ( A ) } ⊂ M s.t. Γ n ( A ) ⊂ M (cid:48) X n ( A ) , lim n →∞ Γ n ( A ) = Ξ( A ) , d ( X n ( A ) , Ξ( A )) ≤ − n ∀ A ∈ Ω } , Π α = {{ Ξ , Ω } ∈ ∆ α | ∃ { Γ n } ∈ T , { X n ( A ) } ⊂ M s.t. Ξ( A ) ⊂ M (cid:48) X n ( A ) , lim n →∞ Γ n ( A ) = Ξ( A ) , d ( X n ( A ) , Γ n ( A )) ≤ − n ∀ A ∈ Ω } , where ⊂ M (cid:48) means inclusion in the metric space M (cid:48) , and { X n ( A ) } is a sequence where X n ( A ) ∈ M depends on A . Moreover, Σ αm +1 = {{ Ξ , Ω } ∈ ∆ αm +2 | ∃ { Γ n m +1 ,...,n } ∈ T , { X n m +1 ( A ) } ⊂ M s.t. Γ n m +1 ( A ) ⊂ M (cid:48) X n m +1 ( A ) , lim n m +1 →∞ Γ n m +1 ( A ) = Ξ( A ) , d ( X n m +1 ( A ) , Ξ( A )) ≤ − n m +1 ∀ A ∈ Ω } , Π αm +1 = {{ Ξ , Ω } ∈ ∆ αm +2 | ∃ { Γ n m +1 ,...,n } ∈ T , { X n m +1 ( A ) } ⊂ M s.t. Ξ( A ) ⊂ M (cid:48) X n m +1 ( A ) , lim n m +1 →∞ Γ n m +1 ( A ) = Ξ( A ) , d ( X n m +1 ( A ) , Γ n m +1 ( A )) ≤ − n m +1 ∀ A ∈ Ω } , where d can be either d H or d AW . Note that to build a Σ algorithm, it is enough by taking subsequences of n to construct Γ n ( A ) such that Γ n ( A ) ⊂ N E n ( A ) (Ξ( A )) with some computable E n ( A ) that converges to zero. The same extension canbe applied to the real line with the usual metric, or { , } with the discrete metric (where we interpret as“Yes”). Definition 5.8 (The SCI Hierarchy (totally ordered set)) . Given the set-up in Definition 5.6 and suppose inaddition that M is a totally ordered set. Define Σ α = Π α = ∆ α , Σ α = {{ Ξ , Ω } ∈ ∆ α | ∃ { Γ n } ∈ T s.t. Γ n ( A ) (cid:37) Ξ( A ) ∀ A ∈ Ω } , Π α = {{ Ξ , Ω } ∈ ∆ α | ∃ { Γ n } ∈ T s.t. Γ n ( A ) (cid:38) Ξ( A ) ∀ A ∈ Ω } , where (cid:37) and (cid:38) denotes convergence from below and above respectively, as well as, for m ∈ N , Σ αm +1 = {{ Ξ , Ω } ∈ ∆ αm +2 | ∃ { Γ n m +1 ,...,n } ∈ T s.t. Γ n m +1 ( A ) (cid:37) Ξ( A ) ∀ A ∈ Ω } , Π αm +1 = {{ Ξ , Ω } ∈ ∆ αm +2 | ∃ { Γ n m +1 ,...,n } ∈ T s.t. Γ n m +1 ( A ) (cid:38) Ξ( A ) ∀ A ∈ Ω } . Remark 5.9 ( ∆ α (cid:40) Σ α (cid:40) ∆ α ) . Note that the inclusions are strict. For example, if Ω K consists of the set ofcompact infinite matrices acting on l ( N ) and Ξ( A ) = Sp( A ) (the spectrum of A ) then { Ξ , Ω K } ∈ ∆ α butnot in Σ α ∪ Π α for α representing either towers of arithmetical or general type (see [8] for a proof). Moreover,as was demonstrated in [43], if ˜Ω is the set of discrete Schr¨odinger operators on l ( Z ) , then { Ξ , Ω } ∈ Σ α but not in ∆ α .Suppose we are given a computational problem { Ξ , Ω , M , Λ } , and that Λ = { f j } j ∈ β , where β is someindex set that can be finite or infinite. However, obtaining f j may be a computational task on its own, whichis exactly the problem in most areas of computational mathematics. In particular, for A ∈ Ω , f j ( A ) could bethe number e πj i for example. Hence, we cannot access f j ( A ) , but rather f j,n ( A ) where f j,n ( A ) → f j ( A ) as n → ∞ . Or, just as for problems that are high up in the SCI hierarchy, it could be that we need several limits, in particular one may need mappings f j,n m ,...,n : Ω → D + i D , where D denotes the dyadic rationalnumbers, such that(5.1) lim n m →∞ . . . lim n →∞ (cid:107){ f j,n m ,...,n ( A ) } j ∈ β − { f j ( A ) } j ∈ β (cid:107) ∞ = 0 ∀ A ∈ Ω . In particular, we may view the problem of obtaining f j ( A ) as a problem in the SCI hierarchy, where ∆ classification would correspond to the existence of mappings f j,n : Ω → D + i D such that(5.2) (cid:107){ f j,n ( A ) } j ∈ β − { f j ( A ) } j ∈ β (cid:107) ∞ ≤ − n ∀ A ∈ Ω . This idea is formalised in the following definition.
Definition 5.10 ( ∆ m -information) . Let { Ξ , Ω , M , Λ } be a computational problem. For m ∈ N we saythat Λ has ∆ m +1 -information if each f j ∈ Λ is not available, however, there are mappings f j,n m ,...,n :Ω → D + i D such that (5.1) holds. Similarly, for m = 0 there are mappings f j,n : Ω → D + i D suchthat (5.2) holds. Finally, if k ∈ N and ˆΛ is a collection of such functions described above such that Λ has ∆ k -information, we say that ˆΛ provides ∆ k -information for Λ . Moreover, we denote the family of all such ˆΛ by L k (Λ) . Note that we want to have algorithms that can handle all computational problems { Ξ , Ω , M , ˆΛ } when ˆΛ ∈ L m (Λ) . In order to formalise this, we define what we mean by a computational problem with ∆ m -information. Definition 5.11 (Computational problem with ∆ m -information) . Given m ∈ N , a computational problemwhere Λ has ∆ m -information is denoted by { Ξ , Ω , M , Λ } ∆ m := { ˜Ξ , ˜Ω , M , ˜Λ } , where ˜Ω = (cid:110) ˜ A = { f j,n m ,...,n ( A ) } j,n m ,...,n ∈ β × N m | A ∈ Ω , { f j } j ∈ β = Λ , f j,n m ,...,n satisfy (*) (cid:111) , and (*) denotes (5.1) if m > and (*) denotes (5.2) if m = 1 . Moreover, ˜Ξ( ˜ A ) = Ξ( A ) , and we have ˜Λ = { ˜ f j,n m ,...,n } j,n m ,...,n ∈ β × N m where ˜ f j,n m ,...,n ( ˜ A ) = f j,n m ,...,n ( A ) . Note that ˜Ξ is well-defined byDefinition 2.1 of a computational problem. The SCI and the SCI hierarchy, given ∆ m -information, is then defined in the standard obvious way. Wewill use the notation { Ξ , Ω , M , Λ } ∆ m ∈ ∆ αk to denote that the computational problem is in ∆ αk given ∆ m -information. When M and Λ are obvious then we will write { Ξ , Ω } ∆ m ∈ ∆ αk for short. Remark 5.12 (Classifications in this paper) . For the problems considered in this paper, the SCI classifi-cations do not change if we consider arithmetic towers with ∆ -information. This is easy to see throughChurch’s thesis and analysis of the stability of our algorithms. For example, we have been careful to restrictall relevant operations to Q rather than R and errors incurred from ∆ -information can be removed in the firstlimit. Explicitly, for the algorithms based on DistSpec (see Appendix A) it is possible to carry out an erroranalysis. We can also bound numerical errors (e.g. using interval arithmetic) and incorporate this uncertaintyfor the estimation of (cid:107) R ( z, A ) (cid:107) − and still gain the same classification of our problems. Similarly, for otheralgorithms based on similar functions. In other words, it does not matter which model of computation oneuses for a definition of ‘algorithm’, from a classification point of view they are equivalent for these spectralproblems. This leads to rigorous Σ αk or Π αk type error control suitable for verifiable numerics. In particular,for Σ α or Π α towers of algorithms, this could be useful for computer-assisted proofs.6. P ROOFS OF T HEOREMS ON U NBOUNDED O PERATORS ON G RAPHS
We will now prove the theorems in § § l ( N ) case. Given the graph G and enumeration v , v , ... of the vertices, consider the induced isomorphism l ( V ( G )) ∼ = l ( N ) . This induces a correspondingoperator on l ( N ) , where the functions α now become matrix values. For the lower bounds, we can consider OUNDATIONS OF SPECTRAL COMPUTATIONS 19 diagonal operators in Ω G (that is, α ( v, w ) = 0 if v (cid:54) = w ) with the trivial choice of S ( n ) = n . Hencelower bounds for Ω D translate to lower bounds for Ω G and Ω G g . For the upper bounds, the constructionof algorithms for l ( N ) will make clear that given the above isomorphism, we can compute a dispersionbounding function f for the induced operator on l ( N ) simply by taking f ( n ) = S ( n ) . This has D f,n ( A ) =0 . Note that any of the functions in Λ for the relevant class of operators on l ( N ) can be computed via theabove isomorphism using functions in Λ for the relevant class of operators on l ( V ( G )) . For instance, toevaluate matrix elements, we use α ( v i , v j ) .There is a useful characterisation of the Attouch–Wets topology. For any closed non-empty sets C and C n , the convergence d AW ( C n , C ) → holds if and only if d K ( C n , C ) → for any compact K ⊂ C where d K ( C , C ) = max (cid:26) sup a ∈ C ∩ K dist( a, C ) , sup b ∈ C ∩ K dist( b, C ) (cid:27) , with the convention that the supremum over the empty set is . This occurs if and only if for any δ > and K , there exists N such that if n > N then C n ∩ K ⊂ C + B δ (0) and C ∩ K ⊂ C n + B δ (0) . Furthermore,it is enough to consider K of the form B m (0) , the closed ball of radius m about the origin for m ∈ N , for m large. Throughout this section we take our metric space ( M , d ) to be (Cl( C ) , d AW ) . Remark 6.1 (A note on the empty set) . There is a slight subtlety regarding the empty set. It could be the casethat the output of our algorithm is the empty set and hence Γ n ( A ) does not map to the required metric space.However, the proofs will make clear that for large n , Γ n ( A ) is non-empty and we gain convergence. Bysuccessively computing Γ n ( A ) and outputting Γ m ( n ) ( A ) , where m ( n ) ≥ n is minimal with Γ m ( n ) ( A ) (cid:54) = ∅ ,we see that this does not matter for the classification, but the algorithm in this case is adaptive.The following lemma is a useful criterion for determining Σ A error control in the Attouch–Wets topologyand will be used in the proofs without further comment. Lemma 6.2.
Suppose that
Ξ : Ω → (Cl( C ) , d AW ) is a problem function and Γ n is a sequence of arithmeticalgorithms with each output a finite set such that lim n →∞ d AW (Γ n ( A ) , Ξ( A )) = 0 , ∀ A ∈ Ω . Suppose also that there is a function E n provided by Γ n (and defined over the output of Γ n ), such that lim n →∞ sup z ∈ Γ n ( A ) ∩ B m (0) E n ( z ) = 0 for all m ∈ N and such that dist( z, Ξ( A )) ≤ E n ( z ) , ∀ z ∈ Γ n ( A ) . Then:(1) For each m ∈ N and given Γ n ( A ) , we can compute in finitely many arithmetic operations andcomparisons a sequence of non-negative numbers a mn → (as n → ∞ ) such that Γ n ( A ) ∩ B m (0) ⊂ Ξ( A ) + B a mn (0) . (2) Given Γ n ( A ) , we can compute in finitely many arithmetic operations and comparisons a sequence ofnon-negative numbers b n → such that Γ n ( A ) ⊂ A n for some A n ∈ Cl( C ) with d AW ( A n , Ξ( A )) ≤ b n .Hence we can convert Γ n to a Σ A tower using the sequence { b n } by taking subsequences if necessary.Proof. For the proof of (1), we may take a mn = sup { E n ( z ) : z ∈ Γ n ( A ) ∩ B m (0) } and the result follows.Note that we need Γ n ( A ) to be finite to be able to compute this number with finitely many arithmeticoperations and comparisons. We next show (2) by defining A mn = (cid:0) (Ξ( A ) + B a mn (0)) ∩ B m (0) (cid:1) ∪ (Γ n ( A ) ∩ { z : | z | ≥ m } ) . It is clear that Γ n ( A ) ⊂ A mn and given Γ n ( A ) we can easily compute a lower bound m such that Ξ( A ) ∩ B m (0) (cid:54) = ∅ . Compute this from Γ ( A ) and then fix it. Suppose that m ≥ m , and suppose that | z | < (cid:98) m/ (cid:99) . Then the points in A mn and Ξ( A ) nearest to z must lie in B m (0) and hence dist( z, A mn ) ≤ dist( z, Ξ( A )) , dist( z, Ξ( A )) ≤ dist( z, A mn ) + a mn . It follows that d AW ( A mn , Ξ( A )) ≤ a mn + 2 −(cid:98) m/ (cid:99) . We now choose a sequence m ( n ) such that setting A n = A m ( n ) n and b n = a m ( n ) n + 2 −(cid:98) m ( n ) / (cid:99) proves theresult. Clearly it is enough to ensure that b n is null. If n < m then set m ( n ) = 4 m , otherwise consider m ≤ k ≤ n . If such a k exists with a kn ≤ − k then let m ( n ) be the maximal such k and finally if no such k exists then set m ( n ) = 4 m . For a fixed m , a mn → as n → ∞ . It follows that for large n , a m ( n ) n ≤ − m ( n ) and that m ( n ) → ∞ . (cid:3) Remark 6.3.
We will only consider algorithms where the output of Γ n ( A ) is at most finite for each n . Hencethe above restriction does not matter in what follows.In order to build our algorithms, we will need to characterise the reciprocal of resolvent norm in terms ofthe injection modulus. For A ∈ C ( l ( N )) define the injection modulus as(6.1) σ ( A ) = inf {(cid:107) Ax (cid:107) : x ∈ D ( A ) , (cid:107) x (cid:107) = 1 } , and define the function γ ( z, A ) = min { σ ( A − zI ) , σ ( A ∗ − ¯ zI ) } . Lemma 6.4.
For A ∈ C ( l ( N )) , γ ( z, A ) = 1 / (cid:107) R ( z, A ) (cid:107) , where R ( z, A ) denotes the resolvent ( A − zI ) − and we adopt the convention that / (cid:107) R ( z, A ) (cid:107) = 0 if z ∈ Sp( A ) .Proof. We deal with the case z / ∈ Sp( A ) first, where we prove σ ( A − zI ) = σ ( A ∗ − ¯ zI ) = 1 / (cid:107) R ( z, A ) (cid:107) .We show this for σ ( A − zI ) and the other case is similar using the fact that R ( z, A ) ∗ = R ( z, A ∗ ) and (cid:107) R ( z, A ) (cid:107) = (cid:107) R ( z, A ) ∗ (cid:107) . Let x ∈ D ( A ) with (cid:107) x (cid:107) = 1 then (cid:107) R ( z, A )( A − zI ) x (cid:107) ≤ (cid:107) R ( z, A ) (cid:107) (cid:107) ( A − zI ) x (cid:107) and hence upon taking infinum, σ ( A − zI ) ≥ / (cid:107) R ( z, A ) (cid:107) . Conversely, let x n ∈ l ( N ) such that (cid:107) x n (cid:107) = 1 and (cid:107) R ( z, A ) x n (cid:107) → (cid:107) R ( z, A ) (cid:107) . It follows that (cid:107) ( A − zI ) R ( z, A ) x n (cid:107) ≥ σ ( A − zI ) (cid:107) R ( z, A ) x n (cid:107) . Letting n → ∞ we get σ ( A − zI ) ≤ / (cid:107) R ( z, A ) (cid:107) .Now suppose that z ∈ Sp( A ) . If at least one of A − zI or A ∗ − ¯ zI is not injective on their respectivedomain then we are done, so assume both are one to one. Suppose also that σ ( A − zI ) , σ ( A ∗ − ¯ zI ) > otherwise we are done. It follows that R ( A − zI ) is dense in l ( N ) by injectivity of A ∗ − ¯ zI since R ( A − zI ) ⊥ = N ( A ∗ − ¯ zI ) . It follows that we can define ( A − zI ) − , bounded on the dense set R ( A − zI ) .We can extend this inverse to a bounded operator on the whole of l ( N ) . Closedness of A now implies that ( A − zI )( A − zI ) − = I . Clearly ( A − zI ) − ( A − zI ) x = x for all x ∈ D ( A ) . Hence, ( A − zI ) − = R ( z, A ) ∈ B ( l ( N )) so that z / ∈ Sp( A ) , a contradiction. (cid:3) Suppose we have a sequence of functions γ n ( z, A ) that converge uniformly to γ ( z, A ) on compact subsetsof C . Define the grid(6.2) Grid ( n ) = 1 n ( Z + i Z ) ∩ B n (0) . For an strictly increasing continuous function g : R ≥ → R ≥ , with g (0) = 0 and diverging at infinity, for n ∈ N and y ∈ R ≥ define(6.3) CompInvg ( n, y, g ) = min { k/n : k ∈ N , g ( k/n ) > y } . OUNDATIONS OF SPECTRAL COMPUTATIONS 21
Note that
CompInvg ( n, y, g ) can be computed from finitely many evaluations of the function g . We nowbuild the algorithm converging to the spectrum step by step using the functions in (3.10). For each z ∈ Grid ( n ) , let Υ n,z = B CompInvg ( n,γ n ( z,A ) ,g (cid:100)| z |(cid:101) ) ( z ) ∩ Grid ( n ) . If γ n ( z, A ) > (cid:16) | z | + 1 (cid:17) − then set M z = ∅ , otherwise set M z = { w ∈ Υ n,z : γ n ( w, A ) = min v ∈ Υ n,z γ n ( v, A ) } . Finally define Γ n ( A ) = ∪ z ∈ Grid ( n ) M z . It is clear that if γ n ( z, A ) can be computed in finitely many arith-metic operations and comparisons from the relevant functions in Λ for each problem, then this defines anarithmetic algorithm. If A ∈ C ( l ( N )) with non-empty spectrum then there exists z ∈ B m (0) with γ ( z, A ) ≤ ( m + 1) − / and, for large n , z n ∈ Grid ( n ) sufficiently close to z with γ ( z n , A ) ≤ ( | z n | + 1) − . Hence,by computing successive Γ n ( A ) , we can assume that Γ n ( A ) (cid:54) = ∅ without loss of generality (see Remark6.1). Proposition 6.5.
Suppose A ∈ C ( l ( N )) with non-empty spectrum and we have a function γ n ( z, A ) thatconverges uniformly to γ ( z, A ) on compact subsets of C . Suppose also that (3.10) holds, namely g m (dist( z, Sp( A ))) ≤ (cid:107) R ( z, A ) (cid:107) − , ∀ z ∈ B m (0) . Then Γ n ( A ) converges in the Attouch–Wets topology to Sp( A ) (assuming Γ n ( A ) (cid:54) = ∅ without loss of gener-ality).Proof. We use the characterisation of the Attouch–Wets topology. Suppose that m ∈ N is large such that B m (0) ∩ Sp( A ) (cid:54) = ∅ . We must show that given δ > , there exists N such that if n > N then Γ n ( A ) ∩ B m (0) ⊂ Sp( A ) + B δ (0) and Sp( A ) ∩ B m (0) ⊂ Γ n ( A ) + B δ (0) . Throughout the rest of the proof we fixsuch an m . Let (cid:15) n = (cid:107) γ n ( · , A ) − γ ( · , A ) (cid:107) ∞ ,B m +1 (0) , where the notation means the supremum norm overthe set B m +1 (0) .We deal with the second inclusion first. Suppose that z ∈ Sp( A ) ∩ B m (0) , then there exists some w ∈ Grid ( n ) such that | w − z | ≤ /n . It follows that γ n ( w, A ) ≤ γ ( w, A ) + (cid:15) n ≤ dist( w, Sp( A )) + (cid:15) n ≤ (cid:15) n + 1 /n. By choosing n large, we can ensure that (cid:15) n < (2 m + 2) − and that /n ≤ (2 m + 2) − so that γ n ( w, A ) < ( | w | + 1) − . It follows that M w is non-empty. If y ∈ M w then | y − z | ≤ | w − z | + | y − w | ≤ /n + 1 /n + g − (cid:100) w (cid:101) ( γ n ( w, A )) . But the g k ’s are non-increasing in k , strictly increasing continuous functions with g k (0) = 0 . Since γ n ( w, A ) ≤ (cid:15) n + 1 /n , it follows that(6.4) | y − z | ≤ /n + g − m +1 ( (cid:15) n + 1 /n ) . There exists N such that if n ≥ N then (6.4) holds and /n + g − m +1 ( (cid:15) n + 1 /n ) ≤ δ and this gives thesecond inclusion.For the first inclusion, suppose for a contradiction that this is false. Then there exists n j → ∞ , δ > and z n j ∈ Γ n j ( A ) ∩ B m (0) such that dist( z n j , Sp( A )) ≥ δ . Then z n j ∈ M w nj for some w n j ∈ Grid ( n j ) . Let I ( j ) = B CompInvg ( n j ,γ nj ( w nj ,A ) ,g (cid:100)| wnj |(cid:101) ) ( w n j ) ∩ Grid ( n j ) , the set that we compute minima of γ n j over. Let y n j ∈ Sp( A ) be of minimal distance to w n j (such a y n j exists since the spectrum restricted to any compact ball is compact). It follows that (cid:12)(cid:12) y n j − w n j (cid:12)(cid:12) ≤ g − (cid:100)| w nj |(cid:101) ( γ ( w n j , A )) . A simple geometrical argument (which also works when we restrict everything to thereal line for self-adjoint operators), shows that there must be a v n j in I ( j ) so that (cid:12)(cid:12) v n j − y n j (cid:12)(cid:12) ≤ n j + g − (cid:100)| w nj |(cid:101) ( γ ( w n j , A )) − g − (cid:100)| w nj |(cid:101) ( γ n j ( w n j , A )) . Since z n j minimises γ n j over I ( j ) and M w nj is non-empty, it follows that γ ( z n j , A ) ≤ γ n j ( z n j , A ) + (cid:15) n j ≤ min (cid:40) (cid:12)(cid:12) w n j (cid:12)(cid:12) + 1 , γ n j ( v n j , A ) (cid:41) + (cid:15) n j . This implies that(6.5) δ ≤ dist( z n j , Sp( A )) ≤ g − m (cid:32) min (cid:40) (cid:12)(cid:12) w n j (cid:12)(cid:12) + 1 , γ n j ( v n j , A ) (cid:41) + (cid:15) n j (cid:33) , where we recall that g − m is continuous. It follows that the w n j must be bounded and hence so are the v n j .Due to the local uniform convergence of γ n to γ , it follows that n j + g − (cid:100)| w nj |(cid:101) ( γ ( w n j , A )) − g − (cid:100)| w nj |(cid:101) ( γ n j ( w n j , A )) → , as n j → ∞ . But then γ ( v n j , A ) ≤ dist( v n j , Sp( A )) ≤ (cid:12)(cid:12) v n j − y n j (cid:12)(cid:12) → . Again the local uniform convergence implies that γ n j ( v n j , A ) → , which contradicts (6.5) and completesthe proof. (cid:3) Next, given such a sequence γ n , we would like to provide an algorithm for computing the pseudospectrum.However, care must be taken in the unbounded case since the resolvent norm can be constant on open subsetsof C [95]. Simply taking Grid ( n ) ∩ { z : γ n ( z, A ) ≤ (cid:15) } is not guaranteed to converge, as can be seen in thecase that γ n is identically γ and A is such that the level set {(cid:107) R ( z, A ) (cid:107) − = (cid:15) } has non-empty interior. Toget around this, we will need an extra assumption on the functions γ n . Lemma 6.6.
Suppose A ∈ C ( l ( N )) with non-empty spectrum and let (cid:15) > . Suppose we have a sequenceof functions γ n ( z, A ) that converge uniformly to (cid:107) R ( z, A ) (cid:107) − on compact subsets of C . Set Γ (cid:15)n ( A ) = Grid ( n ) ∩ { z : γ n ( z, A ) < (cid:15) } . For large n , Γ (cid:15)n ( A ) (cid:54) = ∅ so we can assume this without loss of generality. Suppose also ∃ N ∈ N (pos-sibly dependent on A but independent of z ) such that if n ≥ N then γ n ( z, A ) ≥ (cid:107) R ( z, A ) (cid:107) − . Then d AW (Γ (cid:15)n ( A ) , Sp (cid:15) ( A )) → as n → ∞ .Proof. Since the pseudospectrum is non-empty, for large n , Γ (cid:15)n ( A ) (cid:54) = ∅ so by our usual argument of comput-ing successive Γ (cid:15)n (see Remark 6.1) we may assume that this holds for all n without loss of generality. We usethe characterisation of the Attouch–Wets topology. Suppose that m is large such that B m (0) ∩ Sp (cid:15) ( A ) (cid:54) = ∅ . ∃ N ∈ N such that if n ≥ N then γ n ( z, A ) ≥ (cid:107) R ( z, A ) (cid:107) − and hence Γ (cid:15)n ( A ) ∩ B m (0) ⊂ Sp (cid:15) ( A ) . Hence wemust show that given δ > , there exists N such that if n > N then Sp (cid:15) ( A ) ∩ B m (0) ⊂ Γ (cid:15)n ( A ) + B δ (0) .Suppose for a contradiction that this were false. Then there exists z n j ∈ Sp (cid:15) ( A ) ∩ B m (0) , δ > and n j → ∞ such that dist( z n j , Γ (cid:15)n j ( A )) ≥ δ . Without loss of generality, we can assume that z n j → z ∈ Sp (cid:15) ( A ) ∩ B m (0) .There exists some w with (cid:107) R ( w, A ) (cid:107) − < (cid:15) and | z − w | ≤ δ/ . Assuming n j > m + δ , there exists y n j ∈ Grid ( n j ) with (cid:12)(cid:12) y n j − w (cid:12)(cid:12) ≤ /n j . It follows that γ n j ( y n j , A ) ≤ (cid:12)(cid:12) γ n j ( y n j , A ) − γ ( y n j , A ) (cid:12)(cid:12) + (cid:12)(cid:12) γ ( w, A ) − γ ( y n j , A ) (cid:12)(cid:12) + (cid:107) R ( w, A ) (cid:107) − . But γ is continuous and γ n j converges uniformly to γ on compact subsets. Hence for large n j , it followsthat γ n j ( y n j , A ) < (cid:15) so that y n j ∈ Γ (cid:15)n j ( A ) . But (cid:12)(cid:12) y n j − z (cid:12)(cid:12) ≤ | z − w | + (cid:12)(cid:12) y n j − w (cid:12)(cid:12) ≤ δ/ /n j , which issmaller than δ for large n j . This gives the required contradiction. (cid:3) OUNDATIONS OF SPECTRAL COMPUTATIONS 23
Now suppose that A ∈ ˆΩ and let D f,n ( A ) ≤ c n . The following shows that we can construct the requiredsequence γ n ( z, A ) , each function output requiring finitely many arithmetic operations and comparisons ofthe corresponding input information. Theorem 6.7.
Let A ∈ ˆΩ and define the function ˜ γ n ( z, A ) = min { σ ( P f ( n ) ( A − zI ) | P n ( l ( N )) ) , σ ( P f ( n ) ( A ∗ − ¯ zI ) | P n ( l ( N )) ) } . We can compute ˜ γ n up to precision /n using finitely many arithmetic operations and comparisons. We callthis approximation ˆ γ n and set γ n ( z, A ) = ˆ γ n ( z, A ) + c n + 1 /n. Then γ n ( z, A ) converges uniformly to γ ( z, A ) on compact subsets of C and γ n ( z, A ) ≥ γ ( z, A ) .Proof. We will first prove that σ (( A − zI ) | P n ( l ( N )) ) ↓ σ ( A − zI ) as n → ∞ . It is trivial that σ (( A − zI ) | P n ( l ( N )) ) ≥ σ ( A − zI ) and that σ (( A − zI ) | P n ( l ( N )) ) is non-increasing in n . Using Lemma 6.4, let (cid:15) > and x ∈ D ( A ) such that (cid:107) x (cid:107) = 1 and (cid:107) ( A − zI ) x (cid:107) ≤ σ ( A − zI ) + (cid:15) . Since span { e n : n ∈ N } forms a core of A , AP n j x n j → Ax and P n j x n j → x for some n j → ∞ and some sequence of vectors x n j that we can assume have norm . It follows that for large n j σ (( A − zI ) | P nj ( l ( N )) ) ≤ (cid:13)(cid:13) ( A − zI ) P n j x n j (cid:13)(cid:13)(cid:13)(cid:13) P n j x n j (cid:13)(cid:13) → (cid:107) ( A − zI ) x (cid:107) ≤ σ ( A − zI ) + (cid:15). Since (cid:15) > was arbitrary, this shows the convergence of σ (( A − zI ) | P n ( l ( N )) ) . The fact that span { e n : n ∈ N } forms a core of A ∗ can also be used to show that σ (( A − zI ) ∗ | P n ( l ( N )) ) ↓ σ ( A ∗ − zI ) .Next we will use the assumption of bounded dispersion. For any bounded operators B, C , it holds that | σ ( A ) − σ ( B ) | ≤ (cid:107) A − B (cid:107) . The definition of bounded dispersion now implies that (cid:12)(cid:12) ˜ γ n ( z, A ) − min { σ (( A − zI ) | P n ( l ( N )) ) , σ (( A − zI ) ∗ | P n ( l ( N )) ) } (cid:12)(cid:12) ≤ c n . The monotone convergence of min { σ (( A − zI ) | P n ( l ( N )) ) , σ (( A − zI ) ∗ | P n ( l ( N )) ) } , together with Dini’stheorem, imply that ˜ γ n ( z, A ) converges uniformly to the continuous function γ ( z, A ) on compact subsets of C with ˜ γ n ( z, A ) + c n ≥ γ ( z, A ) .The proof will be complete if we can show that we can compute ˜ γ n ( z, A ) to precision /n using finitelymany arithmetic operations and comparisons. To do this, consider the matrices B n ( z ) = P n ( A − zI ) ∗ P f ( n ) ( A − zI ) P n , C n ( z ) = P n ( A − zI ) P f ( n ) ( A − zI ) ∗ P n . By an interval search routine and Lemma 6.8 below, we can determine the smallest l ∈ N such that at leastone of B n ( z ) − ( l/n ) I or C n ( z ) − ( l/n ) I has a negative eigenvalue. We then output l/n to get the /n bound. (cid:3) Recall that every finite Hermitian matrix B (not necessarily positive definite) has a decomposition P BP T = LDL ∗ , where L is lower triangular with ’s along its diagonal, D is block diagonal with block sizes or and P is a permutation matrix. Furthermore, this decomposition can be computed with finitely many arithmeticoperations and comparisons. Throughout, we will assume without loss of generality that P is the identitymatrix. Lemma 6.8.
Let B ∈ C n be self-adjoint (Hermitian), then we can determine the number of negative eigen-values of B in finitely many arithmetic operations and comparisons (assuming no round-off errors) on thematrix entries of B . Proof.
We can compute the decomposition B = LDL ∗ in finitely many arithmetical operations and compar-isons. By Sylvester’s law of inertia (the Hermitian version), D has the same number of negative eigenvaluesas B . It is then clear that we only need to deal with × matrices corresponding to the maximum block sizeof D . Let λ , λ be the two eigenvalues of such a matrix, then we can determine their sign pattern from thetrace and determinant of the matrix. (cid:3) This lemma has a corollary that will be useful in § Corollary 6.9.
Let B ∈ C n be self-adjoint (Hermitian) and list its eigenvalues in increasing order, includingmultiplicity, as λ ≤ λ ≤ ... ≤ λ n . In exact arithmetic, given (cid:15) > , we can compute λ , λ , ...λ n toprecision (cid:15) using only finitely many arithmetic operations and comparisons.Proof. Consider A ( λ ) = B − λI . We will apply Lemma 6.8 to A ( λ ) for various λ . First by considering thesequences − , − , ... and , , ... we can find m ∈ N such that Sp( B ) ⊂ ( − m , m ) . Now let m ∈ N such that /m < (cid:15) and let a j be the output of Lemma 6.8 applied to A ( j/m ) for − m m ≤ j ≤ m m .Set ˜ λ k = min { j : − m m ≤ j ≤ m m , a j ≥ k } , k = 1 , ..., n. If λ k ∈ [ j/m , ( j + 1) /m ) then ˜ λ k = ( j + 1) /m and hence (cid:12)(cid:12)(cid:12) ˜ λ k − λ k (cid:12)(cid:12)(cid:12) ≤ /m < (cid:15) . (cid:3) Remark 6.10.
Of course, in practice, there are much more computationally efficient ways to numericallycompute eigenvalues or singular values - the above is purely used to show this can be done to any precisionwith finitely many arithmetic operations.Note that by taking successive minima, υ n ( z, A ) = min ≤ j ≤ n γ n ( z, A ) , we can obtain a sequence offunctions υ n that converge uniformly on compact subsets of C to γ ( z, A ) monotonically from above. Hencewithout loss of generality, we will always assume that γ n have this property. We can now prove our mainresult. Proof of Theorem 3.8.
By considering bounded diagonal operators, it is straightforward to see that none ofthe problems lie in ∆ G . We first deal with convergence of height one arithmetical towers. For the spectrum,we use the function γ n described in Theorem 6.7 together with Proposition 6.5 and its described algorithm.For the pseudospectrum, we use the same function γ n described in Theorem 6.7 and convergence followsfrom using the algorithm in Proposition 6.6.We are left with proving that our algorithms have Σ A error control. For any A ∈ ˆΩ , the output of the algo-rithm in Proposition 6.6 is contained in the true pseudospectrum since γ n ( z, A ) ≥ γ ( z, A ) = (cid:107) R ( z, A ) (cid:107) − .Hence we need only show that the algorithm in Proposition 6.5 provides Σ A error control for input A ∈ Ω g .Denote the algorithm by Γ n and set E n ( z ) = CompInvg ( n, γ n ( z, A ) , g − (cid:100)| z |(cid:101) ) on Γ n ( A ) and zero on C \ Γ n ( A ) . Since γ n ( z, A ) ≥ (cid:107) R ( z, A ) (cid:107) − , the assumptions on { g m } imply that dist( z, Sp( A )) ≤ E n ( z ) , ∀ z ∈ Γ n ( A ) . Suppose for a contradiction that E n does not converge uniformly to zero on compact subsets of C . Thenthere exists some compact set K , some (cid:15) > , a sequence n j → ∞ and z n j ∈ K such that E n j ( z n j ) ≥ (cid:15) .It follows that z n j ∈ Γ n j ( A ) . Without loss of generality, z n j → z . By convergence of Γ n j ( A ) , z ∈ Sp( A ) and hence γ n j ( z n j , A ) → γ ( z, A ) = 0 . Now choose M large such that K ⊂ B M (0) . But then E n j ( z n j ) ≤ g − M ( γ n j ( z n j , A )) + 1 n j → , the required contradiction. (cid:3) OUNDATIONS OF SPECTRAL COMPUTATIONS 25
Remark 6.11.
The above makes it clear that E n ( z ) converges uniformly to the function g − (cid:100)| z |(cid:101) ( γ ( z, A )) as n → ∞ on compact subsets of C .Finally, we consider the decision problems Ξ and Ξ . Proof of Theorem 3.9.
It is clearly enough to prove the lower bounds for Ω D × K ( C ) and the existence oftowers for ˆΩ × K ( C ) . The proof of lower bounds for Ω D × K ( C ) can also be trivially adapted to the morerestrictive versions of the problem described in the theorem. Step 1 : { Ξ , Ω D × K ( C ) } (cid:54)∈ ∆ G . Suppose this were false, and Γ n is a height one tower solving theproblem. For every A and n there exists a finite number N ( A, n ) ∈ N such that the evaluations from Λ Γ n ( A ) only take the matrix entries A ij = (cid:104) Ae j , e i (cid:105) with i, j ≤ N ( A, n ) into account. Without loss ofgenerality (by shifting our argument), we assume that K ∩ [0 ,
1] = { } . We will consider the operators A m = diag { , / , ..., /m } ∈ C m × m , B m = diag { , , ..., } ∈ C m × m and C = diag { , , ... } . Set A = (cid:76) ∞ m =1 ( B k m ⊕ A k m ) where we choose an increasing sequence k m inductively as follows.Set k = 1 and suppose that k , ..., k m have been chosen. Sp( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = { , / , ..., /m } and hence Ξ ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = “No” , so there exists some n m ≥ m such that if n ≥ n m then Γ n ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = “No” . Now let k m +1 ≥ max { N ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C, n m ) , k m +1 } . By assumption (iii) in Definition5.1 it follows that Λ Γ nm ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = Λ Γ nm ( A ) and hence by assumption (ii) inthe same definition that Γ n m ( A ) = Γ n m ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = “No”. But ∈ Sp( A ) andso must have lim n →∞ (Γ n ( A )) = “Yes”, a contradiction. Step 2 : { Ξ , Ω D } (cid:54)∈ ∆ G . The same proof as step 1, but replacing A by A + (cid:15)I works in this case. Step 3 : { Ξ , ˆΩ × K ( C ) } ∈ Π A . Recall that we can compute, with finitely many arithmetic operationsand comparisons, a function γ n that converges monotonically down to (cid:107) R ( z, A ) (cid:107) − uniformly on compacts.Set Γ n ,n ( A ) = “Does there exist some z ∈ K n such that γ n ( z, A ) < / n ?” . It is clear that this is an arithmetic algorithm since each K n is finite and that lim n →∞ Γ n ,n ( A ) = “Does there exist some z ∈ K n such that (cid:107) R ( z, A ) (cid:107) − < / n ?” =: Γ n ( A ) . If K ∩ Sp( A ) = ∅ , then (cid:107) R ( z, A ) (cid:107) − is bounded below on the compact set K and hence for large n , Γ n ( A ) = “No”. However, if z ∈ Sp( A ) ∩ K then let z n ∈ K n minimise the distance to z . Then (cid:107) R ( z n , A ) (cid:107) − ≤ dist( z n , Sp( A )) < / n and hence Γ n ( A ) = “Yes” for all n . This also shows the Π A classification. Step 4 : { Ξ , ˆΩ × K ( C ) } ∈ Π A . Set Γ n ,n ( A ) = “Does there exist some z ∈ K n such that γ n ( z, A ) < / n + (cid:15) ?” , then the same argument used in step 3 works in this case. (cid:3) Examples of f used in the computational examples. We end with some examples for the graph case l ( V ( G )) . Suppose our enumeration of the vertices obeys the following pattern. The neighbours of v (including itself) are S = { v , v , ..., v q } for some finite q . The set of neighbours of these vertices is S = { v , ..., v q } for some finite q where we continue the enumeration of S and this process continuesinductively enumerating S m . Example 6.12.
Suppose that the bounded operator A can be written as(6.6) A = (cid:88) v ∼ k w α ( v, w ) | v (cid:105) (cid:104) w | for some k ∈ N (we write v ∼ k w for two vertices v, w ∈ V if there is a path of at most k edges connecting v and w , that is, A only involves k -th nearest neighbour interactions). Suppose also that the vertex degree of G is bounded by M . It holds that v n ∈ S n and { w ∈ V : v ∼ k w } ⊂ S n + k . Inductively | S m | ≤ ( M + 1) m and hence we may take the upper bound S ( n ) = ( M + 1) n + k . Example 6.13.
Consider a nearest neighbour operator ( k = 1 in (6.6)) on l ( Z d ) . It holds that | S m | ∼O ( m d ) whilst | S m +1 − S m | ∼ O ( m d − ) (by considering radial spheres). It is easy to see that we canchoose a suitable S such that S ( n ) − n ∼ O ( n d − d ) , that is, S grows at most linearly.7. P ROOFS OF T HEOREMS ON D IFFERENTIAL O PERATORS ON U NBOUNDED D OMAINS
Here we shall prove Theorems 3.3 and 3.5. The constructed algorithms involve technical error estimateswith parameters depending on these estimates. In the construction of the algorithms, our strategy will beto reduce the problem to one handled by the proofs in §
6. In order to do so, we must first select a suitablebasis and then compute matrix values. Recall that our aim is to compute the spectrum and pseudospectrumfrom the information given to us regarding the functions a k and ˜ a k , with the information we can evaluatemade precise by the mappings Ξ j , Ξ j , Ξ j and Ξ j . We will start by constructing the algorithms used for thepositive results in Theorems 3.3 and 3.5 and then prove the lower bounds.7.1. Construction of algorithms.
We begin with the description for d = 1 and comment how this can easilybe extended to arbitrary dimensions. As an orthonormal basis of L ( R ) we choose the Hermite functions ψ m ( x ) = (2 m m ! √ π ) − / e − x / H m ( x ) , m ∈ Z ≥ , where H n denotes the n -th (physicists’) Hermite polynomial defined by H n ( x ) = ( − n exp( x ) d n dx n exp( − x ) . These obey the recurrence relations ψ (cid:48) m ( x ) = (cid:114) m ψ m − ( x ) − (cid:114) m + 12 ψ m +1 ( x ) xψ m ( x ) = (cid:114) m ψ m − ( x ) + (cid:114) m + 12 ψ m +1 ( x ) . We let C H ( R ) = span { ψ m : m ∈ Z ≥ } . Note that since the Hermite functions decay like e − x / (up topolynomials) and the functions a k and ˜ a k can only grow polynomially, the formal differential operator T and its formal adjoint T ∗ make sense as operators from C H ( R ) to L ( R ) . The next proposition says that wecan use the chosen basis. Proposition 7.1.
Consider an operator T ∈ Ω . Then C H ( R ) forms a core of both T and T ∗ .Proof. Let f ∈ C H ( R ) and choose φ ∈ C ∞ ( R ) (the space of compactly supported smooth functions)bounded by such that φ ( x ) = 1 for all | x | ≤ . It is straightforward using the fact that the a k ’s arepolynomially bounded to show that lim n →∞ φ ( x/n ) f ( x ) = f ( x ) , lim n →∞ T φ ( x/n ) f ( x ) = ( T f )( x ) OUNDATIONS OF SPECTRAL COMPUTATIONS 27 in L ( R ) , where T f is the formal differential operator applied to f . The fact that T is closed implies that f ∈ D ( T ) . Let ˜ T denote the closure of the formal operator T , acting on C H ( R ) , then we have shown that ˜ T exists with ˜ T ⊂ T . Hence to show that C H ( R ) forms a core of T , we must show that C ∞ ( R ) ⊂ D ( ˜ T ) . Let g ∈ C ∞ ( R ) then in the L sense write g = (cid:88) m ≥ b m ψ m . Define g n = (cid:80) nm =0 b m ψ m then, since ˜ T is closed, it is enough to show that ˜ T g n converges as n → ∞ . Let H denote the closure of the operator − d /dx + x with initial domain C ∞ ( R ) then Hψ m = (2 m + 1) ψ m and H is self-adjoint. Note also that g ∈ D ( H n ) for any n ∈ N . But (cid:104) Hg, ψ m (cid:105) = (2 m + 1) (cid:104) g, ψ m (cid:105) =(2 m + 1) b m , so { (2 m + 1) | b m |} is square summable. We can repeat this argument any number of times toget that the coefficients b m decay faster than any inverse polynomial. To prove the required convergence, itis enough to consider one of the terms a k ( x ) ∂ k that defines ˜ T acting on C H ( R ) . The coefficient a k ( x ) ispolynomially bounded almost everywhere, and for some A k and B k (cid:104) a k ∂ k ψ m , a k ∂ k ψ m (cid:105) ≤ A k (cid:90) R (1 + | x | B k ) ∂ k ψ m ( x ) ∂ k ψ m ( x ) dx. But we can use the recurrence relations for the derivatives of the Hermite functions and orthogonality tobound the right hand side by a polynomial in m . The convergence now follows since ˜ T g n is a Cauchysequence due to the rapid decay of the { b m } . Exactly the same argument works for T ∗ . (cid:3) Clearly, all of the above analysis holds in higher dimensions by considering tensor products C H ( R d ) := span { ψ m ⊗ ... ⊗ ψ m d | m , ..., m d ∈ Z ≥ } of Hermite functions. We will abuse notation and write ψ m = ψ m ⊗ ... ⊗ ψ m d . It will be clear from thecontext when we are dealing with the multidimensional case. In order to build the required algorithms with Σ A error control, we need to select an enumeration of Z d ≥ in order to represent T as an operator acting on l ( N ) . A simple way to do this is to consider successive half spheres S n = { m ∈ Z d ≥ : | m | ≤ n } . We list S as { e , ..., e r } and given an enumeration { e , ..., e r n } of S n , we list S n +1 \ S n as { e r n +1 , ..., e r n +1 } . We willthen list our basis functions as e , e , ... with ψ m = e h ( m ) . In practice, it is often more efficient (especiallyfor large d ) to consider other orderings such as the hyperbolic cross [79]. Now that we have a suitable basis,the next question to ask is how to recover the matrix elements of T . In § γ n ( z, T ) , which also converges uniformly from aboveto (cid:107) R ( z, T ) (cid:107) − on compact subsets of C . Such a sequence of functions is given by Ψ n ( z, T ) := min { σ (( T − zI ) | P n ( l ( N )) ) , σ (( T ∗ − ¯ zI ) | P n ( l ( N )) ) } as long as the linear span of the basis forms a core of T and T ∗ . In § γ n ( z, T ) , it suffices to use the following. Lemma 7.2.
Let (cid:15) > and n ∈ N , and suppose that we can compute, with finitely many arithmeticoperations and comparisons, the matrices { W n ( z ) } ij = (cid:104) ( T − zI ) e j , ( T − zI ) e i (cid:105) + E n, ij ( z ) { V n ( z ) } ij = (cid:104) ( T − zI ) ∗ e j , ( T − zI ) ∗ e i (cid:105) + E n, ij ( z ) for ≤ i, j ≤ n where the entrywise errors E n, i,j and E n, i,j have magnitude at most (cid:15) . Then (cid:12)(cid:12) Ψ n ( z, T ) − min { σ ( W n ) , σ ( V n ) } (cid:12)(cid:12) ≤ n(cid:15). It follows that if (cid:15) is known, we can compute Ψ n ( z, T ) to within n(cid:15) . If (cid:15) is unknown, then for any δ > ,we can compute Ψ n ( z, T ) to within n(cid:15) + δ . (In each case with finitely many arithmetic operations andcomparisons.) Proof.
Given { W n ( z ) } ij , ( { W n ( z ) } ij + { W n ( z ) } ji ) / still has an entrywise absolute error bounded by (cid:15) .Hence without loss of generality we can assume that the approximations W n ( z ) and V n ( z ) are self-adjoint.Call the matrices with no errors ˜ W n ( z ) and ˜ V n ( z ) then note that min { σ (( T − zI ) | P n ( l ( N )) ) , σ (( T ∗ − ¯ zI ) | P n ( l ( N )) ) } = min { σ ( ˜ W n ) , σ ( ˜ V n ) } and(7.1) (cid:12)(cid:12)(cid:12) min { σ ( ˜ W n ) , σ ( ˜ V n ) } − min { σ ( W n ) , σ ( V n ) } (cid:12)(cid:12)(cid:12) ≤ max (cid:110)(cid:13)(cid:13)(cid:13) W n − ˜ W n (cid:13)(cid:13)(cid:13) , (cid:13)(cid:13)(cid:13) V n − ˜ V n (cid:13)(cid:13)(cid:13)(cid:111) . But for a finite matrix M , we can bound (cid:107) M (cid:107) by its Frobenius norm (cid:113)(cid:80) | M ij | . Hence the right hand sideof (7.1) is at most n(cid:15) . In order to use finitely many arithmetic operations and comparisons, we note that givena self-adjoint positive semi-definite matrix M , we can compute σ ( M ) to arbitrary precision using finitelymany arithmetic operations and comparisons via the argument in the proof of Theorem 6.7. The lemma nowfollows. (cid:3) Finally, we will need some results from the subject of quasi-Monte Carlo numerical integration, whichwe use to build the algorithm. Note that with either no prior information concerning the coefficients or forlarge d , this is the type of approach one would use in practice. We start with some definitions and theoremswhich we include here for completeness. An excellent reference for these results is [85]. Definition 7.3.
Let { t , ..., t j } be a sequence in [0 , d and let K denote all subsets of [0 , d of the form (cid:81) dk =1 [0 , y k ) for y k ∈ (0 , . Then we define the star discrepancy of { t , ..., t j } to be D ∗ j ( { t , ..., t j } ) = sup K ∈K (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) j j (cid:88) k =1 χ K ( t j ) − | K | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , where χ K denotes the characteristic function of K . Theorem 7.4. If { t k } k ∈ N is the Halton sequence in [0 , d in the pairwise relatively prime bases q , ..., q d ,then D ∗ j ( { t , ..., t j } ) < dj + 1 j d (cid:89) k =1 (cid:18) q k −
12 log( q k ) log( j ) + q k + 12 (cid:19) . Note that given d (and suitable q ,..., q d ), we can easily compute in finitely many arithmetic operationsand comparisons a constant C ( d ) such that the above implies(7.2) D ∗ j ( { t , ..., t j } ) < C ( d ) (log( j ) + 1) d j . The following theorem says why this is useful.
Theorem 7.5 (Koksma–Hlawka inequality [85]) . If f has bounded variation TV [0 , d ( f ) on the hypercube [0 , d then for any t , ..., t j in [0 , d (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) j j (cid:88) k =1 f ( t k ) − (cid:90) [0 , d f ( x ) dx (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ TV [0 , d ( f ) D ∗ j ( { t , ..., t j } ) . By re-scaling, if f has bounded variation TV [ − r,r ] d ( f ) and s k = 2 rt k − ( r, r, ..., r ) T then we obtain (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (2 r ) d j j (cid:88) k =1 f ( s k ) − (cid:90) [ − r,r ] d f ( x ) dx (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (2 r ) d · TV [ − r,r ] d ( f ) D ∗ j ( { t , ..., t j } ) . Finally, in order to deal with our choice of basis, we need the following.
Lemma 7.6.
Consider the tensor product ψ m ( x ) := ψ m ( x ) · ... · ψ m ( x d ) in d dimensions and let r > .Then (7.3) TV [ − r,r ] d ( ψ m ) ≤ (cid:16) r (cid:112) | m | + 1) (cid:17) d − . OUNDATIONS OF SPECTRAL COMPUTATIONS 29
Proof.
We will use an alternative form of the total variation which holds for smooth enough functions andcan be found in [85]: TV [ − r,r ] d ( ψ m ) = d (cid:88) k =1 (cid:88) ≤ i <...
Given T ∈ Ω or T ∈ Ω and (cid:15) > , we can approximate the matrix values (cid:104) ( T − zI ) ψ m , ( T − zI ) ψ n (cid:105) and (cid:104) ( T − zI ) ∗ ψ m , ( T − zI ) ∗ ψ n (cid:105) to within (cid:15) using finitely many arithmetical operations and comparisons of the relevant information (capturedby Ξ j and Ξ j in § Let T ∈ Ω or T ∈ Ω and (cid:15) > . Recall that T = (cid:88) | k |≤ N a k ( x ) ∂ k , T ∗ = (cid:88) | k |≤ N ˜ a k ( x ) ∂ k , so by expanding out the inner products and also considering the case a k = 1 , it is sufficient to approximate (cid:104) a k ∂ k ψ m , a j ∂ j ψ n (cid:105) and (cid:104) ˜ a k ∂ k ψ m , ˜ a j ∂ j ψ n (cid:105) for all relevant k, j, m and n . Due to the symmetry in the assumptions of T and T ∗ , we only need to showthat one can compute the first inner product, the proof for the second one is identical. Note that by the specificchoice of the basis functions ψ m it follows that ∂ k ψ m can be written as a finite linear combination of tensorproducts of Hermite functions using the recurrence relations (the coefficients in the linear combinations arethus recursively defined as a function of k ). Hence, in the inner product, we can assume that there are nopartial derivatives. In doing this, we have assumed that we can compute square roots of integers (which occurin the coefficients) to arbitrary precision (recall we want an arithmetic tower) which can be achieved by asimple interval bisection routine. It follows that we only need to consider approximations of inner productsof the form (cid:104) a k ψ m , a j ψ n (cid:105) . To do so, let
R > then, by H¨older’s inequality and the assumption of polynomially bounded growth onthe coefficients a k , we have (cid:90) | x i |≥ R | a k a j | | ψ m ψ n | dx ≤ A k A j (cid:32)(cid:90) | x i |≥ R (cid:16) | x | B k (cid:17) (cid:16) | x | B j (cid:17) ψ m ( x ) dx (cid:33) / (cid:32)(cid:90) | x i |≥ R ψ n ( x ) dx (cid:33) / . The first integral on the right hand side can be bounded by (cid:90) R d | x | B ψ m ( x ) dx ≤ (cid:90) R d (cid:0) x + ... + x d (cid:1) B ψ m ( x ) dx, for B = 4( B k + B j ) , since we restrict to | x i | ≥ R with R > and | x | ≤ (cid:107) x (cid:107) . B is even so we canexpand out the product ( x + ... + x d ) B/ ψ m using the recurrence relations for the Hermite functions. Inone dimension this gives xψ m ( x ) = (cid:114) m ψ m − ( x ) + (cid:114) m + 12 ψ m +1 ( x ) ,x ψ m ( x ) = (cid:114) m xψ m − ( x ) + (cid:114) m + 12 xψ m +1 ( x ) , = (cid:114) m (cid:32)(cid:114) m − ψ m − ( x ) + (cid:114) m ψ m ( x ) (cid:33) + (cid:114) m + 12 (cid:32)(cid:114) m + 12 ψ m ( x ) + (cid:114) m + 22 ψ m +2 ( x ) (cid:33) , and so on. We can do the same for tensor products of Hermite functions. In particular, multiplying a tensorproduct of Hermite functions, ψ m , by ( x + ... + x d ) induces a linear combination of at most d such tensorproducts, each with a coefficient of magnitude at most ( | m | + 2) and index with l ∞ norm bounded by | m | + 2 (allowing repetitions). It follows that ( x + ... + x d ) B/ ψ m can be written as a linear combinationof at most (4 d ) B/ such tensor products, each with a coefficient of magnitude at most ( | m | + B ) B . Squaringthis and integrating, the orthogonality and normalisation of the tensor product of Hermite functions impliesthat (cid:90) R d ( x + ... + x d ) B ψ m ( x ) dx ≤ d ) B/ ( | m | + B ) B =: p ( | m | ) . For the other integral, define p ( | n | ) := 4 d ( | n | + 2) . We then have (cid:90) | x i |≥ R ψ n dx ≤ R (cid:90) R d | x | ψ n dx ≤ p ( | n | ) R , by using the same argument as above but with B = 2 .So given δ > and n, m, B, A k , A j , (and d ) we can choose r ∈ N large such that (cid:90) | x i |≥ r | a k a j | | ψ m ψ n | dx ≤ A k A j p ( | m | ) / p ( | n | ) / r ≤ δ. We now have to consider the cases T ∈ Ω or T ∈ Ω separately, noting that it is sufficient to approxi-mate the integral (cid:82) | x i |≤ r a k a j ψ m ψ n dx to any given precision. For notational convenience, let L r ( m ) = (cid:20) σ (cid:18)(cid:16) r (cid:112) | m | + 1) (cid:17) d − (cid:19)(cid:21) so that with σ = 3 d + 1 as in the definition of (cid:107)·(cid:107) A r , we have via Lemma 7.6 that (cid:107) ψ m (cid:107) A r ≤ L r ( m ) . Case 1: T ∈ Ω . Given k, j, m, n, δ and r ∈ N as above, choose M large such that(7.4) (2 r ) d · C ( d ) (cid:0) log( M ) + 1 (cid:1) d M · c r · L r ( m ) · L r ( n ) ≤ δ/ , where C ( d ) is as (7.2) and c r controls the total variation as in (3.6). Again, note that such an M can be chosenin finitely many arithmetic operations and comparisons with the given data and assuming that logarithms andsquare roots can be computed to arbitrary precision (say by a power series representation and bound on theremainder). Using the fact that A r is a Banach algebra (in particular we can bound the norms of product offunctions by the product of their norms) and Theorem 7.5, it follows that (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (2 r ) d M M (cid:88) l =1 a k ( s l ) a j ( s l ) ψ m ( s l ) ψ n ( s l ) − (cid:90) | x i |≤ r a k a j ψ m ψ n dx (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ δ/ , where s l = 2 rt l − ( r, r, ..., r ) T are the rescaled Halton points. Hence it is enough to show that each product a k ( s l ) a j ( s l ) ψ m ( s l ) ψ n ( s l ) can be computed to a given accuracy using finitely many arithmetic operationsand comparisons. Since each s l ∈ Q d we can evaluate a k ( s l ) a j ( s l ) . Note that we can compute exp( − x / to arbitrary precision with finitely many arithmetic operations and comparisons (again say by a power seriesrepresentation and bound on the remainder) and that we can compute the coefficients of the polynomials Q m OUNDATIONS OF SPECTRAL COMPUTATIONS 31 with ψ m ( x ) = Q m ( x ) exp( − x / , using the recursion formulae to any given precision, it follows that wecan compute ψ m ( s l ) ψ n ( s l ) to a given accuracy using finitely many arithmetic operations and comparisons.Using the bounds on the a k and a j and Cram´er’s inequality, we can bound the error in the product and hencethe result follows. Case 2: T ∈ Ω . On the compact cube | x i | ≤ r the double series a k ( x ) a j ( x ) = (cid:88) t ∈ ( Z ≥ ) d (cid:88) s ∈ ( Z ≥ ) d a tk a sj x t + s converges uniformly (recall that { a tk } t ∈ ( Z ≥ ) d are the power series coefficients for a k ) so we can exchangethe series and integration to write(7.5) (cid:90) | x i |≤ r a k a j ψ m ψ n dx = (cid:88) t,s ∈ ( Z ≥ ) d a tk a sj (cid:90) | x i |≤ r x s + t ψ m ( x ) ψ n ( x ) dx. But (cid:12)(cid:12)(cid:12)(cid:82) | x i |≤ r x s + t ψ m ( x ) ψ n ( x ) dx (cid:12)(cid:12)(cid:12) is bounded by r | t | + | s | (cid:82) x ∈ R d | ψ m | | ψ n | dx ≤ r | t | + | s | , by H¨older’s inequal-ity. Let τ = r/ ( r + 1) , then using the fact that we know d r in (3.8), we can bound the tail of the series in(7.5) by d r (cid:88) | t | , | s | >M τ | t | + | s | ≤ d r (cid:88) | t | >M τ | t | d + ... + | td | d , using the fact that | x | ≤ ( | x | + ... + | x d | ) /d . We can explicitly sum this series (as the difference of geometricseries) to gain the bound d r (cid:34) − (1 − τ ( M +1) /d ) d (cid:0) − τ /d (cid:1) d (cid:35) . Given r and d r (and d ) we can keep increasing M and evaluating the bound (strictly speaking an upper boundaccurate to /M say), to choose M large such that the tail is smaller than δ/ for any given δ > . It followsthat it is enough to estimate integrals of the form (cid:82) | x i |≤ r x s + t ψ m ( x ) ψ n ( x ) dx. Using the recurrence relationsfor Hermite functions and writing ψ m ( x ) = Q m ( x ) exp( − x / , it is enough to split the multidimensionalintegral up as products and sums of one-dimensional integrals of the form (cid:82) r − r x a exp( − x ) dx, for a ∈ Z ≥ .Again, we have assumed that we can compute the coefficients of the Q m to any given accuracy using finitelymany arithmetic operations and comparisons and using this we can bound the total error of the expression by δ/ . The above integral vanishes unless a is even, so integration by parts (again assuming we can evaluate exp( − x ) to any desired accuracy) reduces this to estimating (cid:82) r − r exp( − x ) dx. Consider the Taylor seriesfor exp( − x ) . The tail can be bounded by (cid:88) k>N r k k ! ≤ r N N ! exp( − r ) . Integrating this estimate over the interval [ − r, r ] , we can bound this by any given η > by choosing N large enough. We can then explicitly compute (cid:82) r − r (cid:80) k ≤ N x k /k ! dx . Keeping track of all the errors iselementary and hence (cid:82) | x i |≤ r a k a j ψ m ψ n dx can be approximated with finitely many arithmetic operationsand comparisons as required. (cid:3) In some cases, we can also directly compute matrix elements without the cut-off argument used in theabove proof. For instance, if each a k ( x ) (and hence ˜ a k ( x ) ) is a polynomial then we can simply integrate thepower series to compute (cid:104) a k ( x ) ψ m , a j ( x ) ψ n (cid:105) and use the recurrence relations for Hermite functions. If weknow a bound on the degree of the polynomials, then clearly we can compute(7.6) (cid:104) ( T − zI ) ψ m , ( T − zI ) ψ n (cid:105) and (cid:104) ( T − zI ) ∗ ψ m , ( T − zI ) ∗ ψ n (cid:105) to within (cid:15) using finitely many arithmetical operations and comparisons directly. Even if we do not know thedegree of the polynomials and are only promised that each a k ( x ) is a polynomial, then we can successively approximate by more terms of the power series and eventually compute (7.6) to within (cid:15) using finitely manyarithmetical operations and comparisons. Though we do not know when the given accuracy has been reached(recall that we only know a finite portion of the coefficients c , c , ... at any one time for T ∈ Ω ).We can now prove the positive parts of Theorems 3.3 and 3.5. Proof of inclusions in Theorems 3.3 and 3.5.
Step 1 : { Ξ , Ω } , { Ξ , Ω } ∈ Σ A . The proof of thissimply strings together the above results. The linear span of { e , e , ... } (the reordered Hermite func-tions) is a core of T and T ∗ by Proposition 7.1. By Proposition 7.7, we can compute the inner products (cid:104) ( T − zI ) e j , ( T − zI ) e i (cid:105) and (cid:104) ( T − zI ) ∗ e j , ( T − zI ) ∗ e i (cid:105) up to arbitrary precision with finitely many arith-metic operations and comparisons. Using Lemma 7.2, given z ∈ C , we can compute some approximation υ n ( z, T ) in finitely many arithmetic operations and comparisons such that (cid:12)(cid:12) υ n ( z, T ) − min { σ (( T − zI ) | P n ( l ( N )) ) , σ (( T ∗ − ¯ zI ) | P n ( l ( N )) ) } (cid:12)(cid:12) ≤ n . We now set(7.7) γ n ( z, T ) = υ n ( z, T ) + 1 /n. Then γ n satisfies the hypotheses of Proposition 6.5. The proof of Theorem 3.8 also makes clear that we haveerror control since γ n ( z, T ) ≥ (cid:107) R ( z, T ) (cid:107) − . Step 2 : { Ξ , Ω } , { Ξ , Ω } ∈ Σ A . Consider the sequence of functions γ n defined by equation (7.7).These converge uniformly to (cid:107) R ( z, T ) (cid:107) − on compact subsets of C and satisfy γ n ( z, T ) ≥ (cid:107) R ( z, T ) (cid:107) − .We can now apply Proposition 6.6. Step 3: { Ξ , Ω } , { Ξ , Ω } ∈ ∆ A . Let T ∈ Ω . Our strategy will be to compute the inner products (cid:104) ( T − zI ) e j , ( T − zI ) e i (cid:105) and (cid:104) ( T − zI ) ∗ e j , ( T − zI ) ∗ e i (cid:105) to an error which decays rapidly enough aswe let the cut-off parameter r tend to ∞ . We follow the proof of Proposition 7.7 closely. Recall that given n, m , we can choose r ∈ N large such that (cid:90) | x i |≥ r | a k a j | | ψ m ψ n | dx ≤ A k A j p ( | m | ) / p ( | n | ) / r , with the crucial difference that now we do not assume we can compute A k , A j , p or p . It follows that thereexists some polynomial p , with coefficients not necessarily computable from the given information, suchthat (cid:90) | x i |≥ r | a k a j | | ψ m ψ n | dx ≤ p ( | m | , | n | ) r , for all | j | , | k | ≤ N . Now we use the sequence b r to bound the error in the integral over the compactcube asymptotically. We assume without loss of generality that b r is increasing monotonically to ∞ with r . Using Halton sequences and the same argument in the proof of Proposition 7.7, we can approximate (cid:82) | x i |≤ r a k a j ψ m ψ n dx, with an error that, asymptotically up to some unknown constant, is bounded by(7.8) r d · (cid:0) log( M ) + 1 (cid:1) d M · b r · L r ( m ) · L r ( n ) , where M is the number of Halton points. We can let M depend on r, n and m such that (7.8) is boundedby a constant times /r . It follows that we can bound the total error in approximating (cid:104) a k ψ m , a j ψ n (cid:105) forany j, k by p ( | m | , | n | ) /r , by making the coefficients of p larger if necessary. We argue similarly for theadjoint and note that (cid:104) ( T − zI ) ψ m , ( T − zI ) ψ n (cid:105) and (cid:104) ( T − zI ) ∗ ψ m , ( T − zI ) ∗ ψ n are both approximatedto within (1 + | z | ) P ( | m | , | n | ) r , for some unknown polynomial P . Hence we can apply Lemma 7.2 (the form where we do not know the errorin inner product estimates), changing the polynomial P to take into account the basis mapping from Z d ≥ to OUNDATIONS OF SPECTRAL COMPUTATIONS 33 N to some polynomial Q , to gain some approximation υ n ( z, T ) in finitely many arithmetic operations andcomparisons such that(7.9) (cid:12)(cid:12) υ n ( z, T ) − min { σ (( T − zI ) | P n ( l ( N )) ) , σ (( T ∗ − ¯ zI ) | P n ( l ( N )) ) } (cid:12)(cid:12) ≤ n (1 + | z | ) Q ( n ) r ( n, z ) + 1 n . We now choose r ( z, n ) larger if necessary such that r ( z, n ) ≥ (1 + | z | ) exp( n ) . We now set γ n ( z, T ) = υ n ( z, T ) + 1 /n. Then γ n satisfies the hypotheses of Proposition 6.5 and Proposition 6.6 since the error in(7.9) decays faster than /n . We can use these propositions to build the required arithmetical algorithm. Step 4: { Ξ , Ω } , { Ξ , Ω } ∈ ∆ A . We argue as in step 3. To control the error in the approximationof the integral over a compact hypercube, choose the cut-off M ( r ) such that (cid:88) | t | , | s | >M ( r ) (cid:18) rr + 1 (cid:19) | t | + | s | ≤ b r r . It follows that there exists some (unknown) constant B such that we can bound the error in approximating (cid:82) | x i |≤ r a k a j ψ m ψ n dx by B/r where we have absorbed the arbitrarily small error that comes from approxi-mating the integral of the truncated power series using finitely many arithmetic operations and comparisons.The rest of the argument is the same as in step 3. (cid:3) Proofs of impossibility results in Theorems 3.3 and 3.5.
We first deal with Theorem 3.3. Recall themaps Ξ j , Ξ j : Ω , Ω (cid:51) T (cid:55)→ Sp( T ) ∈ M AW j = 1Sp (cid:15) ( T ) ∈ M AW j = 2 , We split up the arguments to deal with Ω and then Ω . Proof that { Ξ j , Ω } / ∈ ∆ G . Suppose first for a contradiction that a height one tower, Γ n , exists for theproblem { Ξ , Ω } such that d AW (Γ n ( T ) , Ξ ( T )) ≤ − n . We will deal with the one-dimensional case andhigher dimensions are similar. Let ρ ( x ) be any smooth bump function with maximum value , minimumvalue and support [0 , . Let ρ n denote the translation of ρ to have support [ n, n + 1] . We will consider thetwo (self-adjoint and bounded) operators ( T u )( x ) = 0 , ( T m u )( x ) = ρ m ( x ) u ( x ) , which have spectra { } and [0 , respectively. For these we can take the polynomial bound (the { A k } and { B k } ) to be and the total variation bound to be c r = 1 + σ TV [0 , ( ρ ) . When we compute Γ ( T ) , we onlyuse finitely many evaluations of the coefficient function a ( x ) = 0 (as well as the other given information).We can then choose m large such that the support of ρ m does not intersect the points of evaluation. Byassumptions (ii) and (iii) in Definition 5.1, Γ ( T m ) = Γ ( T ) . But this contradicts the triangle inequalitysince d AW ( { } , [0 , ≥ To argue for the pseudospectrum let (cid:15) > and note that (cid:15) / ∈ Sp (cid:15) ( T ) but (cid:15) ∈ Sp (cid:15) ( (cid:15)T m ) . We now alterthe given c r to (cid:15) (1 + σ TV [0 , ( ρ )) and the polynomial bound to (cid:15) . The argument is now exactly as before.Namely, we choose n large such that d AW (Γ n ( T ) , [ − (cid:15), (cid:15) ]) > − n then choose m large such that Γ n ( T ) = Γ n ( (cid:15)T m ) . (cid:3) Proof that { Ξ j , Ω } / ∈ Σ G ∪ Π G . Suppose first of all that a Σ G tower, Γ n , exists for { Ξ , Ω } . We willdeal with the one-dimensional case and higher dimensions are similar. Consider the operators ( T u )( x ) = 0 , ( T u )( x ) = f ( x ) u ( x ) , where we define f in terms of Γ n as follows. We will ensure that f ( x ) = 1 except for finitely many valuesof x where it takes the value and hence T and T have spectra { } and { } respectively and are both self-adjoint. Note that once the zeros of f are fixed, this choice ensures that f has total variation boundedby a constant on any hypercube and hence we may take b r = 1 for all r ∈ N . There exists some n such that Γ n ( T ) contains z n ∈ B / (0) with a guaranteed error estimate of dist( z n , Sp( T )) ≤ / . But Γ n ( T ) can only depend on finitely many evaluations of (as well as b r = 1 and the trivial choice of g j ( x ) = x ).We choose f to be zero at precisely these evaluation points. By assumptions (ii) and (iii) in Definition 5.1, Γ n ( T ) = Γ n ( T ) , including the given error estimates, which is the required contradiction.For { Ξ , Ω } / ∈ Σ G , given (cid:15) > we replace f by (cid:15)f in the above argument and keep all otherinputs the same. Hence T and T have (cid:15) -pseudospectra [ − (cid:15), (cid:15) ] and [2 (cid:15), (cid:15) ] respectively. We note thatagain there exists some n such that Γ n ( T ) contains z n ∈ B (cid:15)/ (0) with a guaranteed error estimate of dist( z n , Sp (cid:15) ( T )) ≤ (cid:15)/ . But Γ n ( T ) can only depend on finitely many evaluations of (as well as b r = 1 and the trivial choice of g j ( x ) = x ). We choose f to be zero at precisely these evaluation points. Byassumptions (ii) and (iii) in Definition 5.1, Γ n ( T ) = Γ n ( T ) , including the given error estimates, which isthe required contradiction.To argue that neither problem lies in Π G , we can use the same arguments in the proof that { Ξ j , Ω } / ∈ ∆ G . The only change now is that the algorithm, Γ n , used to derive the contradiction provides Π G informationrather than ∆ G . For the spectrum, we consider the operators ( T u )( x ) = 0 and ( T m u )( x ) = ρ m ( x ) u ( x ) , and choose n large such that Γ n ( T ) produces the guarantee Sp( T ) ∩ B / (0) c = ∅ . For m sufficientlylarge, we argue as before to get Γ n ( T m ) = Γ n ( T ) , including the guarantee, the required contradiction.Again a similar argument works for the pseudospectrum by rescaling T m to (cid:15)T m . (cid:3) We now deal with the impossibility results in Theorem 3.5 where Ξ j , Ξ j : Ω , Ω (cid:55)→ Sp( T ) ∈ M AW j = 1Sp (cid:15) ( T ) ∈ M AW j = 2 . Proof that { Ξ j , Ω } / ∈ ∆ G . Suppose for a contradiction that a height one tower, Γ n , exists for { Ξ , Ω } such that d AW (Γ n ( T ) , Ξ ( A )) ≤ − n . Now consider the two (self-adjoint and bounded) operators ( T u )( x ) = 0 and ( T u )( x ) = x k exp( − x ) u ( x ) /s k , where k is even and will be chosen later. We choose s k such that the range of the function is x k exp( − x ) /s k is [0 , and hence T has spectrum [0 , . We can take the polynomial bounding function to be the constant for both operators and must show that we can use the same d r for both operators in (3.8), independent of k . Simple calculus yields that s k = ( k/ (2 e )) k/ . It follows that such a d r must satisfy(7.10) (cid:18) ek (cid:19) k/ ( r + 1) m + k m ! ≤ d r , ∀ k ∈ N , m ∈ Z ≥ . Hence it suffices to show that the function on the left hand side of (7.10) is bounded (as a function of m, k for all r ∈ N ). Using Stirling’s approximation (explicitly the bounds on m ! ) this will follow if we can show r m + k k k/ m m +1 / ≤ (cid:18) r √ k (cid:19) k (cid:18) r √ m (cid:19) m is bounded for all r ∈ N (now with m > . But this is obvious.We can now choose k (which depends on the algorithm Γ n ) to gain a contradiction. Since Sp( T ) = { } and ∈ Sp( T ) for all even k , there exists n such that dist(1 , Γ n ( T )) > / but dist(1 , Γ n ( T )) < / .However, Γ n ( T ) can only depend on finitely many of the coefficients { c j } , say c , ..., c ˜ N ( T,n ) , of T (aswell as the other given information). By assumption (iii) in Definition 5.1, we can choose k such that thecoefficient corresponding to x k , call it c l k , has l k > ˜ N ( T , n ) and get Γ n ( T ) = Γ n ( T ) , the requiredcontradiction. OUNDATIONS OF SPECTRAL COMPUTATIONS 35
To show { Ξ , Ω } / ∈ ∆ G uses exactly the same argument as above. In order to gain the neces-sary separation (cid:15) / ∈ Sp (cid:15) ( T ) but (cid:15) ∈ Sp (cid:15) ( T ) , we rescale T to (cid:15)T . Then there exists n such that dist(3 (cid:15), Γ n ( T )) > (cid:15)/ but dist(3 (cid:15), Γ n ( T )) < (cid:15)/ . The rest of the contradiction follows. (cid:3) Proof that { Ξ j , Ω } , { Ξ j , Ω p } / ∈ Σ G ∪ Π G . Since Ω p ⊂ Ω , it is enough to show the results for Ω p .Suppose for a contradiction that there exists a Σ G algorithm, Γ n , for { Ξ , Ω p } . Consider ( T u )( x ) = xu ( x ) and ( T u )( x ) = ( x − x k ) u ( x ) , where k is even and chosen later. ( T j ± iI ) C ∞ ( R ) are dense in L ( R ) with T j initially defined on C ∞ ( R ) symmetric. It follows that the closure of T j | C ∞ ( R ) is self-adjoint and hence that T j ∈ Ω p . Note that Sp( T ) = R but Sp( T ) ⊂ ( −∞ , . Now choose n such that Γ n ( T ) contains a point z n ∈ B / (2) with aguaranteed error estimate of dist( z n , Sp( T )) ≤ / . However, Γ n ( T ) can only depend on the first ˜ N ( T, n ) coefficients, c , ..., c ˜ N ( T,n ) , of T (as well as the trivial choice g j ( x ) = x and the numbers b n = n ! ). Byassumption (iii) in Definition 5.1, we can choose k such that the coefficient corresponding to x k , call it c r k ,has r k > ˜ N ( T , n ) and get Γ n ( T ) = Γ n ( T ) , the required contradiction. Similarly by rescaling as above,we get { Ξ , Ω p } / ∈ Σ G .To show { Ξ , Ω p } / ∈ Π G we argue the same way, but now set ( T u )( x ) = 0 and ( T u )( x ) = x k u ( x ) . Asbefore, T j ∈ Ω p , but now Sp( T ) = { } and ∈ Sp( T ) . Choose n such that Γ n ( T ) produces the guarantee Sp( T ) ∩ B / (0) c = ∅ . Again, choose k such that c r k has r k > ˜ N ( T , n ) and get Γ n ( T ) = Γ n ( T ) , therequired contradiction. Rescaling and using the same argument shows { Ξ , Ω p } / ∈ Π G . (cid:3)
8. P
ROOFS OF T HEOREMS ON D ISCRETE S PECTRA
Here we prove our results related to the discrete spectrum. We need some results on finite section approxi-mations to the discrete spectrum of a Hermitian operator below the essential spectrum. There are two cases toconsider; either there are infinitely many eigenvalues below the essential spectrum, or there are only finitelymany. The following are well-known and follow from the ‘min-max’ theorem characterising eigenvalues.
Lemma 8.1.
Let B ∈ B ( l ( N )) be self-adjoint with eigenvalues λ ≤ λ ≤ ... (infinitely many, countedaccording to multiplicity) below the essential spectrum. Consider the finite section approximates B n = P n BP n ∈ C n and list the eigenvalues of B n as µ n ≤ µ n ≤ ... ≤ µ nn . Then the following hold:(1) λ j ≤ µ nj for j = 1 , ..., n ,(2) for any j ∈ N , µ nj ↓ λ j as n → ∞ ( n ≥ j so that µ nj makes sense). Lemma 8.2.
Let B ∈ B ( l ( N )) be self-adjoint with finitely many eigenvalues λ ≤ λ ≤ ... ≤ λ m (countedaccording to multiplicity) below the essential spectrum and let a = inf { x : x ∈ Sp ess ( B ) } . For j > m weset λ j = a . Consider the finite section approximates B n = P n BP n ∈ C n and list the eigenvalues of B n as µ n ≤ µ n ≤ ... ≤ µ nn . Then the following hold:(1) λ j ≤ µ nj for j = 1 , ..., n ,(2) for any j ≤ m , µ nj ↓ λ j as n → ∞ ( n ≥ j so that µ nj makes sense),(3) given (cid:15) > and k ∈ N , there exists N such that for all n ≥ N , µ nk ≤ a + (cid:15) .Proof of Theorem 3.13 for Ξ d . Step 1 : { Ξ d , Ω d D } / ∈ ∆ G . Suppose this were false and that there exists someheight one tower Γ n solving the problem. Consider the matrix operators A m = diag { , , ..., , } ∈ C m × m and C = diag { , , ... } and set A = diag { , } ⊕ ∞ (cid:77) m =1 A k m , where we choose an increasing sequence k m inductively as follows. Set k = 1 and suppose that k , ..., k m have been chosen. Sp d (diag { , } ⊕ A k ⊕ A k ⊕ ... ⊕ A k m ⊕ C ) = { , } is closed and so there exists some n m ≥ m such that if n ≥ n m then(8.1) dist(2 , Γ n (diag { , } ⊕ A k ⊕ ... ⊕ A k m ⊕ C ) ≤ . Now let k m +1 ≥ max { N (diag { , } ⊕ A k ⊕ ... ⊕ A k m ⊕ C, n m ) , k m + 1 } . Arguing as in the proof ofTheorem 3.9, it follows that Γ n m ( A ) = Γ n m (diag { , } ⊕ A k ⊕ ... ⊕ A k m ⊕ C ) . But Γ n m ( A ) convergesto Sp d ( A ) = { } , contradicting (8.1). Step 2 : { Ξ d , Ω d N } ∈ Σ A . We now construct an arithmetic height two tower for Ξ d and the class Ω d N . To dothis, we recall that a height two tower ˜Γ n ,n for the essential spectrum of operators in Ω d N was constructedin [8]. For completeness, we write out the algorithm here. Let P n be the usual projection onto the first n basis elements and set Q n = I − P n . Define µ m,n ( A ) := min { σ ( P f ( n ) ( A − zI ) | Q m P n ( l ( N )) ) , σ ( P f ( n ) ( A − zI ) ∗ | Q m P n ( l ( N )) ) } ,G n := min (cid:26) s + it n : s, t ∈ {− n , ..., n } (cid:27) , Υ m ( z ) := z + { w ∈ C : | Re( w ) | , | Im( w ) | ≤ − ( m +1) } . We then define the following sets for n > m : S m,n ( z ) := { j = m + 1 , ..., n : ∃ w ∈ Υ m ( z ) ∩ G j with µ m,i ( w ) ≤ /m } ,T m,n ( z ) := { j = m + 1 , ..., n : ∃ w ∈ Υ m ( z ) ∩ G j with µ m,i ( w ) ≤ / ( m + 1) } ,E m,n ( z ) := | S m,n ( z ) | + | T m,n ( z ) | − n,I m,n := (cid:26) z ∈ (cid:26) s + it m : s, t ∈ Z (cid:27) : E m,n ( z ) > (cid:27) . Finally we define for n > n ˜Γ n ,n ( A ) = (cid:91) z ∈ I n ,n Υ n ( z ) , and set ˜Γ n ,n ( A ) = { } if n ≤ n . Furthermore, the tower has the following desirable properties:(1) For fixed n , the sequence ˜Γ n ,n ( A ) is eventually constant as we increase n ,(2) The sets lim n →∞ ˜Γ n ,n ( A ) =: ˜Γ n ( A ) are nested, converging down to Sp ess ( A ) .We also need the height one tower, ˆΓ n , for the spectrum of operators in Ω d N discussed in § §
6. Notethat ˆΓ n ( A ) is a finite set for all n . For z ∈ ˆΓ n ( z ) , this also outputs an error control E ( n, z ) such that dist( z, Sp( A )) ≤ E ( n, z ) and such that E ( n, z ) converges to the true distance to the spectrum uniformlyon compact subsets of C (with the choice of g ( x ) = x since the operator is normal). We now fit the piecestogether and initially define ζ n ,n ( A ) = { z ∈ ˆΓ n ( A ) : E ( n , z ) < dist( z, ˜Γ n ,n ( A ) + B /n (0)) } . We must show that this defines an arithmetic tower in the sense of Definitions 5.1 and 5.3. Given z ∈ ˆΓ n ( A ) and using Pythagoras’ theorem, along with the fact that ˜Γ n ,n ( A ) consists of finitely many squaresin the complex plane aligned with the real and imaginary axes, we can compute dist( z, ˜Γ n ,n ( A )) infinitely many arithmetic operations and comparisons. We can compute ( E ( n , z ) + 1 /n ) and check if thisis less than dist( z, ˜Γ n ,n ( A )) . Hence ζ n ,n ( A ) can be computed with finitely many arithmetic operationsand comparisons. There are now two cases to consider: Case 1: Sp d ( A ) ∩ (˜Γ n ( A ) + B /n (0)) c = ∅ . For large n , ˜Γ n ( A ) = ˜Γ n ,n ( A ) and this set containsthe essential spectrum. It follows, for large n , since E ( n , z ) ≥ dist( z, ˜Γ n ,n ( A )) for all z ∈ ˆΓ n ( A ) , that ζ n ,n ( A ) = ∅ . Case 2: Sp d ( A ) ∩ (˜Γ n ( A ) + B /n (0)) c (cid:54) = ∅ . In this case, this set is a finite subset of Sp d ( A ) , { ˆ z , ..., ˆ z m ( n ) } , separated from the closed set ˜Γ n ( A ) + B /n (0) (we need the + B /n (0) for this to The actual algorithm is slightly more complicated to avoid the empty set, but its listed properties still hold.
OUNDATIONS OF SPECTRAL COMPUTATIONS 37 be true to avoid accumulation points of the discrete spectrum). There exists some δ n > such that theballs B δ n (ˆ z j ) for j = 1 , ..., m ( n ) are pairwise disjoint and such that their union does note intersect ˜Γ n ( A ) + B /n (0) .Using the convergence of ˆΓ n ( A ) to Sp( A ) and E ( n, z ) ≥ dist( z, Sp( A )) , it followsthat for large n that(8.2) ζ n ,n ( A ) ⊂ m ( n ) (cid:91) j =1 B δ n (ˆ z j ) , is non-empty and that ζ n ,n ( A ) converges to Sp d ( A ) ∩ (˜Γ n ( A ) + B /n (0)) c (cid:54) = ∅ in the Hausdorff metric.Suppose that ζ n ,n ( A ) is non-empty. Recall that we only want one output per eigenvalue in the discretespectrum. To do this, we partition the finite set ζ n ,n ( A ) into equivalence classes as follows. For z, w ∈ ζ n ,n ( A ) , we say that z ∼ n w if there exists a finite sequence z = z , z , ..., z n = w ∈ ζ n ,n ( A ) such that B E ( n ,z j ) ( z j ) and B E ( n ,z j +1 ) ( z j +1 ) intersect. The idea is that equivalence classes correspond to clustersof points in ζ n ,n ( A ) . Given any z ∈ ζ n ,n ( A ) we can compute its equivalence class using finitely manyarithmetic operations and comparisons. Let S be the set { z } and given S n , let S n +1 be the union of any w ∈ ζ n ,n ( A ) such that B E ( n ,w ) ( w ) and B E ( n ,v ) ( v ) intersect for some v ∈ S n . Given S n , we cancompute S n +1 using finitely many arithmetic operations and comparisons. The equivalence class is any S n where S n = S n +1 which must happen since ζ n ,n ( A ) is finite. We let Φ n ,n consist of one element of eachequivalence class that minimises E ( n , · ) over its respective equivalence class. By the above comments itis clear that Φ n ,n can be computed in finitely many arithmetic operations and comparisons from the givendata. Furthermore, due to (8.2) which holds for large n , the separation of the B δ n (ˆ z j ) and the fact that E ( n , · ) converges uniformly on compact subsets to the distance to Sp( A ) , it follows that for large n thereis exactly one point in each intersection B δ n (ˆ z j ) ∩ Φ n ,n ( A ) . But we can shrink δ n and apply the sameargument to see that Φ n ,n ( A ) converges to Sp d ( A ) ∩ (˜Γ n ( A ) + B /n (0)) c (cid:54) = ∅ in the Hausdorff metric.Now suppose that ζ n ,n ( A ) is non-empty and z , z ∈ Φ n ,n ( A ) and both lie in B (cid:15) ( z ) for some z ∈ Sp d ( A ) and (cid:15) > with Sp( A ) ∩ B (cid:15) ( z ) = { z } . It follows that z minimises the distance to the spectrumfrom both z and z . Hence, B E ( n ,z ) ( z ) and B E ( n ,z ) ( z ) both contain the point z so that z ∼ n z .But then at most one of z , z can lie in Φ n ,n ( A ) and hence z = z .To finish, we must alter Φ n ,n ( A ) to take care of the case when ζ n ,n ( A ) = ∅ and to produce a Σ A algorithm. In the case that ζ n ,n ( A ) = ∅ , set Φ n ,n ( A ) = ∅ . Let N ( A ) ∈ N be minimal such that Sp d ( A ) ∩ (˜Γ N ( A ) + B /N (0)) c (cid:54) = ∅ (recall the discrete spectrum is non-empty for our class of operators). If n > n then set Γ n ,n ( A ) = { } , otherwise consider Φ k,n ( A ) for n ≤ k ≤ n . If all of these are emptythen set Γ n ,n ( A ) = { } , otherwise choose minimal k with Φ k,n ( A ) (cid:54) = ∅ and let Γ n ,n ( A ) = Φ k,n ( A ) .Note that this defines an arithmetic tower of algorithms, with Γ n ,n ( A ) non-empty. By the above caseanalysis, for large n it holds that Γ n ,n ( A ) = Φ n ∨ N ( A ) ,n ( A ) and it follows that(8.3) lim n →∞ Γ n ,n ( A ) =: Γ n ( A ) = Sp d ( A ) ∩ (˜Γ n ∨ N ( A ) ( A ) + B /n ∨ N ( A ) (0)) c . Hence Γ n ( A ) ⊂ Sp d ( A ) and Γ n ( A ) converges up to Sp d ( A ) in the Hausdorff metric. Step 3 : Multiplicities.
Suppose that z n ,n ∈ Γ n ,n ( A ) converges as n → ∞ to some z n = z ∈ Γ n ( A ) ⊂ Sp d ( A ) , where Γ n is the first limit of the height two tower constructed in step 2. Considerthe following operator, viewed as a finite matrix acting on C n , A n = P n ( A − zI ) ∗ ( A − zI ) P n . This is atruncation of the operator ( A − zI ) ∗ ( A − zI ) . The key observation is that lies in the discrete spectrum of ( A − zI ) ∗ ( A − zI ) with h (( A − zI ) ∗ ( A − zI ) ,
0) = h ( A, z ) , the multiplicity of the eigenvalue z . To see this, note that ker( A − zI ) = ker(( A − zI ) ∗ ( A − zI )) and that if (cid:107) x (cid:107) = 1 then (cid:107) ( A − zI ) x (cid:107) ≤ (cid:112) (cid:107) ( A − zI ) ∗ ( A − zI ) x (cid:107) . Since ( A − zI ) is bounded below on ker( A − zI ) ⊥ , the same must be true for ( A − zI ) ∗ ( A − zI ) . Now set h n ,n ( A, z n ,n ) = min { n , (cid:12)(cid:12) { w ∈ Sp( P n ( A − z n ,n I ) ∗ P f ( n ) ( A − z n ,n I ) P n ) : | w | < /n − d n } (cid:12)(cid:12) } , where d n is some non-negative sequence converging to that we define below. As usual we consider therelevant operator as a matrix acting on C n and we count eigenvalues according to their multiplicity. Viashifting by (1 /n − d n ) I and assuming d n can be computed with finitely many arithmetic operations andcomparisons, Lemma 6.8 shows that this is a general algorithm and can be computed with finitely manyarithmetic operations and comparisons. Consider the similar function (that we cannot necessarily computesince we do not know z ), q n ,n ( A, z ) = min { n , |{ w ∈ Sp( A n ) : | w | < /n }|} , where A n = P n ( A − zI ) ∗ ( A − zI ) P n . We set B = ( A − zI ) ∗ ( A − zI ) and list λ ≤ λ ≤ ... as in Lemmas 8.1 and 8.2, then lim n →∞ q n ,n ( A, z ) = min { n , | λ j : λ j < /n |} . It is then clear from the same lemmas that lim n →∞ lim n →∞ q n ,n ( A, z ) = h (( A − zI ) ∗ ( A − zI ) ,
0) = h ( A, z ) . We will have completed the proof if we can choose d n such that lim n →∞ | h n ,n ( A, z n ,n ) − q n ,n ( A, z ) | = 0 . It is straightforward to show that (cid:13)(cid:13) A n − P n ( A − z n ,n I ) ∗ P f ( n ) ( A − z n ,n I ) P n (cid:13)(cid:13) ≤ (cid:0) | z − z n ,n | + c n (cid:1)(cid:0) (cid:107) AP n (cid:107) + | z − z n ,n | + | z n ,n | + (cid:13)(cid:13) P f ( n ) ( A − z n ,n I ) P n (cid:13)(cid:13) (cid:1) ≤ (cid:0) | z − z n ,n | + c n (cid:1)(cid:0) (cid:13)(cid:13) P f ( n ) AP n (cid:13)(cid:13) + | z − z n ,n | + 2 | z n ,n | + c n (cid:1) , where D f,m ( A ) ≤ c m is the dispersion bound. Choose d n = (cid:0) E ( n , z n ,n ) + c n (cid:1)(cid:0) E ( n , z n ,n ) + 2 | z n ,n | + 2 k n + c n (cid:1) , where k n overestimates (cid:13)(cid:13) P f ( n ) AP n (cid:13)(cid:13) by at most . k n can be computed using a similar positive defi-niteness test as in DistSpec (see Appendix A). Since z n ,n converges to z ∈ Sp d ( A ) , it is clear that (cid:13)(cid:13) A n − P n ( A − z n ,n I ) ∗ P f ( n ) ( A − z n ,n I ) P n (cid:13)(cid:13) ≤ d n eventually and that d n converges to 0. Weyl’s inequality for eigenvalue perturbations of Hermitian matricesimplies the needed convergence. (cid:3) Proof of Theorem 3.13 for Ξ d . Since Ω D ⊂ Ω f N , its suffices to show that { Ξ d , Ω f N } ∈ Σ A and { Ξ d , Ω D } / ∈ ∆ G . Step 1 : { Ξ d , Ω D } / ∈ ∆ G . The proof is almost identical to step 1 in the proof of Theorem 3.13 for Ξ d .Suppose there exists some height one tower Γ n solving the problem. Consider the matrix operators A m =diag { , , ..., , } ∈ C m × m and C = diag { , , ... } and set A = ∞ (cid:77) m =1 A k m , OUNDATIONS OF SPECTRAL COMPUTATIONS 39 where we choose an increasing sequence k m inductively as follows. Set k = 1 and suppose that k , ..., k m have been chosen. Sp d ( A k ⊕ A k ⊕ ... ⊕ A k m ⊕ C ) = { } so there exists some n m ≥ m such that if n ≥ n m then Γ n ( A k ⊕ ... ⊕ A k m ⊕ C ) = 1 . Now let k m +1 ≥ max { N (diag { , } ⊕ A k ⊕ ... ⊕ A k m ⊕ C, n m ) , k m + 1 } . Arguing as in the proof ofTheorem 3.9, it follows that Γ n m ( A ) = Γ n m ( A k ⊕ ... ⊕ A k m ⊕ C ) . But Γ n m ( A ) converges to as A hasno discrete spectrum and this contradiction finishes this step. Step 2 : { Ξ d , Ω f N } ∈ Σ A . Consider the height two tower, ζ n ,n , defined in step 2 of the proof of Theorem3.13 for Ξ d . Let A ∈ Ω f N and if ζ n ,n ( A ) = ∅ , define ρ n ,n ( A ) = 0 , otherwise define ρ n ,n ( A ) = 1 . Thediscussion in the proof of Theorem 3.13 for Ξ d shows that lim n →∞ ρ n ,n ( A ) =: ρ n ( A ) = , if Sp d ( A ) ∩ (˜Γ n ( A ) + B /n (0)) c = ∅ , otherwise.Since Sp d ( A ) ∩ (˜Γ n ( A ) + B /n (0)) c increases to Sp d ( A ) , it follows that lim n →∞ ρ n ( A ) = Ξ d ( A ) andthat if ρ n ( A ) = 1 , then Ξ d ( A ) = 1 . Hence, ρ n ,n provides a Σ A tower for { Ξ d , Ω f N } . (cid:3) Proof of Theorem 3.14.
Let (cid:15) = ( E ( n , z n ) + δ ) and consider the matrix B = P n ( A − z n I ) ∗ P f ( n ) ( A − z n I ) P n − (cid:15)I n ∈ C n × n , where I n is the n × n identity matrix. B is a Hermitian matrix and is not positive semi-definite. It followsthat B can be put into the form P BP T = LDL ∗ , where L is lower triangular with ’s along its diagonal, D is block diagonal with block sizes × or × and P is a permutation matrix. This can be computed in finitely many arithmetic operations and comparisons.Without loss of generality we assume that P = I . Let x be an eigenvector of B with negative eigenvaluethen set y = L ∗ x . Such an x exists by assumption. Note that (cid:104) y, Dy (cid:105) = (cid:104) L ∗ x, DL ∗ x (cid:105) = (cid:104) x, Bx (cid:105) < . It follows that there exists a unit vector y n with (cid:104) y n , Dy n (cid:105) < . Such a vector is easy to spot if a value inone of the × blocks of D is negative. If not then we need to consider × blocks. Using the argumentin the proof of Lemma 6.8 we can find a × block with a negative eigenvalue by computing the trace anddeterminant. Without loss of generality we assume that this block is the upper × portion of D . It followsthat there exist real numbers a, b , not both equal to , such that y n = ( a, b, , ..., T has (cid:104) y n , Dy n (cid:105) < .If D , < then we can take a = 0 , b = 1 . Otherwise, a (cid:54) = 0 so set a = 1 . We then note that there is anopen interval J such that if b ∈ J then y n = ( a, b, , ..., T has (cid:104) y n , Dy n (cid:105) < . We can now perform asearch routine on R with finer and finer spacing to find such a b . L ∗ is invertible and upper triangular so we can efficiently solve for ˜ x n = ( L ∗ ) − y n using finitelymany arithmetic operations and comparisons. We then approximately normalise ˜ x n by computing (cid:107) ˜ x n (cid:107) ≈ t n ( ρ ) > to precision ρ > using arithmetic operations and comparisons. If we set x n = ˜ x n /t n ( ρ ) then − ρt n ( ρ ) = t n ( ρ ) − ρt n ( ρ ) ≤ (cid:107) x n (cid:107) ≤ t n ( ρ ) + ρt n ( ρ ) = 1 + ρt n ( ρ ) . So we successively choose ρ smaller until we reach ρ n such that ρ n /t n ( ρ n ) < δ . This is always possiblesince lim ρ ↓ t n ( ρ ) = (cid:107) ˜ x n (cid:107) > . Let t n = t n ( ρ n ) , then (cid:104) x n , Bx n (cid:105) = t − n (cid:104) L ∗ ˜ x n , DL ∗ ˜ x n (cid:105) = t − n (cid:104) y n , Dy n (cid:105) < . Note that (cid:13)(cid:13) P f ( n ) ( A − z n I ) x n (cid:13)(cid:13) = (cid:104) x n , Bx n (cid:105) + (cid:107) x n (cid:107) (cid:15) < (cid:107) x n (cid:107) (cid:15). Taking square roots and recalling that D f,n ( A ) ≤ c n and the definition of D f,n finishes the proof. (cid:3) Note that even in the finite-dimensional case this type of error control is the best possible owing to nu-merical errors due to round-off and finite precision. This method is quick and can easily be efficientlyimplemented.
Proof of Theorem 3.15.
Step 1: { Ξ d , Ω d } / ∈ ∆ G . Suppose for a contradiction that Γ n ,n is a height twotower solving this problem. For this proof we shall use one of the decision problems in [41] that were provento have SCI G = 3 . Let ( M , d ) be the discrete space { , } , let Ω (cid:48) denote the collection of all infinite matrices { a i,j } i,j ∈ N with entries a i,j ∈ { , } and consider the problem function Ξ (cid:48) ( { a i,j } ) : “Does { a i,j } have only finitely many columns containing only finitely may non-zero entries?”We will gain a contradiction by using the supposed height two tower for { Ξ d , Ω d } , Γ n ,n , to solve { Ξ (cid:48) , Ω (cid:48) } .Without loss of generality, identify B ( l ( N )) with B ( X ) where X = C ⊕ (cid:76) ∞ j =1 X j in the l -sense with X j = l ( N ) . Now let { a i,j } ∈ Ω (cid:48) and for the j th column define B j ∈ B ( X j ) with the following matrixrepresentation: B j = M j (cid:77) r =1 A l jr , A m := . . .
01 1 ∈ C m × m , where if M j is finite then l jM j = ∞ with A ∞ = diag(1 , , , ... ) . The l jr are defined such that(8.4) (cid:80) mi =1 a i,j (cid:88) r =1 l jr = m + m (cid:88) i =1 a i,j . Define the self-adjoint operator A = diag { , } ⊕ ∞ (cid:77) j =1 B j . Note that no matter what the choices of l jr are, ∈ Sp d ( A ) and hence A ∈ Ω d . Note also that the spectrumof A is contained in { , , , } . If Ξ (cid:48) ( { a i,j } ) = 1 then is an isolated eigenvalue of finite multiplicity andhence in Sp d ( A ) . But if Ξ (cid:48) ( { a i,j } ) = 0 then is an isolated eigenvalue of infinite multiplicity so does notlie in the discrete spectrum and hence Sp d ( A ) ⊂ { , , } .Consider the intervals J = [0 , / , and J = [3 / , ∞ ) . Set α n ,n = dist(1 , Γ n ,n ( A )) . Let k ( n , n ) ≤ n be maximal such that α n ,k ( A ) ∈ J ∪ J . If no such k exists or α n ,k ( A ) ∈ J then set ˜Γ n ,n ( { a i,j } ) = 1 . Otherwise set ˜Γ n ,n ( { a i,j } ) = 0 . It is clear from (8.4) that this defines a generalisedalgorithm. In particular, given N we can evaluate { A k,l : k, l ≤ N } using only finitely many evaluations of { a i,j } , where we can use a suitable bijection between bases of l ( N ) and C ⊕ (cid:76) ∞ j =1 X j to view A as actingon l ( N ) . The point of the intervals J , J is that we can show lim n →∞ ˜Γ n ,n ( { a i,j } ) = ˜Γ n ( { a i,j } ) ex-ists. If Ξ (cid:48) ( { a i,j } ) = 1 , then, for large n , lim n →∞ α n ,k ( A ) < / and hence lim n →∞ ˜Γ n ( { a i,j } ) = 1 .Similarly, if Ξ (cid:48) ( { a i,j } ) = 0 , then, for large n , lim n →∞ α n ,k ( A ) > / and hence it follows that lim n →∞ ˜Γ n ( { a i,j } ) = 0 . Hence ˜Γ n ,n is a height two tower of general algorithms solving { Ξ (cid:48) , Ω (cid:48) } ,a contradiction. Step 2: { Ξ d , Ω d } / ∈ ∆ G . To prove this we can use a slight alteration of the argument in step 1. Replace X by X = l ( N ) ⊕ (cid:76) ∞ j =1 X j and A by A = diag { , , , , , ... } ⊕ ∞ (cid:77) j =1 B j . It is then clear that Ξ d ( A ) = 1 if and only if Ξ (cid:48) ( { a i,j } ) = 1 . OUNDATIONS OF SPECTRAL COMPUTATIONS 41
Step 3: { Ξ d , Ω d } ∈ Σ A . For this we argue similarly to the proof of Theorem3.13 for Ξ d step 2. It wasshown in [8] that there exists a height three arithmetic tower ˜Γ n ,n ,n for the essential spectrum of operatorsin Ω d such that • Each ˜Γ n ,n ,n ( A ) consists of a finite collection of points in the complex plane. • For large n , ˜Γ n ,n ,n ( A ) is eventually constant and equal to ˜Γ n ,n ( A ) . • ˜Γ n ,n ( A ) is increasing with n with limit ˜Γ n ( A ) containing the essential spectrum. The limit ˜Γ n ( A ) is also decreasing with n .Furthermore, it was proven in [8] that for operators in Ω d , there exists a height two arithmetic tower ˆΓ n ,n for computing the spectrum such that • ˆΓ n ,n ( A ) is constant for large n . • For any z ∈ ˆΓ n ( A ) , dist( z, Sp( A )) ≤ − n .Using these, we initially define ζ n ,n ,n ( A ) = { z ∈ ˆΓ n ,n ( A ) : 2 − n − − n ≤ dist( z, ˜Γ n ,n ,n ( A )) } . The arguments in the proof of Theorem 3.13 for Ξ d show that this can be computed in finitely many arith-metic operations and comparisons using the relevant evaluation functions. Note that for large n ζ n ,n ,n ( A ) = { z ∈ ˆΓ n ( A ) : 2 − n − − n ≤ dist( z, ˜Γ n ,n ( A )) } =: ζ n ,n ( A ) . There are now two cases to consider (we use D η ( z ) to denote the open ball of radius η about a point z ): Case 1: Sp d ( A ) ∩ (˜Γ n ( A ) + D − n (0)) c = ∅ . Suppose, for a contradiction, in this case that there exists z m j ∈ ζ n ,m j ( A ) with m j → ∞ . Then, without loss of generality, z m j → z ∈ Sp( A ) . We also have that dist( z m j , ˜Γ n ,m j ( A )) ≥ − n − − m j , which implies that dist( z, ˜Γ n ( A )) ≥ − n and hence z ∈ Sp d ( A ) ∩ (˜Γ n ( A ) + D − n (0)) c , the requiredcontradiction. It follows that ζ n ,n ( A ) is empty for large n . Case 2: Sp d ( A ) ∩ (˜Γ n ( A ) + D − n (0)) c (cid:54) = ∅ . In this case, this set is a finite subset of Sp d ( A ) , { ˆ z , ..., ˆ z m ( n ) } . Each of these points is an isolated point of the spectrum. It follows that there exists z n ∈ ˆΓ n ( A ) with z n → ˆ z and | z n − ˆ z | ≤ − n for large n . Since the ˜Γ n ,n ( A ) are increasing, thisimplies that dist( z n , ˜Γ n ,n ( A )) ≥ dist( z n , ˜Γ n ( A )) ≥ dist(ˆ z , ˜Γ n ( A )) − − n ≥ − n − − n , so that z n ∈ ζ n ,n ( A ) . The same argument holds for points converging to all of { ˆ z , ..., ˆ z m ( n ) } . On theother hand, the argument used in Case 1 shows that any limit points of ζ n ,n ( A ) as n → ∞ are containedin Sp d ( A ) ∩ (˜Γ n ( A ) + D − n (0)) c . It follows that in this case ζ n ,n ( A ) converges to Sp d ( A ) ∩ (˜Γ n ( A ) + B /n (0)) c (cid:54) = ∅ in the Hausdorff metric as n → ∞ .Let N ( A ) ∈ N be minimal such that Sp d ( A ) ∩ (˜Γ N ( A ) + D − N (0)) c (cid:54) = ∅ (recall the discrete spectrumis non-empty for our class of operators). If n > n then set Γ n ,n ,n ( A ) = { } , otherwise consider ζ k,n ,n ( A ) for n ≤ k ≤ n . If all of these are empty then set Γ n ,n ,n ( A ) = { } , otherwise chooseminimal k with ζ k,n ,n ( A ) (cid:54) = ∅ and let Γ n ,n ,n ( A ) = ζ k,n ,n ( A ) . Note that this defines an arithmetictower of algorithms, with Γ n ,n ,n ( A ) non-empty. Since we consider finitely many of the sets ζ k,n ,n ( A ) ,and these are constant for large n , it follows that Γ n ,n ,n ( A ) is constant for large n and constructed inthe same manner with replacing ζ k,n ,n ( A ) by ζ k,n ( A ) . Call this limit Γ n ,n ( A ) .For large n , Γ n ,n ( A ) = ζ n ∨ N ( A ) ,n ( A ) and it follows that lim n →∞ Γ n ,n ( A ) =: Γ n ( A ) = Sp d ( A ) ∩ (˜Γ n ∨ N ( A ) ( A ) + D − n ∨ N ( A ) (0)) c . Hence Γ n ( A ) ⊂ Sp d ( A ) and Γ n ( A ) converges up to Sp d ( A ) in the Hausdorff metric. Step 4: { Ξ d , Ω d } ∈ Σ A . Consider the height three tower, ζ n ,n ,n , defined in step 3. Let A ∈ Ω d andif ζ n ,n ,n ( A ) = ∅ , define ρ n ,n ,n ( A ) = 0 , otherwise define ρ n ,n ,n ( A ) = 1 . The discussion in step 3shows that lim n →∞ lim n →∞ ρ n ,n ,n ( A ) =: ρ n ( A ) = , if Sp d ( A ) ∩ (˜Γ n ( A ) + D − n (0)) c = ∅ , otherwise.Since Sp d ( A ) ∩ (˜Γ n ( A ) + D − n (0)) c increases to Sp d ( A ) , it follows that lim n →∞ ρ n ( A ) = Ξ d ( A ) andthat if ρ n ( A ) = 1 , then Ξ d ( A ) = 1 . Hence, ρ n ,n ,n provides a Σ A tower for { Ξ d , Ω d } . (cid:3)
9. P
ROOF OF T HEOREM ON THE S PECTRAL G AP AND S PECTRAL C LASSIFICATION
Proof of Theorem 3.11 for Ξ gap . Step 1 : { Ξ gap , (cid:98) Ω SA } ∈ Σ A . Let A ∈ (cid:98) Ω SA . Using Corollary 6.9 wecan compute all n eigenvalues of P n AP n to arbitrary precision in finitely many arithmetic operations andcomparisons. Note that it is not completely straightforward to deduce this with the QR algorithm, as one hasto deal with halting criteria in order to achieve the correct precision. Moreover, one must approximate rootsin order to extract the approximate eigenvalues from a potential × matrix block. Thus we use Corollary6.9 instead. In the notation of Lemmas 8.1, and 8.2 (whose analogous results also hold for the possiblyunbounded A ∈ (cid:98) Ω SA ), consider an approximation ≤ l n := µ n − µ n + (cid:15) n , n ≥ , where we have computed µ n − µ n to accuracy | (cid:15) n | ≤ /n using Corollary 6.9 with B = P n AP n . Recallthat for A ∈ (cid:98) Ω SA we restricted the class so that either the bottom of the spectrum is in the discrete spectrumwith multiplicity one, or there is a closed interval in the spectrum of positive measure with the bottom of thespectrum as its left end-point. It follows that l n converges to zero if and only if Ξ gap ( A ) = 0 , otherwise itconverges to some positive number. If n = 1 then set Γ n ,n ( A ) = 1 , otherwise consider the following.Let J n = [0 , / (2 n )] and J n = (1 /n , ∞ ) . Given n ∈ N , consider l k for k ≤ n . If no such k existswith l k ∈ J n ∪ J n then set Γ n ,n ( A ) = 0 . Otherwise, consider k maximal with l k ∈ J n ∪ J n and set Γ n ,n ( A ) = 0 if l k ∈ J n and Γ n ,n ( A ) = 1 if l k ∈ J n . The sequence l n → c ≥ for some number c . The separation of the intervals J n and J n , ensures that l n cannot be in both intervals infinitely often as n → ∞ and hence the first limit Γ n ( A ) := lim n →∞ Γ n ,n ( A ) exists. If c = 0 , then Γ n ( A ) = 0 butif c > then there exists n with /n < c and hence for large n , l n ∈ J n . It follows in this case that Γ n ( A ) = 1 and we also see that if Γ n ( A ) = 1 then Ξ gap ( A ) = 1 . Hence Γ n ,n provides a Σ A tower. Step 2 : { Ξ gap , (cid:98) Ω D } / ∈ ∆ G . We argue by contradiction and assume the existence of a height one tower, Γ n converging to Ξ gap . The method of proof follows the same lines as before. For every A and n thereexists a finite number N ( A, n ) ∈ N such that the evaluations from Λ Γ n ( A ) only take the matrix entries A ij = (cid:104) Ae j , e i (cid:105) with i, j ≤ N ( A, n ) into account. List the rationals in (0 , without repetition as d , d , ... .We consider the operators A m = diag { d , d , ..., d m } ∈ C m × m , B m = diag { , , ..., } ∈ C m × m and C = diag { , , ... } . Let A = ∞ (cid:77) m =1 ( B k m ⊕ A k m ) , where we choose an increasing sequence k m inductively as follows. In what follows, all operators consideredare easily seen to be in (cid:98) Ω D .Set k = 1 and suppose that k , ..., k m have been chosen with the property that upon defining ζ p := min { d r : 1 ≤ r ≤ k p } , OUNDATIONS OF SPECTRAL COMPUTATIONS 43 we have ζ p > ζ p +1 for p = 1 , ..., m − . Sp( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = { d , d , ..., d m , } has ζ m the minimum of its spectrum and an isolated eigenvalue of multiplicity , hence Ξ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = “Yes” . It follows that there exists some n m ≥ m such that if n ≥ n m then Γ n ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = “Yes” . Now let k m +1 ≥ max { N ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C, n m ) , k m + 1 } with ζ m > ζ m +1 . The sameargument used in the proof of Theorem 3.9 shows that Γ n m ( A ) = Γ n m ( B k ⊕ A k ⊕ ... ⊕ B k m ⊕ A k m ⊕ C ) = “Yes”. But Sp( A ) = [0 , is gappless and so must have lim n →∞ (Γ n ( A )) = “No”, a contradiction. (cid:3) Proof of Theorem 3.11 for Ξ class . By restricting (cid:101) Ω D to (cid:98) Ω D and composing with the map ρ : { , , , } → { , } ,ρ (1) = 1 , ρ (2) = ρ (3) = ρ (4) = 0 , it is clear that the result for Ξ gap implies { Ξ class , (cid:101) Ω f SA } , { Ξ class , (cid:101) Ω D } / ∈ ∆ G . Since (cid:101) Ω D ⊂ (cid:101) Ω f SA , we need only construct a Π A tower for { Ξ class , (cid:101) Ω f SA } .Let A ∈ (cid:101) Ω f SA . For a given n , set B n = P n AP n and in the notation of Lemmas 8.2 and 8.1, let ≤ l jn := µ nj +1 − µ n + (cid:15) jn , for j < n. where we again have computed µ nj +1 − µ n to accuracy (cid:12)(cid:12) (cid:15) jn (cid:12)(cid:12) ≤ /n using only finitely many arithmeticoperations and comparisons by Corollary 6.9. Ξ class ( A ) = 1 if and only if l n converges to a positiveconstant as n → ∞ and Ξ class ( A ) = 2 if and only if l n converges to zero as n → ∞ but there exists j with l jn convergent to a positive constant.Note that we can use the algorithm, denoted ˆΓ n , to compute the spectrum presented in §
6, with errorfunction denoted by E ( n, · ) converging uniformly on compact subsets of C to the true error from above(again with the choice of g ( x ) = x since the operator is normal). Setting a n ( A ) = min x ∈ ˆΓ n ( A ) { x + E ( n, x ) } , we see that a n ( A ) ≥ a ( A ) := inf x ∈ Sp( A ) { x } and that a n ( A ) → a ( A ) . Now consider b n ,n ( A ) = min { E ( k, a k ( A ) + 1 /n ) + 1 /k : 1 ≤ k ≤ n } then b n ,n ( A ) is positive and decreasing in n so converges to some limit b n ( A ) . Lemma 9.1.
Let A ∈ (cid:101) Ω f SA and c n ,n ( A ) = E ( n , a n ( A ) + 1 /n ) + 1 /n , then lim n →∞ c n ,n ( A ) =: c n ( A ) = dist( a + 1 /n , Sp( A )) . Furthermore, if Ξ class ( A ) (cid:54) = 4 then for large n it follows that c n ( A ) = b n ( A ) = 1 /n .Proof of Lemma 9.1. We know that a n ( A ) + 1 /n converges to a ( A ) + 1 /n as n → ∞ . Furthermore, dist( z, Sp( A )) is continuous in z and E ( n , z ) converges uniformly to dist( z, Sp( A )) on compact subsets of C . Hence, the limit c n ( A ) exists and is equal to dist( a ( A )+1 /n , Sp( A )) . It is clear that b n ( A ) ≤ c n ( A ) .Suppose now that Ξ class ( A ) (cid:54) = 4 , then for large n , say bigger than some N , and for large enough n , E ( n , a n ( A ) + 1 /n ) ≥ dist( a n ( A ) + 1 /n , Sp( A ))= | a n ( A ) + 1 /n − a ( A ) |≥ /n = dist( a ( A ) + 1 /n , Sp( A )) . Now choose n large such that the above inequality holds and /n ≤ /N . Then b n ,n ( A ) ≥ /n .Taking limits finishes the proof. (cid:3) If n ≥ n then set Γ n ,n ( A ) = 1 . Otherwise, for ≤ j ≤ n , let k jn ,n be maximal with ≤ k jn ,n < n such that l jk jn ,n ∈ J n ∪ J n if such k jn ,n exist, where J n and J n are as in the proof for Ξ gap . If k n ,n exists with l k n ,n ∈ J n then set Γ n ,n ( A ) = 1 . Otherwise, if any of k mn ,n exists with l mk mn ,n ∈ J n for ≤ m ≤ n then set Γ n ,n ( A ) = 2 . Suppose that neither of these two cases hold. In thiscase compute b n ,n ( A ) . If b n ,n ( A ) ≥ /n then set Γ n ,n ( A ) = 3 , otherwise set Γ n ,n ( A ) = 4 . Wenow must show this provides a Π A tower solving our problem.First we show convergence of the first limit. Fix n and consider n large. The separation of the intervals J n and J n ensures that each sequence { l jn } n ∈ N cannot visit each interval infinitely often. Since b n ,n ( A ) is non-increasing in n , we also see that the question whether b n ,n ( A ) ≥ /n eventually has a constantanswer. These observations ensure convergence of the first limit Γ n ( A ) = lim n →∞ Γ n ,n ( A ) .If Ξ class ( A ) = 1 then for large n , l n must eventually be in J n and hence Γ n ( A ) = 1 . It is also clearthat if Γ n ( A ) = 1 then l n converges to a positive constant, which implies Ξ class ( A ) = 1 . If Ξ class ( A ) = 2 then for large n , l mn eventually lies in J n for some ≤ m ≤ n , but l n eventually in J n . It follows that Γ n ( A ) = 2 . If Γ n ( A ) = 2 , then we know that there exists some l mn convergent to l ≥ /n and hence weknow Ξ class ( A ) is either or .Now suppose that Ξ class ( A ) = 3 , then for fixed n and any ≤ m ≤ n , l mn eventually lies in J n andhence our lowest level of the tower must eventually depend on whether b n ,n ( A ) ≥ /n . From Lemma9.1, b n ( A ) = c n ( A ) = 1 /n for large n . It follows that for large n , b n ( A ) ≥ /n for all n and Γ n ( A ) = 3 . Furthermore, if Γ n ( A ) = 3 then we know that c n ( A ) ≥ b n ( A ) ≥ /n , which implies Ξ class ( A ) (cid:54) = 4 . Finally, note that if Ξ class ( A ) = 4 but there exists n with Γ n ( A ) (cid:54) = 4 then the aboveimplies the contradiction Ξ class ( A ) (cid:54) = 4 . The partial converses proven above imply Γ n ,n realises the Π A classification. (cid:3)
10. C
OMPUTATIONAL E XAMPLES
In this section, we demonstrate that, as well as being optimal from a foundations point of view, the algo-rithms constructed in this paper are usable and can be efficiently implemented for large scale computations.The algorithms constructed in this paper have desirable convergence properties, with some converging mono-tonically or being eventually constant as captured by the Σ / Π classification. They are also completely localand hence parallelisable . Pseudocodes are given in Appendix A. Remark 10.1.
Although we only stated results for the graph case, l ( V ( G )) , in § § § § H . Once a basis is chosen (so that matrix elements make sense) we can introduceconcepts like bounded dispersion etc.10.1. Polynomial coefficients.
We first consider computing spectra and pseudospectra for partial differen-tial operators of the form(10.1) T = P ( x , ..., x d , ∂ , ..., ∂ d ) for a polynomial P . In this case, the algorithms (which use a Hermite basis in the proof) reduce to computingthe spectrum/pseudospectrum of infinite matrices A acting on l ( N ) . From the comments in Example 6.13and recurrence relations for Hermite functions, we can choose a basis such that f ( n ) − n ∼ Cn ( d − /d ,where f is the dispersion function and C ( d ) a constant, such that f also describes the off-diagonal sparsitystructure of A in the sense that A n,k = A k,n = 0 if k > f ( n ) . Hence, this section also showcases thealgorithms presented in § verified rigorously with interval arithmetic . OUNDATIONS OF SPECTRAL COMPUTATIONS 45
Potential Exact n = 500 n = 1000 V − − ± × − − ± − V − − ± × − − ± − V .
375 0 . ± . × − . ± × − V .
125 1 . ± . × − . ± . × − T ABLE
2. Test run of algorithm on some potentials with known eigenvalues. Note that wequickly converge to the eigenvalue with error bounds computed by the algorithm (through
DistSpec ) and using interval arithmetic.10.1.1.
Anharmonic oscillators.
First, consider operators of the form H = − ∆ + V ( x ) = − ∆ + d (cid:88) j a j x j + d (cid:88) j,k =1 b j,k x j x k + (cid:88) α ∈ Z d ≥ , | α |≤ M c ( α ) x α , where a j , b j,k , c ( α ) ∈ R and the multi-indices α are chosen such that (cid:80) | α |≤ M c ( α ) x α is bounded frombelow. The fact that H is essentially self-adjoint follows from the Faris–Lavine theorem [92, TheoremX.28]. Anharmonic oscillators have attracted interest in quantum research for over three decades [12, 13,64, 111] and amongst their uses are approximations of potentials near stationary points. The problem ofdeveloping efficient algorithms to compute their spectra has received renewed interest due to advances inasymptotic analysis and symbolic computing algebra [6, 66, 107]. The methods in the cited works are richand diverse, but lack uniformity. We show that we can obtain error control for general anharmonic operatorsin a computationally efficient manner.The algorithm Γ n ( A ) for the spectrum is described by the routine CompSpecUB , shown as pseudocodein Appendix A. This relies on the approximation to (cid:107) R ( z, A ) (cid:107) − in Theorem 6.7 given by the routine DistSpec . Throughout, we have used the fact that
DistSpec can be implemented using only finitelymany arithmetic operations and comparisons.We begin with comparisons to some known results in one dimension, calculated using super-symmetricquantum mechanics [38]: V ( x ) = x − x + x E = − ,V ( x ) = 4 x − x + x E = − ,V ( x ) = (105 / x − (43 / x + x − x + x E = 3 / ,V ( x ) = (169 / x − (59 / x + x − x + x E = 9 / . Following the physicists’ convention, if the spectrum is discrete and bounded below, we list the energylevels as E ≤ E ≤ E ≤ ... . Note that other methods such as finite section (of the correspondingmatrices A constructed using Hermite functions) will converge in this case, but do not provide the sharp Σ classification. We found that the grid resolution of the search routine and the search accuracy for thesmallest singular values, not the matrix size, were the main deciding factors in the error bound. Clearly, oncewe know roughly where the eigenvalues are, we can speed up computations using the fact that the algorithmis local. Furthermore, the search routine’s computational time only grows logarithmically in its precision.Hence we set the grid spacing and the spacing of the search routine to be n . Table 2 shows the results;all values were computed rapidly using a local search grid. Note that we quickly gain convergence and thatthe error bounds become the precision of the search routine in DistSpec (namely, n ). In this simpleexample, the output happens to agree precisely with the eigenvalues since they lie on the search grid.Next, we demonstrate how the algorithm can be used in more than one dimension. We consider H = − ∆ + x x , n = 1000 DistSpec n = 2000 z n = 5000 DistSpec
Finite SectionTrue Eigenvalues -16 -14 -12 -10 -8 -6 -4 -2 F IGURE
2. Left: The convergence of the algorithm (shown as
DistSpec ) and finitesection to the true eigenvalues on the interval [0 , . Note that points with reliable finitesection eigenvalues correspond to points where the estimate of the resolvent norm is wellresolved. Right: Error bounds computed using DistSpec and interval arithmetic (for anadaptive grid spacing). -5 0 5 10 15 20 25 30
Re( z ) -15-10-5051015 I m ( z ) . . . . . . . Resolvent NormEigenvalues
Re( z ) -2-1.5-1-0.500.511.52 I m ( z )
22 55 1010 2020 5050 100100200200 5005001000100020002000 50005000100001000020000200005000050000100000100000200000200000 500000500000 10000001000000 20000002000000 5000000500000010000000100000002000000020000000
Resolvent Norm F IGURE
3. Left: Calculated pseudospectrum for the imaginary cubic oscillator. Notethe clear presence of eigenvalues. Right: Calculated pseudospectrum for imaginary Airyoperator. Both figures were produced with n = 1000 .which is a classic example of a potential that does not blow up at ∞ in every direction, yet still inducesan operator with compact resolvent [99]. Figure 2 shows the convergence of the estimate of (cid:107) R ( z, H ) (cid:107) − from above, as well as finite section estimates. As expected from variational methods, the finite sectionmethod produces eigenvalues converging to the true eigenvalues from above (there is no essential spectrumand the operator is positive). Furthermore, the areas where DistSpec has converged correspond to areaswhere finite section has converged. We also show rigorous error bounds computed using
DistSpec fordifferent n for the first five eigenvalues. These are computed using an adaptive grid spacing to resolve thelocal minima of the approximation of (cid:107) R ( z, H ) (cid:107) − using rectangular truncations.10.1.2. Pseudospectra and PT symmetry. We now turn to pseudospectra and consider PT -symmetric non-self-adjoint operators T (and for which compactly supported smooth functions form a core of T and T ∗ [54]).The first example is the imaginary cubic oscillator defined formally (in one dimension) by H = − d /dx + ix . OUNDATIONS OF SPECTRAL COMPUTATIONS 47
Potential
V E E E E E cos( x ) 1 . . . . . x ) 0 . . . . . − x ) 1 . . . . . x ) − . . . . . T ABLE
3. Computed eigenvalues for different potentials (first five shown). Each eigen-value E n , computed with an error bound at most − via DistSpec , is a shift of theharmonic oscillator eigenvalue n + 1 This operator is the most studied example of a PT -symmetric operator [10, 11], as well as appearing instatistical physics and quantum field theory [65]. It is known that the resolvent is compact [37] with alleigenvalues simple and residing in R ≥ [52,105]. The eigenvectors are complete but do not form a Riesz ba-sis [98]. Figure 3 (left) shows the computed pseudospectrum computed using n = 1000 . This demonstratesthe instability of the spectrum of the operator.Next, we consider the imaginary Airy operator H = − d /dx + ix, since this is known to have empty spectrum [74], demonstrating that the algorithm is effective in this casealso. Note that any finite section method will overestimate the pseudospectrum due to the presence of falseeigenvalues. H is PT -symmetric and has compact resolvent. The resolvent norm (cid:107) R ( z, H ) (cid:107) only de-pends on the real part of z and blows up exponentially as Re( z ) → + ∞ . We have shown the computedpseudospectrum for n = 1000 in Figure 3 (right).We do not need to discretise anything to apply the above method. Up to numerical errors in the testingof positive definiteness, all computed pseudospectra are guaranteed to be inside the correct pseudospectra.In fact, in our case, we checked results using interval arithmetic, and hence the output is 100% reliable.This is in contrast to the numerical experiments conducted in [47], where the operator is discretised. Itis also easy to construct examples where discretisations fail dramatically, either not capturing the wholespectrum or suffering from spectral pollution (even without spectral pollution - figuring out which parts ofcomputations are trustworthy can be very difficult for finite section and related methods [113]). Algorithmslike PseudoSpec are a useful tool to test the reliability of such outputs.10.2.
Partial differential operators with general coefficients.
Perturbed harmonic oscillator.
As a first set of examples, we consider T = − ∆ + x + V ( x ) , on L ( R ) , where V is a bounded potential. Such operators have discrete spectra, however, the perturbation V causes the eigenvalues to shift relative to the classical harmonic oscillator (whose spectrum is the set ofodd positive integers). Table 3 shows the first five eigenvalues for a range of potentials. Each entry in thetable is computed with an error bound at most − provided by DistSpec .10.2.2.
Fourth-order operator.
Next, we consider the operator T λ = d dx + (cid:18) − i ddx + x (cid:19) + 2 λx + λ x , on L ( R ) , as an example with gaps in the essential spectrum. Figure 4 shows a portion of the spectrum,as well as the output of finite section, using basis functions. The maximum error bound provided by Meaning [ H , PT ] = 0 with ( P f )( x ) = f ( − x ) and ( T f )( x ) = f ( x ) . Finite Section
2 4 6
Spectral Estimate Spectral P oll u tion F IGURE
4. Left: Output of finite section. Right: Output of
CompSpecUB . Both use basis functions and the spectrum of T is shown as red stars. DistSpec is bounded by − . Finite section produces heavy spectral pollution in the gaps of the essentialspectrum. The spectrum for λ = 0 is shown as red stars, and consists of isolated eigenvalues of infinitemultiplicity. As λ increases, these fan out to produce the essential spectra shown.10.3. Example for discrete Spectra.
Although it is hard to analyse the convergence of a height two tower,we can take advantage of the extra structure in this problem. The routine
DiscreteSpec in AppendixA computes Γ n ,n ( A ) such that lim n →∞ Γ n ,n ( A ) is a finite subset of Sp d ( A ) . Furthermore, for each z ∈ Sp d ( A ) , there is at most one point in z n ∈ Γ n ,n ( A ) approximating z . We can use the routine DistSpec to gain an error bound of dist( z n , Sp( A )) , which, for large n , will be equal to | z − z n | since z is isolated. As we increase n , more and more of the discrete spectrum (in general portions nearer theessential spectrum) are approximated.Our example is the Almost Mathieu Operator on l ( Z ) given by ( H α x ) n = x n − + x n +1 + 2 λ cos(2 πnα ) x n , where we set λ = 1 (critical coupling). For rational choices of α , the operator is periodic and its spectrumis purely absolutely continuous. For irrational α the spectrum is a Cantor set (Ten Martini Problem). Togenerate a discrete spectrum, we add a perturbation of the potential of the form(10.2) V ( n ) = V n / ( | n | + 1) , where V n are independent and uniformly distributed in [ − , . The perturbation is compact so preservesthe essential spectrum. This type of problem is well studied in the more general setting of Jacobi operators[75, 106].Figure 5 shows a typical result for a realisation of the random potential. The figure shows the output offinite section and our algorithm (with a uniform error bound of − ) for computing the total spectrum. Wehave also shown the output of DiscreteSpec , which separates the discrete spectrum from the essentialspectrum. For each α we took n large enough for an expected limit inclusion Sp d ( A ) ⊂ Γ n ( A ) + B . (0) (obtained by comparing with the output of the height two tower for computing the essential spectrum). Tak-ing n larger caused sharper inclusion bounds. Additionally, we confirmed the accuracy of the results usinga height one tower to compute the spectrum with and without the random potential. Note that it is difficultto detect spectral pollution when using finite section with the additional perturbation (10.2). In contrast, DiscreteSpec computes the discrete spectrum without spectral pollution and allows us to separate thediscrete spectrum from the essential spectrum.
OUNDATIONS OF SPECTRAL COMPUTATIONS 49 F IGURE
5. Top: Output of finite section. Spectral pollution detected by our algorithmis shown as red crosses. Bottom: Output of
DiscreteSpec and the splitting into theessential spectrum and the discrete spectrum. The output captures the discrete spectrumdown to a distance ≈ . away from the essential spectrum, which can be made smallerfor larger n .The error bounds provided by DistSpec (applied to the output of
DiscreteSpec ) can also be trans-lated into computing approximates of the eigenvectors of an operator A corresponding to the discrete spec-trum with an error bound in the following manner. The routine ApproxEigenvector in Appendix Acomputes a vector x n of norm ≈ such that (in this case taking δ ↓ , c n = 0 ) (cid:107) ( A − z n I ) x n (cid:107) ≤ DistSpec ( A, n , f ( n ) , z n ) . We write x n = x dn + y n , where x dn is an eigenvector of A with eigenvalue z , and y n is perpendicular tothe eigenspace associated with z and z n → z . It follows that (cid:107) ( A − zI ) y n (cid:107) ≤ | z − z n | + DistSpec ( A, n , f ( n ) , z n ) ≤ × DistSpec ( A, n , f ( n ) , z n ) , for large n . But A − zI is bounded below on the orthogonal complement of the eigenspace, with lowerbound dist( z, Sp( A ) \{ z } ) . Hence, (cid:107) y n (cid:107) ≤ × DistSpec ( A, n , f ( n ) , z n )dist( z, Sp( A ) \{ z } ) for large n . This also bounds the l distance of x n to the eigenspace and can be estimated by approximatingthe spectrum of A . It is also straightforward to adjust this procedure to eigenvalues of multiplicity greaterthan and approximate the whole eigenspace. We note that for this example, all eigenvalues were found tohave multiplicity as expected for a random perturbation. Finally, the method of computing eigenvectorsand error bounds can also be used for the unbounded case when z lies in the discrete spectrum.R EFERENCES[1] W. Arveson. Discretized CCR algebras.
Journal of Operator Theory , 26(2):225–239, 1991.[2] W. Arveson. Improper filtrations for C ∗ -algebras: spectra of unilateral tridiagonal operators. Acta Sci. Math. (Szeged) , 57(1-4):11–24, 1993.[3] W. Arveson. Noncommutative spheres and numerical quantum mechanics. In
Operator algebras, mathematical physics, andlow-dimensional topology , volume 5 of
Res. Notes Math. , pages 1–10. A K Peters, Wellesley, MA, 1993.[4] W. Arveson. C ∗ -algebras and numerical linear algebra. Journal of Functional Analysis , 122(2):333–360, 1994.[5] W. Arveson. The role of C ∗ -algebras in infinite-dimensional numerical linear algebra. In C ∗ -algebras: 1943–1993 (San Antonio,TX, 1993) , volume 167 of Contemp. Math. , pages 114–129. Amer. Math. Soc., Providence, RI, 1994.[6] T. Barakat. The asymptotic iteration method for the eigenenergies of the anharmonic oscillator potential V ( x ) = Ax α + Bx . Physics Letters A , 344(6):411–417, 2005.[7] A. Bastounis, A. C. Hansen, and V. Vlacic. On computational barriers and paradoxes in estimation, regularisation and learning.
Preprint , 2018.[8] J. Ben-Artzi, M. J. Colbrook, A. C. Hansen, O. Nevanlinna, and M. Seidel. Computing Spectra – On the Solvability ComplexityIndex hierarchy and towers of algorithms. arXiv:1508.03280v5 , 2020.[9] J. Ben-Artzi, M. Marletta, and F. R¨osler. Computing scattering resonances. arXiv:2006.03368 , 2020.[10] C. M. Bender and S. Boettcher. Real spectra in non-Hermitian Hamiltonians having
P T symmetry.
Physical Review Letters ,80(24):5243, 1998.[11] C. M. Bender, D. C. Brody, and H. F. Jones. Complex extension of quantum mechanics.
Physical Review Letters , 89(27):270401,2002.[12] C. M. Bender and S. A. Orszag.
Advanced mathematical methods for scientists and engineers I: Asymptotic methods and pertur-bation theory . Springer Science & Business Media, 2013.[13] C. M. Bender and T. T. Wu. Anharmonic oscillator. II. A study of perturbation theory in large order.
Physical Review D , 7(6):1620,1973.[14] L. Blum, F. Cucker, M. Shub, and S. Smale.
Complexity and Real Computation . Springer-Verlag New York, Inc., Secaucus, NJ,USA, 1998.[15] L. Blum, M. Shub, and S. Smale. On a theory of computation and complexity over the real numbers: NP-completeness, recursivefunctions and universal machines.
American Mathematical Society. Bulletin. , 21(1):1–46, 1989.[16] M. Bl¨umlinger and R. F. Tichy. Topological algebras of functions of bounded variation I.
Manuscripta Mathematica , 65(2):245–255, 1989.[17] D. Boffi, F. Brezzi, and L. Gastaldi. On the problem of spurious eigenvalues in the approximation of linear elliptic problems inmixed form.
Mathematics of Computation , 69(229):121–140, 2000.[18] D. Boffi, R. G. Duran, and L. Gastaldi. A remark on spurious eigenvalues in a square.
Appl. Math. Lett. , 12(3):107–114, 1999.[19] S. B¨ogli, B. M. Brown, M. Marletta, C. Tretter, and M. Wagenhofer. Guaranteed resonance enclosures and exclosures for atomsand molecules.
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , 470(2171):20140488,2014.[20] E. Bombieri, S. Cook, P. Deligne, C. Fefferman, J. Gray, A. Jaffe, J. Milnor, A. Wiles, and E. Witten. The Millennium PrizeProblems.
CMI/AMS , 2006.[21] A. B¨ottcher. Pseudospectra and singular values of large convolution operators.
Journal of Integral Equations and Applications ,6(3):267–301, 1994.[22] A. B¨ottcher. Infinite matrices and projection methods. In
Lectures on operator theory and its applications (Waterloo, ON, 1994) ,volume 3 of
Fields Inst. Monogr. , pages 1–72. Amer. Math. Soc., Providence, RI, 1996.[23] A. B¨ottcher, H. Brunner, A. Iserles, and S. P. Nørsett. On the singular values and eigenvalues of the Fox-Li and related operators.
New York J. Math. , 16:539–561, 2010.
OUNDATIONS OF SPECTRAL COMPUTATIONS 51 [24] A. B¨ottcher, S. Grudsky, and A. Iserles. Spectral theory of large Wiener–Hopf operators with complex-symmetric kernels andrational symbols.
Mathematical Proceedings of the Cambridge Philosophical Society , 151(1):161–191, 2011.[25] A. B¨ottcher and B. Silbermann. The finite section method for Toeplitz operators on the quarter-plane with piecewise continuoussymbols.
Mathematische Nachrichten , 110:279–291, 1983.[26] A. B¨ottcher and B. Silbermann.
Introduction to large truncated Toeplitz matrices . Universitext. Springer-Verlag, New York,1999.[27] A. B¨ottcher and B. Silbermann.
Analysis of Toeplitz operators . Springer Monographs in Mathematics. Springer-Verlag, Berlin,second edition, 2006.[28] N. Brown. Invariant means and finite representation theory of C*-algebras.
Memoirs of the American Mathematical Society , 184,05 2003.[29] N. P. Brown. AF embeddings and the numerical computation of spectra in irrational rotation algebras.
Numer. Funct. Anal.Optim. , 27(5-6):517–528, 2006.[30] N. P. Brown. Quasi-diagonality and the finite section method.
Mathematics of Computation , 76(257):339–360, 2007.[31] N. P. Brown, K. Dykema, and D. Shlyakhtenko. Topological entropy of free product automorphisms.
Acta Mathematica ,189(1):1–35, 2002.[32] H. Brunner, A. Iserles, and S. P. Nørsett. The spectral problem for a class of highly oscillatory Fredholm integral operators.
IMAJournal of Numerical Analysis , 30(1):108–130, 12 2008.[33] H. Brunner, A. Iserles, and S. P. Nørsett. The computation of the spectra of highly oscillatory Fredholm integral operators.
J.Integral Equations Applications , 23(4):467–519, 12 2011.[34] A. Buffa, P. Houston, and I. Perugia. Discontinuous Galerkin computation of the Maxwell eigenvalues on simplicial meshes.
Journal of Computational and Applied Mathematics , 204(2):317–333, 2007.[35] A. Buffa and I. Perugia. Discontinuous Galerkin approximation of the Maxwell eigenproblem.
SIAM Journal on NumericalAnalysis , 44(5):2198–2226, 2006.[36] A. Buffa, I. Perugia, and T. Warburton. The mortar-discontinuous Galerkin method for the 2D Maxwell eigenproblem.
J. Sci.Comput. , 40(1-3):86–114, 2009.[37] E. Caliceti, S. Graffi, and M. Maioli. Perturbation theory of odd anharmonic oscillators.
Communications in MathematicalPhysics , 75(1):51–66, 1980.[38] R. Chaudhuri and M. Mondal. Improved Hill determinant method: General approach to the solution of quantum anharmonicoscillators.
Physical Review A , 43(7):3241, 1991.[39] S. H. Christiansen and R. Winther. On variational eigenvalue approximation of semidefinite operators.
IMA J. Numer. Anal. ,33(1):164–189, 2013.[40] M. J. Colbrook. Computing spectral measures and spectral types. arXiv:1908.06721v2 , 2019.[41] M. J. Colbrook. The foundations of spectral computations via the solvability complexity index hierarchy: Part II. arXiv:1908.09598 , 2019.[42] M. J. Colbrook and A. C. Hansen. On the infinite-dimensional QR algorithm.
Numerische Mathematik , 143(1):17–83, 2019.[43] M. J. Colbrook, B. Roman, and A. C. Hansen. How to compute spectra with error control.
Physical Review Letters ,122(25):250201, 2019.[44] T. S. Cubitt, D. Perez-Garcia, and M. M. Wolf. Undecidability of the spectral gap.
Nature , 528(7581):207, 2015.[45] F. Cucker. The arithmetical hierarchy over the reals.
J. Logic Comput. , 2(3):375–395, 1992.[46] E. B. Davies. Spectral enclosures and complex resonances for general self-adjoint operators.
LMS J. Comput. Math. , 1:42–74,1998.[47] E. B. Davies. Pseudo-spectra, the harmonic oscillator and complex resonances.
R. Soc. Lond. Proc. Ser. A Math. Phys. Eng. Sci. ,455(1982):585–599, 1999.[48] E. B. Davies. A hierarchical method for obtaining eigenvalue enclosures.
Math. Comp. , 69(232):1435–1455, 2000.[49] P. Deift, J. Demmel, L.-C. Li, and C. Tomei. The bidiagonal singular value decomposition and Hamiltonian mechanics.
SIAMjournal on numerical analysis , 28(5):1463–1516, 1991.[50] P. Deift, L. C. Li, and C. Tomei. Toda flows with infinitely many variables.
Journal of Functional Analysis , 64(3):358–402, 1985.[51] T. Digernes, V. S. Varadarajan, and S. S. Varadhan. Finite approximations to quantum systems.
Rev. Math. Phys. , 6(4):621–648,1994.[52] P. Dorey, C. Dunning, and R. Tateo. Spectral equivalences, Bethe ansatz equations, and reality properties in
P T -symmetricquantum mechanics.
Journal of Physics A: Mathematical and General , 34(28):5679, 2001.[53] P. Doyle and C. McMullen. Solving the quintic by iteration.
Acta Mathematica , 163(3-4):151–180, 1989.[54] D. E. Edmunds and W. D. Evans.
Spectral theory and differential operators . Oxford Mathematical Monographs. The ClarendonPress, Oxford University Press, New York, 1987.[55] C. Fefferman and L. Seco. On the energy of a large atom.
Bull. Amer. Math. Soc. (N.S.) , 23(2):525–530, 1990.[56] C. Fefferman and L. Seco. Eigenvalues and eigenfunctions of ordinary differential operators.
Adv. Math. , 95(2):145–305, 1992. [57] C. Fefferman and L. Seco. Aperiodicity of the Hamiltonian flow in the Thomas-Fermi potential.
Rev. Mat. Iberoamericana ,9(3):409–551, 1993.[58] C. Fefferman and L. Seco. The density in a one-dimensional potential.
Adv. Math. , 107(2):187–364, 1994.[59] C. Fefferman and L. Seco. The eigenvalue sum for a one-dimensional potential.
Adv. Math. , 108(2):263–335, 1994.[60] C. Fefferman and L. Seco. On the Dirac and Schwinger corrections to the ground-state energy of an atom.
Adv. Math. , 107(1):1–185, 1994.[61] C. Fefferman and L. Seco. The density in a three-dimensional radial potential.
Adv. Math. , 111(1):88–161, 1995.[62] C. Fefferman and L. Seco. The eigenvalue sum for a three-dimensional radial potential.
Adv. Math. , 119(1):26–116, 1996.[63] C. Fefferman and L. Seco. Interval arithmetic in quantum mechanics. In
Applications of interval computations (El Paso, TX,1995) , volume 3 of
Appl. Optim. , pages 145–167. Kluwer Acad. Publ., Dordrecht, 1996.[64] F. M. Fern´andez, Q. Ma, and R. Tipping. Tight upper and lower bounds for energy eigenvalues of the Schr¨odinger equation.
Physical Review A , 39(4):1605, 1989.[65] M. E. Fisher. Yang-Lee edge singularity and φ field theory. Physical Review Letters , 40(25):1610, 1978.[66] P. J. Gaudreau, R. M. Slevinsky, and H. Safouhi. Computing energy eigenvalues of anharmonic oscillators using the doubleexponential Sinc collocation method.
Annals of Physics , 360:520–538, 2015.[67] H. H. Goldstine, F. J. Murray, and J. von Neumann. The Jacobi method for real symmetric matrices.
J. ACM , 6(1):59–96, Jan.1959.[68] O. Golinelli, T. Jolicoeur, and R. Lacaze. Finite-lattice extrapolations for a Haldane-gap antiferromagnet.
Physical Review B ,50(5):3037, 1994.[69] T. Hales, M. Adams, G. Bauer, T. D. Dang, J. Harrison, L. T. Hoang, C. Kaliszyk, V. Magron, S. McLaughlin, T. T. Nguyen,Q. T. Nguyen, T. Nipkow, S. Obua, J. Pleso, J. Rute, A. Solovyev, T. H. A. Ta, N. T. Tran, T. D. Trieu, J. Urban, K. Vu, andR. Zumkeller. A formal proof of the Kepler conjecture.
Forum Math. Pi , 5:e2, 29, 2017.[70] T. C. Hales. A proof of the Kepler conjecture.
Annals of Mathematics (2) , 162(3):1065–1185, 2005.[71] A. C. Hansen. On the approximation of spectra of linear operators on Hilbert spaces.
Journal of Functional Analysis ,254(8):2092–2126, 2008.[72] A. C. Hansen. Infinite-dimensional numerical linear algebra: theory and applications.
Proc. R. Soc. Lond. Ser. A Math. Phys.Eng. Sci. , 466(2124):3539–3559, 2010.[73] A. C. Hansen. On the solvability complexity index, the n -pseudospectrum and approximations of spectra of operators. Journalof the American Mathematical Society , 24(1):81–124, 2011.[74] B. Helffer.
Spectral theory and its applications , volume 139 of
Cambridge Studies in Advanced Mathematics . Cambridge Uni-versity Press, Cambridge, 2013.[75] D. Hundertmark and B. Simon. Lieb-Thirring inequalities for Jacobi matrices.
Journal of Approximation Theory , 118(1):106–130, 2002.[76] J. Indritz. An inequality for Hermite polynomials.
Proc. Amer. Math. Soc. , 12:981–983, 1961.[77] T. Kato. On the upper and lower bounds of eigenvalues.
Journal of the Physical Society of Japan , 4(4-6):334–339, 1949.[78] A. Laptev and Y. Safarov. Szeg˝o type limit theorems.
Journal of Functional Analysis , 138(2):544–559, 1996.[79] C. Lubich.
From quantum to classical molecular dynamics: reduced models and numerical analysis . Zurich Lectures in AdvancedMathematics. European Mathematical Society (EMS), Z¨urich, 2008.[80] B. Malcolm Brown, M. Langer, M. Marletta, C. Tretter, and M. Wagenhofer. Eigenvalue enclosures and exclosures for non-self-adjoint problems in hydrodynamics.
LMS Journal of Computation and Mathematics , 13:65–81, 2010.[81] M. Marletta. Neumann-Dirichlet maps and analysis of spectral pollution for non-self-adjoint elliptic PDEs with real essentialspectrum.
IMA J. Numer. Anal. , 30(4):917–939, 2010.[82] M. Marletta and R. Scheichl. Eigenvalues in spectral gaps of differential operators.
Journal of Spectral Theory , 2(3):293–320,2012.[83] C. McMullen. Families of rational maps and iterative root-finding algorithms.
Annals of Mathematics (2) , 125(3):467–493, 1987.[84] C. McMullen. Braiding of the attractor and the failure of iterative algorithms.
Invent. Math. , 91(2):259–272, 1988.[85] H. Niederreiter.
Random number generation and quasi-Monte Carlo methods , volume 63 of
CBMS-NSF Regional ConferenceSeries in Applied Mathematics . Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1992.[86] S. Olver. ApproxFun.jl v0.8. github (online) https://github.com/JuliaApproximation/ApproxFun.jl , 2018.[87] S. Olver and A. Townsend. A fast and well-conditioned spectral method.
SIAM Review , 55(3):462–489, 2013.[88] S. Olver and A. Townsend. A Practical Framework for Infinite-dimensional Linear Algebra. In
Proceedings of the 1st FirstWorkshop for High Performance Technical Computing in Dynamic Languages , HPTCDL ’14, pages 57–62, Piscataway, NJ,USA, 2014. IEEE Press.[89] S. Olver and M. Webb. SpectralMeasures.jl. github (online) https://github.com/JuliaApproximation/SpectralMeasures.jl , 2018.[90] J. Rappaz. Approximation of the spectrum of a non-compact operator given by the magnetohydrodynamic stability of a plasma.
Numerische Mathematik , 28(1):15–24, 1977.
OUNDATIONS OF SPECTRAL COMPUTATIONS 53 [91] J. Rappaz, J. Sanchez Hubert, E. Sanchez Palencia, and D. Vassiliev. On spectral pollution in the finite element approximation ofthin elastic “membrane” shells.
Numerische Mathematik , 75(4):473–500, 1997.[92] M. Reed and B. Simon.
Methods of modern mathematical physics. II. Fourier analysis, self-adjointness . Academic Press [Har-court Brace Jovanovich, Publishers], New York-London, 1975.[93] E. Schr¨odinger. A method of determining quantum-mechanical eigenvalues and eigenfunctions.
Proc. Roy. Irish Acad. Sect. A. ,46:9–16, 1940.[94] J. Schwinger. Unitary operator bases.
Proc. Nat. Acad. Sci. U.S.A. , 46:570–579, 1960.[95] E. Shargorodsky. On the level sets of the resolvent norm of a linear operator.
Bull. Lond. Math. Soc. , 40(3):493–504, 2008.[96] E. Shargorodsky. On the limit behaviour of second order relative spectra of self-adjoint operators.
Journal of Spectral Theory ,3(4):535–552, 2013.[97] M. Shub and S. Smale. On the intractability of Hilbert’s Nullstellensatz and an algebraic version of NP (cid:54) = P ? Duke Mathemat-ical Journal , 81(1):47–54, 1995.[98] P. Siegl and D. Krejˇciˇr´ık. On the metric operator for the imaginary cubic oscillator.
Physical Review D , 86(12):121702, 2012.[99] B. Simon. Some quantum operators with discrete spectrum but classically continuous spectrum.
Annals of Physics , 146(1):209–220, 1983.[100] J. Sj¨ostrand and M. Zworski. Asymptotic distribution of resonances for convex obstacles.
Acta Mathematica , 183(2):191–253,1999.[101] S. Smale. The fundamental theorem of algebra and complexity theory.
American Mathematical Society. Bulletin. , 4(1):1–36,1981.[102] S. Smale. On the efficiency of algorithms of analysis.
Bull. Amer. Math. Soc. (N.S.) , 13(2):87–121, 1985.[103] S. Smale. Complexity theory and numerical analysis. In
Acta numerica, 1997 , volume 6 of
Acta Numer. , pages 523–551. Cam-bridge Univ. Press, Cambridge, 1997.[104] G. Szeg˝o. Beitr¨age zur Theorie der Toeplitzschen Formen.
Mathematische Zeitschrift , 6(3-4):167–202, 1920.[105] T. D. Tai. On the simpleness of zeros of Stokes multipliers.
Journal of Differential Equations , 223(2):351–366, 2006.[106] G. Teschl.
Jacobi operators and completely integrable nonlinear lattices , volume 72 of
Mathematical Surveys and Monographs .American Mathematical Society, Providence, RI, 2000.[107] A. V. Turbiner. Double well potential: perturbation theory, tunneling, WKB (beyond instantons).
International Journal of ModernPhysics A , 25(02n03):647–658, 2010.[108] A. M. Turing. On Computable Numbers, with an Application to the Entscheidungsproblem.
Proc. London Math. Soc. (2) ,42(3):230–265, 1936.[109] M. Webb and S. Olver. Spectra of Jacobi operators via connection coefficient matrices. arXiv:1702.03095 , 2017.[110] S. Weinberger.
Computers, Rigidity, and Moduli: The Large-Scale Fractal Geometry of Riemannian Moduli Space . PrincetonUniversity Press, USA, 2004.[111] E. J. Weniger. A convergent renormalized strong coupling perturbation expansion for the ground state energy of the quartic,sextic, and octic anharmonic oscillator.
Annals of Physics , 246(1):133–165, 1996.[112] H. Weyl.
The theory of groups and quantum mechanics . Dover Publications, Inc., New York, 1950.[113] Z. Zhang. How many numerical eigenvalues can we trust?
Journal of Scientific Computing , 65(2):455–466, 2015.[114] S. Zhao. On the spurious solutions in the high-order finite difference methods for eigenvalue problems.
Computer methods inapplied mechanics and engineering , 196(49-52):5031–5046, 2007.[115] M. Zworski. Resonances in physics and geometry.
Notices Amer. Math. Soc. , 46(3):319–328, 1999.[116] M. Zworski. Scattering resonances as viscosity limits. In
Algebraic and Analytic Microlocal Analysis , pages 635–654. Springer,2013. A PPENDIX
A. C
OMPUTATIONAL R OUTINES
We provide pseudocode for the algorithms of this paper, all of which provide sharp classifications in the SCIhierarchy.
Algorithm 1:
The routine
CompSpecUB computes spectra of unbounded operators on l ( N ) (or, more gen-erally, graphs) using the subroutines CompInvg and
DistSpec described above, and provides Σ errorcontrol. The subroutine IsPosDef checks whether a matrix is positive definite and is a standard routine thatcan be implemented in a myriad of ways. In practice, the while loop in
DistSpec is replaced by a muchmore efficient interval bisection method. An alternative method for sparse matrices (which, however, doesnot rigorously guarantee an error bound on the smallest singular values) is to compute the smallest singularvalues of the rectangular matrices using iterative methods.
Function
CompInvg( n , y , g ) Input : n ∈ N , y ∈ R + , g : R + → R + Output: m ∈ R + , an approximation to g − ( y ) m = min { k/n : k ∈ N , g ( k/n ) > y } endFunction DistSpec( A , n , z , f ( n ) ) Input : n ∈ N , f ( n ) ∈ N , matrix A , z ∈ C Output: y ∈ R + , an approximation to the function z (cid:55)→ (cid:107) R ( z, A ) (cid:107) − B = ( A − zI )(1 : f ( n ) , n ) C = ( A − zI ) ∗ (1 : f ( n ) , n ) S = B ∗ BT = C ∗ Cν = 1 , l = 0 while ν = 1 do l = l + 1 p = IsPosDef ( S − l n ) q = IsPosDef ( T − l n ) ν = min( p, q ) end y = ln endFunction CompSpecUB( A , n , { g m } , f ( n ) , c n ) Input : n ∈ N , f ( n ) ∈ N , c n ∈ R + (bound on dispersion), g m : R + → R + , A ∈ Ω g Output: Γ n ( A ) ⊂ C , an approximation to Sp( A ) , E n ( A ) ∈ R + , the error estimate G = Grid( n ) (see (6.2)) for z ∈ G do F ( z ) = DistSpec( A , n , z , f ( n ) ) if F ( z ) ≤ ( | z | + 1) − thenfor w j ∈ B CompInvg( n , F ( z ) , g (cid:100)| z |(cid:101) ) ( z ) ∩ G = { w , ..., w k } do F j = DistSpec( A , n , w j , f ( n ) ) end M z = { w j : F j = min q { F q }} else M z = ∅ endend Γ n ( A ) = ∪ z ∈ G M z E n ( A ) = max z ∈ Γ n ( A ) { CompInvg( n , DistSpec( A , n , z , f ( n ) ) + c n , g (cid:100)| z |(cid:101) ) } end OUNDATIONS OF SPECTRAL COMPUTATIONS 55
Algorithm 2:
PseudoSpecUB computes Γ n ( A ) ⊂ Sp (cid:15) ( A ) with lim n →∞ Γ n ( A ) = Sp (cid:15) ( A ) . Function
PseudoSpecUB( A , n , f ( n ) , c n , (cid:15) ) Input : n ∈ N , f ( n ) ∈ N , c n ∈ R + , A ∈ ˆΩ , (cid:15) > Output: Γ n ( A ) ⊂ C , an approximation to Sp (cid:15) ( A ) G = Grid( n ) for z ∈ G do F ( z ) = DistSpec( A , n , z , f ( n ) ) + c n end Γ n ( A ) = (cid:83) { z ∈ G | F ( z ) < (cid:15) } end Algorithm 3:
TestSpec solves { Ξ , ˆΩ × K ( C ) } (does K , compact, intersect Sp( A ) ) with input K n andaccess to γ n ( z, A ) (e.g. via DistSpec ). Similarly,
TestPseudoSpec solves { Ξ , ˆΩ × K ( C ) } (does K ,compact, intersect Sp (cid:15) ( A ) ) with input K n , (cid:15) > and access to γ n ( z, A ) . Function
TestSpec( n , n , K n , γ n ( z, A ) ) Input : n , n ∈ N , K n an approximation to K , access to evaluation of γ n ( z, A ) . Output: Γ n ,n ( A ) , an approximation of Ξ ( A ) . Γ n ,n ( A ) = “Does there exist some z ∈ K n such that γ n ( z, A ) < / n ?” endFunction TestPseudoSpec( n , n , K n , γ n ( z, A ) , (cid:15) ) Input : n , n ∈ N , K n an approximation to K , access to evaluation of γ n ( z, A ) , (cid:15) > . Output: Γ n ,n ( A ) , an approximation of Ξ ( A ) . Γ n ,n ( A ) = “Does there exist some z ∈ K n such that γ n ( z, A ) < / n + (cid:15) ?” end Algorithm 4:
SpecGap solves the spectral gap problem in Theorem 3.11, and requires an eigenvalue solverto implement Corollary 6.9 to compute all n eigenvalues of P n AP n to arbitrary precision. Function
SpecGap( n , n , P n AP n ) Input : n , n ∈ N , P n AP n the square truncation of the matrix A Output: Γ n ,n ( A ) , an approximation to Ξ gap ( A ) . if n = 1 then Set Γ n ,n ( A ) = 1 else for k ∈ { , ..., n } do Compute l k = µ k − µ k + (cid:15) k , | (cid:15) k | ≤ /k ,using Corollary 6.9 and notation of Lemmas 8.1, and 8.2 applied to P k AP k . end Set J n = [0 , / (2 n )] and J n = (1 /n , ∞ ) if { l k : k ∈ { , ..., n } ∩ ( J n ∪ J n ) } = ∅ then Set Γ n ,n ( A ) = 0 else Let ˜ k ≤ n be maximal with l ˜ k ∈ J n ∪ J n . if l ˜ k ∈ J n then Set Γ n ,n ( A ) = 0 else Set Γ n ,n ( A ) = 1 . endendendend Algorithm 5:
SpecClass solves spectral classification problem in Theorem 3.11. As well as an eigenvaluesolver to implement Corollary 6.9, we need the algorithm
CompSpecUB denoted by ˆΓ n , which computes thespectrum together with an error bound E ( n, · ) on the output. Function
SpecClass( n , n , A , f ) Input : n , n ∈ N , A ∈ (cid:101) Ω f SA , f the dispersion bounding function Output: Γ n ,n ( A ) , an approximation to Ξ class ( A ) . if n ≤ n then Set Γ n ,n ( A ) = 1 else for n ∈ { , ..., n } and j ∈ { , ..., n − } do Compute l jn = µ nj +1 − µ n + (cid:15) jn , (cid:12)(cid:12) (cid:15) jn (cid:12)(cid:12) ≤ /n ,using Corollary 6.9 and notation of Lemmas 8.1, and 8.2 applied to P n AP n . end Set J n = [0 , / (2 n )] and J n = (1 /n , ∞ ) for j ∈ { , ..., n } do Let k jn ,n be maximal with ≤ k jn ,n < n such that l jk jn ,n ∈ J n ∪ J n if such k jn ,n exists. endif k n ,n exists with l k n ,n ∈ J n then Set Γ n ,n ( A ) = 1 else if any k mn ,n exists with l mk mn ,n ∈ J n for ≤ m ≤ n then Set Γ n ,n ( A ) = 2 else for k ∈ { , ..., n } do Set a k ( A ) = min x ∈ ˆΓ k ( A ) { x + E ( k, x ) } ,Set q k = E ( k, a k ( A ) + 1 /n ) + 1 /k . end Set b n ,n ( A ) = min { q k : 1 ≤ k ≤ n } . if b n ,n ( A ) ≥ /n then Set Γ n ,n ( A ) = 3 else Set Γ n ,n ( A ) = 4 . endendendendend OUNDATIONS OF SPECTRAL COMPUTATIONS 57
Algorithm 6:
DiscreteSpec computes the closure of the discrete spectrum of A . The approximation ofthe essential spectrum, ˜Γ n ,n ( A ) , is described in the proof of convergence and was given in [8]. Moreover, lim n →∞ Γ n ,n ( A ) ⊂ Sp d ( A ) (see (8.3)), and converges up to Sp d ( A ) as n → ∞ . Given z n → z , Multiplicity computes the multiplicity, h ( A, z ) , of the eigenvalue z . Function
DiscreteSpec( n , n , ˆΓ n ( A ) , E ( n , · ) , ˜Γ n ,n ( A ) ) Input : n , n ∈ N , ˆΓ n ( A ) an approximation to Sp( A ) , error estimate E ( n , · ) over ˆΓ n ( A ) , ˜Γ n ,n ( A ) anapproximation to Sp ess ( A ) Output: Γ n ,n ( A ) , an approximation to Sp d ( A ) . if n ≤ n thenfor n ≤ k ≤ n do ζ k,n ( A ) = { z ∈ ˆΓ n ( A ) : E ( n , z ) < dist( z, ˜Γ k,n ( A )) − /k } for z, w ∈ ζ k,n ( A ) do z ∼ n w if and only if B E ( n ,w j ) ( w j ) ∩ B E ( n ,w j +1 ) ( w j +1 ) (cid:54) = ∅ for some z = w , w , ..., w n = w ∈ ζ k,n ( A ) end This gives equivalence classes [ z ] , ..., [ z m ] for j ∈ { , ..., m } do Choose z k j ∈ [ z j ] of minimal E ( n , · ) endif ∪ j ∈{ ,...,m } { z k j } (cid:54) = ∅ then Φ k,n ( A ) = ∪ j ∈{ ,...,m } { z k j } else Φ k,n ( A ) = ∅ . endendif At least one of Φ k,n ( A ) (cid:54) = ∅ then Γ n ,n ( A ) = Φ k,n ( A ) = ∅ with k minimal such that ζ k,n ( A ) (cid:54) = ∅ else Γ n ,n ( A ) = { } endelse Γ n ,n ( A ) = { } . endendFunction Multiplicity( A , n , n , f ( n ) , z n , d n ) Input : n , n ∈ N , f ( n ) ∈ N , A ∈ Ω d N , z n ∈ C , d n Output: h n ,n ( A, z n ) , an integer approximation to h ( A, z ) where z n → z . B = [( A − z n I )(1 : f ( n ) , n )] ∗ [( A − z n I )(1 : f ( n ) , n )] − (1 /n − d n ) I [ L, D ] = ldl ( B ) if D is diagonal then Find J the set of j with D ( j, j ) < h n ,n ( A, z n ) = | J | else Find J the set of j with size block D ( j, j ) < Find J the number of negative eigenvalues corresponding to size blocks by looking at trace anddeterminant h n ,n ( A, z n ) = | J | + | J | endend Algorithm 7:
DiscSpecEmpty computes Ξ d ( A ) (is the discrete spectrum non-empty) in two limits for A ∈ Ω f N , the class of bounded normal operators with known dispersion bounding function. The inputsare the algorithm ˆΓ n computing the spectrum (for example, CompSpecUB ), the error control E ( n, z ) (thatconverges to the true error uniformly on compact subsets of C ) and the height two tower ˜Γ n ,n presentedin [8] to compute the essential spectrum. Function
DiscSpecEmpty( n , n , ˆΓ n ( A ) , E ( n , · ) , ˜Γ n ,n ( A ) ) Input : n , n ∈ N , ˆΓ n ( A ) an approximation to Sp( A ) , error estimate E ( n , · ) over ˆΓ n ( A ) , ˜Γ n ,n ( A ) anapproximation to Sp ess ( A ) Output: Γ n ,n ( A ) , an approximation to Ξ d ( A ) . ζ n ,n ( A ) = { z ∈ ˆΓ n ( A ) : E ( n , z ) < dist( z, ˜Γ n ,n ( A )) − /n } if ζ n ,n ( A ) (cid:54) = ∅ then Γ n ,n ( A ) = 1 else Γ n ,n ( A ) = 0 . endend Algorithm 8:
ApproxEigenvector takes as input A , n , f ( n ) , z n and the bound E ( n, z n ) where σ ( P f ( n ) ( A − z n I ) | P n ( l ( N )) ) ≤ E ( n, z n ) . Given δ > , it computes an approximate eigenvector x n (offinite support) satisfying (cid:107) ( A − z n I ) x n (cid:107) ≤ (cid:107) x n (cid:107) ( E ( n, z n ) + c n + δ ) and − δ < (cid:107) x n (cid:107) < δ. Function
ApproxEigenvector( A , n , f ( n ) , z n , E ( n, z n ) , δ ) Input : n ∈ N , f ( n ) ∈ N , A , z n ∈ C , error bound E ( n, z n ) and tolerance δ > Output: x n ∈ C n , a vector satisfying (cid:107) ( A − z n I ) x n (cid:107) ≤ (cid:107) x n (cid:107) ( E ( n, z n ) + c n + δ ) (cid:15) = ( E ( n, z n ) + δ ) B = [( A − z n I )(1 : f ( n ) , n )] ∗ [( A − z n I )(1 : f ( n ) , n )] − (cid:15)I [ L, D ] = ldl ( B ) if D is diagonal then Find i with D ( i, i ) < y = e i else Find y eigenvector of D with eigenvalue < end Solve upper triangular system for x n with y = L ∗ x n then normalise to precision δ ..