Computing Spectral Measures and Spectral Types

Spectral measures arise in numerous applications such as quantum mechanics, signal processing, resonance phenomena, and fluid stability analysis. Similarly, spectral decompositions (into pure point, absolutely continuous and singular continuous parts) often characterise relevant physical properties such as the long-time dynamics of quantum systems. Despite new results on computing spectra, there remains no general method able to compute spectral measures or spectral decompositions of infinite-dimensional normal operators. Previous efforts have focused on specific examples where analytical formulae are available (or perturbations thereof) or on classes of operators that carry a lot of structure. Hence the general computational problem is predominantly open. We solve this problem by providing the first set of general algorithms that compute spectral measures and decompositions of a wide class of operators. Given a matrix representation of a self-adjoint or unitary operator, such that each column decays at infinity at a known asymptotic rate, we show how to compute spectral measures and decompositions. We discuss how these methods allow the computation of objects such as the functional calculus, and how they generalise to a large class of partial differential operators, allowing, for example, solutions to evolution PDEs such as the linear Schrödinger equation on L2(Rd)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2({\mathbb {R}}^d)$$\end{document}. Computational spectral problems in infinite dimensions have led to the Solvability Complexity Index (SCI) hierarchy, which classifies the difficulty of computational problems. We classify the computation of measures, measure decompositions, types of spectra, functional calculus, and Radon–Nikodym derivatives in the SCI hierarchy. The new algorithms are demonstrated to be efficient on examples taken from orthogonal polynomials on the real line and the unit circle (giving, for example, computational realisations of Favard’s theorem and Verblunsky’s theorem, respectively), and are applied to evolution equations on a two-dimensional quasicrystal.


Introduction
The analysis and computation of spectral properties of operators form core parts of many branches of science and mathematics, arising in diverse fields such as differential and integral equations, orthogonal polynomials, quantum mechanics, statistical mechanics, integrable systems and optics [15,43,44,69,110,129,138]. Correspondingly, the problem of numerically computing the spectrum, σ (T ), of an operator T acting on the canonical separable Hilbert space l 2 (N) has attracted a large amount of interest over the last 60 years or so [5][6][7][8]21,22,24,26,27,34,38,40,46,49,73,[79][80][81]94,98,99,116,120,[122][123][124]. However, the richness, beauty and difficulties that are encountered in infinite dimensions lie not just in the set σ (T ) ⊂ C, but also in the generalisation of projections onto eigenspaces and the possibility of different spectral types. Specifically, given a normal operator T , there is an associated projection-valued measure (resolution of the identity), which we denote by E T , whose existence is guaranteed by the spectral theorem and whose support is σ (T ) [86,87,113]. This allows the representation of the operator T as an integral over σ (T ), analogous to the finite-dimensional case of diagonalisation: (1.1) where D(T ) denotes the domain of T . For example, if T is compact, then E T corresponds to projections onto eigenspaces, familiar from the finite-dimensional setting. However, in general, the situation is much richer and more complicated, with different types of spectra (pure point, absolutely continuous and singular continuous). An excellent and readable introduction can be found in Halmos' article [77]. The computation of E T , along with its various decompositions and their supports, is of great interest, both theoretically and for practical applications. For example, spectral measures are intimately related to correlation functions in signal processing, resonance phenomena in scattering theory, and stability analysis for fluids. Moreover, the computation of E T allows one to compute many additional objects, which we provide the first general algorithms for in this paper, such as the functional calculus (Theorem 4.1), the Radon-Nikodym derivative of the absolutely continuous component of the measure (Theorem 4.2), and the spectral measures and spectral set decompositions (Theorem 3.2 and Theorem 5.1). For instance, in Sect. 1.2.1 we show how our results allow the computation of spectral measures and the functional calculus of almost arbitrary selfadjoint partial differential operators on L 2 (R d ). An important class of examples is given by solutions of evolution equations such as the Schrödinger equation on L 2 (R d ) with a potential of locally bounded total variation. We prove that this is computationally possible even when an algorithm is only allowed to point sample the potential. A numerical example of fractional diffusion for a discrete quasicrystal is also given in Sect. 6.4. Despite its importance, there has been no general method able to compute spectral measures of normal operators. Although there is a rich literature on the theory of spectral measures, most of the efforts to develop computational tools have focused on specific examples where analytical formulae are available (or perturbations thereof) or on classes of operators that carry a lot of structure. Indeed, apart from the work of Webb and Olver [151] (which deals with compact perturbations of tridiagonal Toeplitz operators) and methods for computing spectral density functions of Sturm-Liouville problems and other highly structured operators, there has been limited success in computing the measure E T . 1 In some sense, this is not surprising given the difficulty in rigorously computing spectra. One can consider computing spectral measures/projections as the infinite-dimensional analogue of computing projections onto eigenspaces. 2 Thus, from a numerical/computational point of view, the current state of affairs in infinite-dimensional spectral computations is highly deficient, analogous in finite dimensions to being able to compute the location of eigenvalues but not eigenvectors! It has been unknown if the general computation of spectral measures is possible, even for simple subclasses such as discrete Schrödinger operators. In other words, the computational problem of "diagonalisation" through computing spectral measures remains an important and predominantly open problem.
In this paper, we solve this problem by providing the first set of algorithms for the computation of spectral measures for a large class of self-adjoint and unitary operators, namely, those whose matrix columns decay at a known asymptotic rate. This is a very weak assumption and covers the majority of operators, even unbounded, found in applications. In particular, those whose representation is sparse (such as many of the graph or lattice operators typically dealt with in physics) and also partial differential operators, once a suitable basis has been selected (see Theorem 1.1 and Appendix B). We also show how to compute spectral measure decompositions, the functional calculus, the density of the absolutely continuous part of the measure (Radon-Nikodym derivative) and different types of spectra (pure point, absolutely continuous and singular continuous-these sets often characterise different physical properties in quantum mechanics [2,41,42,53,61,93,119,127]). A central ingredient of these new algorithms is the computation of the resolvent operator with error control through appropriate rectangular truncations (Theorem 2.1). Furthermore, we demonstrate the applicability of our algorithms. The algorithms are efficient and parallelisable, allowing large scale computations.

The solvability complexity index and classification of problems.
A surprise thrown up by the infinite-dimensional spectral problem, which turns out to be quite generic, is the Solvability Complexity Index (SCI) [81]. The SCI provides a hierarchy for classifying the difficulty of computational problems. In classical numerical analysis, when computing spectra, one hopes to construct an algorithm, n , with one limit such that for an operator T , n (T ) → σ (T ), as n → ∞, (1.2) preferably with some form of error control or rate of convergence. However, this is not always possible. For example, when considering the class of bounded operators, the best possible alternative is an algorithm depending on three indices n 1 , n 2 , n 3 such that Any algorithm with fewer than three limits will fail on some bounded operator, and neither error control nor convergence rates on any of the limits are possible since these would reduce the required number of limits. However, for self-adjoint operators, it is possible to reduce the number of limits to two, but not one [12,81]. With more structure (such as sparsity or column decay) it is possible to compute the spectrum in one limit with a certain type of error control [40]. Hence, the only way to characterise the computational spectral problem is through a hierarchy, classifying the difficulty of computing spectral properties of different subclasses of operators. The SCI classifies difficulty by considering the minimum number of limits that one must take to calculate the quantity of interest (see Appendix A for a full definition). The SCI has roots in the work of Smale [132,134], and his program on the foundations of computational mathematics and scientific computing, though it is quite distinct. The notions of Turing computability [145] and computability in the Blum-Shub-Smale (BSS) [17] sense become special cases, and impossibility results that are proven in the SCI hierarchy hold in all models of computation. The phenomenon of needing several limits also covers general numerical analysis problems, such as Smale's question on the existence of purely iterative algorithms for polynomial root finding [133]. As demonstrated by McMullen [101,102] and Doyle and McMullen [51], this is a case where several limits are needed in the computation, and their results become special cases of classification in the SCI hierarchy [12]. Extensions of the hierarchy to error control [33,36,37] also have potential applications in the growing field of computer-assisted proofs, where one must perform a computation with absolute certainty. See, for example, the work of Fefferman and Seco on the Dirac-Schwinger conjecture [55,56] and Hales on Kepler's conjecture (Hilbert's 18th problem) [75], both of which implicitly provide classifications in the SCI hierarchy. As well as spectral problems [12,33,[36][37][38][39][40]81], the SCI hierarchy has recently been used to solve problems related to computing resonances [13,14], computing solutions of semigroups and evolution PDEs [10,35], and computational barriers for stable and accurate neural networks [4]. An informal definition of the SCI hierarchy is as follows, with a detailed summary contained in Appendix A. The SCI hierarchy is based on the concept of a computational problem. This is described by a function : → M that we want to compute, where is some domain, and (M, d) is a metric space. For example, we could take (T ) = σ (T ) (the spectrum) for some class of operators and M the collection of non-empty closed subsets of C equipped with the Attouch-Wets metric. The SCI of a computational problem is the smallest number of limits needed in order to compute the solution. For a given set of evaluation functions (the information our algorithm is allowed to read-in our case, 1 or 2 defined in (1.18)), class of objects (in our case, subclasses of operators acting on l 2 (N)) and model of computation α (in this paper general, G, or arithmetic, A) we define: (i) α 0 is the set of problems that can be computed in finite time, the SCI = 0. (ii) α 1 is the set of problems that can be computed using one limit (the SCI = 1) with control of the error, i.e. ∃ a sequence of algorithms { n } such that d( n (A), (A)) ≤ 2 −n , ∀A ∈ . (iii) α 2 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 ∈ . (iv) α 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 1 } such that The class 1 is, of course, a highly desired class; however, non-trivial spectral problems are higher up in the hierarchy. For example, the following classifications are known [12,81]: The compact spectral problem is in 2 \ 1 .
Here, the notation \ indicates the standard "setminus". Hence, the computational spectral problem becomes an infinite classification theory to characterise the above hierarchy. In order to do so, there will, necessarily, have to be many different types of algorithms. Characterising the hierarchy will yield a myriad of different approaches, as different structures on the various classes of operators will require specific algorithms.
This paper provides classifications of spectral problems associated with E T (such as decompositions of the measure and spectrum) in the SCI hierarchy, some of which can be computed in one limit. We provide algorithms for these problems, and one of the main tools used is the computation of the resolvent operator R(z, T ) := (T − z I ) −1 with error control (Theorem 2.1). This is possible through appropriate rectangular truncations of the infinite-dimensional operator. This approach differs from finite-dimensional techniques, which typically consider square truncations.

Remark 1. (Recursivity and independence of the model of computation)
The constructive inclusion results we provide hold for arithmetic algorithms and the impossibility results hold for general algorithms. We refer the reader to Appendix A for a detailed explanation. Put simply, this means that the algorithms constructed can be recursively implemented with inexact input and restrictions to arithmetic operations over the rationals (it is also straightforward to implement them using interval arithmetic), whereas the impossibility results hold in any model of computation (such as the Turing or BSS models).

Summary of the main results
1.2.1. Partial differential operators For N ∈ N, consider the operator formally defined on L 2 (R d ) by where throughout we use multi-index notation with |k| = max{|k 1 | , . . . , We will assume that the coefficients a k (x) are complex-valued measurable functions on R d and that L is self-adjoint. For dimension d and r > 0 consider the space Let PDE consist of all such L such that the following assumptions hold: (1) The set C ∞ 0 (R d ) of smooth, compactly supported functions forms a core of L. (2) For each of the functions a k (x), there exists a positive constant A k and an integer B k (both possibly unknown) such that almost everywhere on R d , that is, we have at most polynomial growth of the coefficients.
We consider the case where our algorithms can do the following: 1. Evaluate any coefficient a k (x) to a given precision at x ∈ Q d , where Q denotes the field of rationals, and output an approximation in Q + iQ.
2. For each n, evaluate a positive number b n (L) such that the sequence {b n (L)} n∈N satisfies In Appendix B, we prove (and state a more precise version of) the following theorem.
Theorem 1.1 (Spectral properties of self-adjoint partial differential operators can be computed). Given the above set-up, there exist sequences of arithmetic algorithms that compute spectral measures, the functional calculus, and Radon-Nikodym derivatives of the absolutely continuous part of the measure over the class PDE . In other words, these objects can be computed in one limit with SCI = 1.

Remark 2.
As noted in Appendix B, we can extend this result to computing the decompositions into pure point, absolutely continuous and singular continuous parts of measures and spectra (with SCI > 1).
The above properties characterising PDE are deliberately very weak, and hence the class PDE is very large. For example, Schrödinger operators L = − + V with polynomially bounded potentials of locally bounded total variation are a subclass of PDE . Hence, in this case, the theorem says that we can compute the spectral properties (measures, functional calculus etc.) of L by point sampling the potential V , if we have an asymptotic bound on the total variation of V over finite rectangles. In particular, we can solve the Schrödinger equation by computing exp(−it L)u 0 with guaranteed convergence. Applications of some of our results to semigroups (where error control can also be provided) are given in [35]. The proof of Theorem 1.1 can also be extended to other domains such as the half-line, and can be adapted to cope with other types of coefficients that are not of locally bounded total variation (for instance, Coulombic potentials for Dirac or Schrödinger operators). In order to prove Theorem 1.1, we prove theorems for operators on l 2 (N).

Operators on l 2 (N)
We consider self-adjoint operators given as an infinite matrix T whose columns decay at a known asymptotic rate: for a sequence α n ↓ 0 and function f : N → N, where P n denotes the orthogonal projection onto the linear span of the first n canonical basis vectors. This set-up is adapted in Appendix B to prove Theorem 1.1, where we make use of the following theorems proven for operators on l 2 (N).
We consider the problem of computing spectral measures. Specifically, for operators of the form (1.7) we develop the first algorithms and SCI classifications for: • Theorem 3.1: A 2 classification for the projection-valued spectral measure and, through taking inner products, the computation of the scalar spectral measures defined via This is done for open sets U and can be extended to other types of sets such as closed intervals or singletons (Theorem 5.2 shows that the problem / ∈ G 1 in certain cases). These scalar-valued spectral measures play an important role in, for example, quantum mechanics, where they correspond to the distribution of the physical observable modelled by a Hamiltonian T . • Theorem 3.2: A 3 \ G 2 classification for the decompositions of the projection-valued and scalar-valued spectral measures into absolutely continuous, singular continuous and pure point parts. These decompositions often characterise different physical properties in quantum mechanics [2,41,42,53,61,93,119,127] • Theorem 4.1: A 2 classification for the functional calculus of operators, i.e. the computation of F(T ) for suitable functions F. For brevity, we consider functions that are bounded and continuous on the spectrum of T . However, the proof makes clear that wider classes can be dealt with. In some cases, A 1 error control is possible, for instance, when considering the holomorphic functional calculus (see Sect. 6.4).
A key application of this result is the computation of solutions of many evolution PDEs, such as the Schrödinger equation through the function F(z) = exp(−i zt). • Theorem 4.2: A 2 classification for the Radon-Nikodym derivatives (densities) of the absolutely continuous parts of the scalar spectral measures with convergence in the L 1 sense on an open set. This requires a certain separation condition, without which our algorithm converges (Lebesgue) almost everywhere.
We also consider the computation of spectra as sets in the complex plane. Convergence is measured using the Hausdorff metric in the bounded case and using the Attouch-Wets metric in the unbounded case (i.e. uniform convergence on compact subsets of C). Specifically, we provide the first algorithms computing these quantities and prove in Theorem 5.1 that: • The absolutely continuous spectrum σ ac (T ) can be computed in two limits but not one limit ( A 3 \ G 2 classification). • The pure point spectrum σ pp (T ) can be computed in two limits but not one limit ( A 3 \ G 2 classification). • The singular continuous spectrum σ sc (T ) can be computed in three limits ( A 4 classification). If f (n) − n ≥ √ 2n + 1/2, then the computation cannot be done in two limits ( / ∈ G 3 classification). That is, if the local asymptotic bandwidth is allowed to grow sufficiently rapidly, three limits are needed, and this computational problem is exceedingly difficult. We do not know whether this growth condition on f can be dropped. However, without it, the problem still requires at least two limits ( / ∈ G 2 classification).
The main tool used to prove the above results is • Theorem 2.1 and Corollary 2.2: The action of the resolvent x → R(z, T )x can be computed with error control. This also opens up potential applications in computerassisted proofs.
We demonstrate that the "one-limit" algorithms constructed in this paper are implementable and efficient. These provide the first set of algorithms addressing these problems, and we have provided extensive numerical experiments in Sect. 6. This includes orthogonal polynomials on the real line and unit circle (where we also discuss acceleration through extrapolation), as well as fractional diffusion for a two-dimensional quasicrystal. Since writing the initial version of this paper, our algorithms have also been implemented in the software package SpecSolve of [39], which uses the machinery of high-order rational kernels to accelerate computation of the Radon-Nikodym derivative of regular enough measures. In the current paper, we also show that in the case where the measure is regular enough, a global collocation method and local Richardson extrapolation for the computation of the Radon-Nikodym derivative can both accelerate convergence.
Finally, some brief remarks are in order.
(i) The impossibility results hold in general, even when restricted to tridiagonal operators. Furthermore, many of the impossibility results hold for structured operators such as bounded discrete Schrödinger operators. Our results (constructive algorithms and impossibility results) also carry over to a large class of normal operators, including unitary operators or skew-adjoint operators, both of which are important in applications. For the sake of clarity, we have stuck to the self-adjoint case in the statement of theorems and proofs. Numerical examples for unitary operators are given in Sect. 6.3. (ii) The difficulty encountered when computing the singular continuous spectrum is partly due to the negative definition of the singular continuous part of a measure. It is the "leftover" part of the measure, the part that is not continuous with respect to Lebesgue measure and does not contain atoms. The challenge of studying σ sc analytically also reflects this difficulty-singular continuous spectra were once thought to be rather rare or exotic. However, they are quite generic; see, for example, [128] for a precise topological statement to this effect. (iii) One might at first expect computational results to be independent of the function f due to tridiagonalisation. However, the infinite-dimensional case is much more subtle than the finite-dimensional case. Using Householder transformations on a bounded sparse self-adjoint operator T leads to a tridiagonal operator, but, in general, this operator is T restricted to a strict subspace of l 2 (N). Part of the operator may be lost in the strong operator limit. Instead, one must consider a sum of possibly infinitely many tridiagonal operators (see [78,Ch. 2 & 8]). Hence some spectral problems may have different classifications for different f .

Open problems.
The results of this paper form part of a program [33,34,[36][37][38][39][40] on determining the foundations of computation (boundaries of what is possible) for infinitedimensional spectral problems. As such, there remain many interesting open problems related to this paper, such as: • Computation of spectral measures for more general normal operators: Proposition 2.4 demonstrates how the techniques of this paper can be generalised to certain classes of normal operators whose spectrum does not necessarily lie along a curve. However, it is not obvious how to extend the techniques to completely general normal operators. For example, the method of integrating the resolvent along a contour cannot be easily extended to interior points of the spectrum. • Brown measure for non-normal operators: There is an extension of spectral measures to certain non-normal operators known as the Brown measure [25,71,72].
Computing spectra of non-normal operators is generally harder (in a sense made precise by the SCI hierarchy) than for normal operators (issues such as numerical stability are also present, even in the finite-dimensional case). It would be interesting to see if this phenomenon is also present for computing the Brown measure. Some works on approximating the Brown measure from matrix elements include [7,11,81].

A motivating example.
As a motivating example, consider the case of a Jacobi operator with matrix where a j , b j ∈ R and a j > 0. An enormous amount of work exists on the study of these operators, and the correspondence between bounded Jacobi matrices and probability measures with compact support [45,143]. The entries in the matrix provide the coefficients in the recurrence relation for the associated orthonormal polynomials. To study the canonical measure μ J , one usually considers the principal resolvent function, which is defined on C\σ (J ) via and then takes z close to the real axis. The function G is also known in the differential equations and Schrödinger communities as the Weyl m-function [64,143] and one can develop the discrete analogue of what is known as Weyl-Titchmarsh-Kodaira theory for Sturm-Liouville operators. Going back to the work of Stieltjes [136] (see also [1,149]), there is a representation of G as a continued fraction: One can also approximate G via finite truncated matrices [143]. However, there are two major obstacles to overcome when using (1.9) and its variants as a means to compute measures. First of all, this representation of the principal resolvent function is structurally dependent. For example, (1.9) is valid for the restricted case of Jacobi operators and hence one is led to seeking different methods for different operators (such as tight-binding Hamiltonians on two-dimensional lattices which have a growing bandwidth when represented as operators on l 2 (N)). Second, this would seem to give the wrong classification of the difficulty of the problem in the SCI hierarchy, giving rise to a tower of algorithms with two limits. One first takes a truncation parameter n to infinity to compute G(z) for Im(z) > 0, and then a second limit as z approaches the real axis. One of the main messages of this paper is that both of these issues can be overcome. Measures can be computed in one limit via an algorithm n and for a large class of operators. The only restriction is a known asymptotic decay rate of the off-diagonal entries. As a by-product, we compute the m-function of such operators with error control. Specific cases where this can be written explicitly do exist, such as periodic Jacobi matrices or perturbations of Toeplitz operators [52] (see also Sect. 1.7). However, there has been no general method proposed to compute the resolvent with error control. This consideration is crucial to allow the computation of measures in one limit.
To see how we might compute the measure using the resolvent, consider the Poisson kernels for the half-plane and the unit disk, defined respectively by P H (x, y) = 1 π y x 2 + y 2 and P D (x, y) = 1 2π where (r, θ) denote the usual polar coordinates. Let T be a normal operator, then for z / ∈ σ (T ), we have from the functional calculus that For self-adjoint T , z = u + iv ∈ C\R (u, v ∈ R) and x ∈ l 2 (N), we define We change variables λ = exp(iθ) and, with an abuse of notation, write d E T (λ) = i exp(iθ)d E T (θ ). A simple calculation then gives (1.13) Returning to our example, we see that the computation of the resolvent with error control allows the computation of G(z) with error control through taking inner products. By considering G(z) − G(z), this allows the computation of the convolution of the measure μ J with the Poisson kernel P H . In other words, we can compute a smoothed version of the measure μ J with error control. We can then take the smoothing parameter to zero to recover the measure (one can show that these smoothed approximations converge weakly in the sense of measures). Fig. 1 demonstrates this for a typical example.

Functional analytic setup.
We consider the canonical 3 separable Hilbert space H = l 2 (N), the set of square summable sequences with canonical basis {e n } n∈N . Let C(l 2 (N)) be the set of closed densely defined linear operators T such that span{e n : n ∈ N} forms a core of T and T * . The spectrum of T ∈ C(l 2 (N)) will be denoted by σ (T ) and the point spectrum (the set of eigenvalues) by σ p (T ). The latter set is not always closed and This paper focusses on the subclass N ⊂ C(l 2 (N)) of normal operators, that is, operators for which D(T ) = D(T * ) and T x = T * x for all x ∈ D(T ). The subclasses ⊂ N of self-adjoint and unitary operators will be denoted by SA and U respectively. For T ∈ SA and T ∈ U , σ (T ) ⊂ R and σ (T ) ⊂ T respectively, where T denotes the unit circle. Given T ∈ N and a Borel set B, E T B will denote the projection E T (B). Given x, y ∈ l 2 (N), we can define a bounded (complex-valued) measure μ T x,y via the formula (1.14) Via the Lebesgue decomposition theorem [76], the spectral measure μ T x,y can be decomposed into three parts the absolutely continuous part of the measure (with respect to the Lebesgue measure), the singular continuous part (singular with respect to the Lebesgue measure and atomless) and the pure point part. When considering SA , we will consider Lebesgue measure on R and let the Radon-Nikodym derivative of μ T x,y,ac with respect to Lebesgue measure. Of course this can be extended to the unitary (and, more generally, normal) case. This naturally gives a decomposition of the Hilbert space H = l 2 (N). For I = ac, sc and pp, we let H I consist of vectors x whose measure μ T x,x is absolutely continuous, singular continuous and pure point respectively. This gives rise to the orthogonal (invariant subspace) decomposition whose associated projections will be denoted by P T ac , P T sc and P T pp respectively. These projections commute with T and the projections obtained through the projection-valued measure. Of particular interest is the spectrum of T restricted to each H I , which will be denoted by σ I (T ). These different sets and subspaces often, but not always, characterise different physical properties in quantum mechanics (such as the famous RAGE theorem [2,53,119]), where a system is modelled by some Hamiltonian T ∈ SA [41,42,61,93]. For example, pure point spectrum implies the absence of ballistic motion for many Schrödinger operators [127].
1.6. Algorithmic setup. Given an operator T ∈ C(l 2 (N)), we can view it as an infinite matrix through the inner products 4 t i j = T e j , e i . All of the algorithms constructed can also be adapted to operators on l 2 (Z), either through the use of a suitable re-ordering of the basis, or though considering truncations of matrices in two directions, which is useful numerically since it preserves bandwidth. To be precise about the information needed to compute spectral properties, we define two classes of evaluation functions as These can be understood as different sets of information our algorithms are allowed to access (see Appendix A for a precise meaning). All the results proven in this paper can be easily extended to the case of inexact input. This means replacing the evaluation functions by i, j,m : C(l 2 (N)) → Q + iQ such that | f (1) i, j,m (T ) − T e j , e i | ≤ 2 −m and | f (2) i, j,m (T ) − T * e j , T * e i | ≤ 2 −m . Hence, the existence results carry over to algorithms that are only allowed to perform arithmetic operations over Q. This could be useful for rigorous bounds using interval arithmetic and computer-assisted proofs (for those familiar with the term, our algorithms are Turing recursive), though for brevity, we stick to 1 and 2 throughout. For discrete operators, the above information is often given to us, for example, in tight-binding models in physics or as a discretisation of a PDE, and hence it is natural to seek to compute spectral properties from matrix values. The set 2 is motivated via variational problems. For partial differential operators, such information is often given through inner products with a suitable basis, and in this case, the inexact input model is needed due to approximating the integrals (see Appendix B). For the classes considered in this paper, the evaluation sets 1 and 2 are in general different, yet the classifications in the SCI remain the same.
We will be concerned operators whose matrix representation has a known asymptotic rate of column/off-diagonal decay. Namely, let f : N → N with f (n) > n and let α = {α n } n∈N , β = {β n } n∈N be null sequences 5 of non-negative real numbers. We then define for X = SA or X = U, where P n denotes the orthogonal projection onto span{e 1 , . . . , e n }. In using the O(·) notation, the hidden constant is allowed to depend on the operator T or the vector x. We will also use When discussing SA f,α,β and SA f,α we will use the notation f,α,β and f,α . The collection of vectors in l 2 (N) satisfying P n x − x = O(β n ) will be denoted by V β . Finally, when α n ≡ 0, we will abuse notation slightly in requiring the stronger condition that Thus f,0 is the class of self-adjoint operators whose matrix sparsity structure is captured by the function f . For example, if f (n) = n + 1, we recover the class of self-adjoint tridiagonal matrices, the most studied class of infinite-dimensional operators. When discussing classes that include vectors x, we extend i to include pointwise evaluations of the coefficients of x. Other additions are sometimes needed, such as data regarding open sets as inputs for computations of measures, but this will always be made clear. When considering the general case of f,α , the function f and sequence α can also be considered as inputs to the algorithm-in other words, the same algorithm works for each class.

Connections with previous work.
We have mentioned the literature on infinitedimensional spectral problems and the SCI hierarchy. Computationally, our point of view in this paper is closest to the work of Olver, Townsend and Webb on practical infinite-dimensional linear algebra [104][105][106][107]151]. Their work includes efficient codes, such as the infinite-dimensional QL (IQL) algorithm [150] (see also [38] for the IQR algorithm, which has its roots in the work of Deift, Li and Tomei [46]), as well as theoretical results. A PDE version of the FEAST algorithm based on contour integration of the resolvent has recently been proposed by Horning and Townsend in [84], which computes discrete spectra. The set of algorithms we provide can be considered as a new member within the growing family of infinite-dimensional techniques.
A similar, though different, object studied in the mathematical physics literature, particularly when considering random Schrödinger operators, is the density of states [29,89,91], which we mention for completeness and also to avoid potential confusion. This object is defined via the "thermodynamic limit", where instead of considering the infinitedimensional operator T , one considers finite truncations, say P n T P n , and the limit n → ∞ of the measure x j ∈σ (P n T P n ) δ x j /n. This measure is subtly different from the spectral measure of T (for instance, on the discrete spectrum). The density of states is an important quantity in quantum mechanics, and there is a large literature on its computation. We refer the reader to the excellent review article [95], which discusses the most common methods. The idea of using the resolvent to approximate the density of states of finite matrices can be found in the method of [82], which approximates the imaginary part of Trace [R(z, P n T P n )] for Im(z) > 0. Similarly, in the random matrix literature, the connection is made through the Stieltjes transformation (see, for example, [9]). There are three immediate differences between our algorithms and those that compute the density of states. First, we seek to deal with the full, infinite-dimensional, operator directly to compute the spectral measure (and not the limit of increasing system sizes). Second, the object we are computing contains more refined spectral information of the operator and does not involve an averaging procedure. The density of states does not capture the full spectral information, such as the contribution of eigenvalues in the discrete spectrum, whereas the projection-valued spectral measure does. Third, there is a subtlety regarding the limits as Im(z) goes to zero and the truncation parameter goes to infinity (a similar trade-off also occurs in random matrix literature when theoretically analysing the density of states-see [54]). In our case, appropriate rectangular truncations of the infinite-dimensional operator are required to compute the resolvent with error control (see Theorem 2.1). This approach differs from finite-dimensional techniques, which typically consider square truncations.
The study of spectral measures also has a rich history in the theory of orthogonal polynomials and quadrature rules for numerical integration going back to the work of Szegő [50,139], briefly touched upon in Sect. 1.4. In certain cases, one can recover a distribution function for the associated measure of the Jacobi operator as a limit of functions constructed using Gaussian quadrature [31,Ch. 2]. Our examples in Sect. 6.1 can be considered as a computational realisation of Favard's theorem.
There are several results in the literature considering the computation of spectral density functions for Sturm-Liouville problems. In the case of Sturm-Liouville problems, the spectral density function corresponds to the multiplicative version of the spectral theorem. This is subtly different from the measures we compute, which arise from the projection-valued measure version of the spectral theorem. A common approach to approximate spectral density functions associated with Sturm-Liouville operators on unbounded domains is to truncate the domain and use the Levitan-Levinson formula, as implemented in the software package SLEDGE [59,60,111]. This approach can be computationally expensive since the eigenvalues cluster as the domain size increases; often, hundreds of thousands of eigenvalues and eigenvectors need to be computed. More sophisticated methods avoiding domain truncation are considered for special cases in [57,58], and an application in plasma physics can be found in [152]. These make use of the additional structure present in Sturm-Liouville problems using results analogous to (1.9) in the continuous case. Our results, particularly Theorem 1.1, hold for partial differential operators much more complicated than Sturm-Liouville operators (see Appendix B).
Finally, we wish to highlight the work of Webb and Olver [151], which is of particular relevance to the present study. There the authors studied, through connection coefficients, Jacobi operators that arise as compact perturbations of Toeplitz operators.
Similar perturbations (where a stronger exponential decay of the perturbation is crucial for analyticity properties of the resolvent) were studied in the work of Bilman and Trogdon [16] in connection with the inverse scattering transform for the Toda lattice. (See also the work of Trogdon, Olver and Deconinck [144] for computations of spectral measures for inverse scattering for the KdV equation.) The results proven in [151] can be stated in terms of the SCI hierarchy: • If the perturbation is finite rank (and known), the computation of σ pp lies in G 1 , and the computation of the μ ac lies in G 0 (note that σ ac is known analytically). • If the perturbation is compact with a known rate of decay at infinity, then the computation of the full spectrum σ lies in G 1 . The current paper complements the work of [151] by; considering operators much more general than tridiagonal compact perturbations of Toeplitz operators (we deal with arbitrary self-adjoint operators and assume we know f such that (1.20) holds) and partial differential operators, allowing operators to be unbounded, building algorithms that are arithmetic and can cope with inexact input, and considering computation of a wider range of spectral information. At the price of this greater generality, the objects we study are generally not computable with error control (unless one has local regularity assumptions on the measure-see [33,Ch. 4]), and some lead to computational problems higher up in the SCI hierarchy, though still computationally useful as we shall demonstrate. Our methods are also entirely different and rely on estimating the resolvent operator with error control (Theorem 2.1).

1.8.
Organisation of the paper. The paper is organised as follows. In Sect. 2 we consider the computation of the resolvent with error control and generalisations of Stone's formula. The computation of measures, their various decompositions and projections are discussed in Sect. 3. We then discuss the functional calculus and density of measures in Sect. 4. The computation of the different types of spectra as sets in the complex plane is discussed in Sect. 5. We run extensive numerical tests in Sect. 6, where we also introduce a new collocation method for the computation of the Radon-Nikodym derivative. We find that increased rates of convergence can also be obtained through iterations of Richardson extrapolation. A summary of the SCI hierarchy can be found in Appendix A and a proof of Theorem 1.1 in Appendix B. Throughout, our theorems and proofs use the notation introduced in Sect. 1.5 and Sect. 1.6.

Preliminary Results
The algorithms built in this paper rely on the computation of the action of the resolvent operator R(z, T ) = (T − z) −1 for z / ∈ σ (T ) with (asymptotic) error control. Given this, one can compute the action of the projections E T S for a wide range of sets S (Theorem 3.1 and its generalisations), and hence the measures μ T x,y . This section discusses the computation of the resolvent with error control and generalisations of Stone's formula, which relate the resolvent to the projection-valued measures.

Approximating the resolvent operator.
The key theorem for computing the action of the resolvent operator is the following, where we use σ 1 to denote the injection modulus of an operator defined as The proof of the following theorem boils down to a careful computation of a least-squares solution of a rectangular linear system. Theorem 2.1. Let α = {α n } n∈N and β = {β n } n∈N be null sequences, and f : where, for notational convenience, we drop the (T, x) dependence in the notation for C 1 and C 2 . Then there exists a sequence of arithmetical algorithms each of which use the evaluation functions in 1 (defined in (1.18)), such that each vector n (T, x, z) has finite support with respect to the canonical basis for each n, and Moreover, the following error bound holds If a bound on C 1 and C 2 are known, this error bound can be computed to arbitrary accuracy using finitely many arithmetic operations and comparisons.
Suppose that n is large enough so that is a least-squares solution of the optimisation problem argmin y P f (n) (T − z I )P n y − x . The linear space span{e n : n ∈ N} forms a core of T and hence also of T − z I . It follows by invertibility of T − z I that given any > 0, there exists an m = m( ) and a y = y( ) with P m y = y such that It follows that for all n ≥ m, This implies that In particular, since α and β are null, this implies that n (T, x, z) is uniformly bounded in n. Since > 0 was arbitrary, we also see that n (T, x, z) converges to R(z, T )x. Define the matrices Given the evaluation functions in 1 , we can compute the entries of these matrices to any given accuracy and hence also to arbitrary accuracy in the operator norm using finitely many arithmetic operations and comparisons (using the error in the Frobenius norm to bound the error in the operator norm). Denote approximations of B n and C n by B n and C n respectively and assume that for null sequences {u n }, {v n }. Note that B −1 n can be computed using finitely many arithmetic operations and comparisons. So long as u n is small enough, the resolvent identity implies that By taking u n and v n smaller if necessary (so that the algorithm is adaptive and it is straightforward to bound the norm of a finite matrix from above), we can ensure that From Proposition A.7 and a simple search routine, we can also compute σ 1 (P n (T * − z I )P f (n) (T − z I )P n ) to arbitrary accuracy using finitely many arithmetic operations and comparisons. Suppose this is done to an accuracy 1/n 2 and denote the approximation via τ n . We then define It follows that n (T, x, z) can be computed using finitely many arithmetic operations and, for large n, so that n (T, x, z) converges to R(z, T )x. By construction, n (T, x, z) has finite support with respect to the canonical basis.
Furthermore, the following error bound holds (which also holds if τ n ≤ 1/n) This bound converges to 0 as n → ∞. If the C 1 and C 2 are known, the bound can be approximated to arbitrary accuracy using finitely many arithmetic operations and comparisons.

Remark 3.
If T corresponds to a choice of basis in a space of functions (for example when using a spectral method), there is often a link between the regularity of the functions x and the decay of the terms β n . The bound (2.1) can then often be adapted to include such asymptotics, and hence indicate how large n needs to be to gain a given approximation.
Of course, a vast literature exists on computing R(z, T ), especially for infinite matrices with structure (such as being banded) and we refer the reader to [70,96,112,121] for a small sample. Note that if T is banded with bandwidth m, then we can take f (n) = n + m and the above computation can be done in O(nm 2 ) operations [66]. The following corollary of Theorem 2.1 will be used repeatedly in the following proofs.

Corollary 2.2. There exists a sequence of arithmetic algorithms
with the following properties: 1. For all (T, x) ∈ f,α,β and z ∈ C\R, n (T, x, z) has finite support with respect to the canonical basis and converges to R(z, T )x in l 2 (N) as n → ∞.

For any
where k are the algorithms from the statement of Theorem 2.1 and m(n, T, x, z) is a subsequence diverging to infinity as n → ∞. Clearly statement (1) holds so we must show how to choose the sequence m(n, T, x, z) such that (2) holds (and hence our algorithms will be adaptive). From (2.1), it is enough to show that m = m(n, T, x, z) can be chosen such that The left-hand side can be approximated to arbitrary accuracy using finitely many arithmetic operations and comparisons. Hence by repeatedly computing approximations to within α n + β n , we can choose the minimal m such that these approximate bounds are at most 2(α n + β n ).

Stone's formula and Poisson kernels.
Here we briefly discuss Stone's famous formula [32,113,137], which relates the convolution of spectral measures with Poisson kernels to the pointwise action of the projection-valued measures associated with an operator T ∈ SA as ↓ 0 (see Sect. 1.4). Stone's formula can also be generalised to unitary operators and a much larger class of normal operators (see Proposition 2.4). We include a short (and standard) proof of Proposition 2.3 for the benefit of the reader.
The following boundary limits hold: Proof. To prove (i), we can apply Fubini's theorem to interchange the order of integration and arrive at is bounded and converges pointwise as ↓ 0 to χ (a,b) (λ)+χ {a,b} (λ)/2, where χ S denotes the indicator function of a set S. Part (i) now follows from the dominated convergence theorem.
To prove (ii), we apply Fubini's theorem again, now noting that dψ.
We can split the interval into small intervals of width ρ (where 0 < ρ < 1) around each point where cos(ψ) = 1, and a finite union of intervals on which 1 − cos(ψ) is positive, bounded away from 0. On these later intervals, the limit vanishes as ↓ 0. Hence by periodicity and considering odd and even parts, we are left with considering dψ, Explicit integration yields I 2 ( , ρ) = O( log( )) and hence the contribution vanishes in the limit. We also have This converges to π as ↓ 0. Considering the contributions of I 1 and I 2 in (2.2), we see that (2.2) converges pointwise as ↓ 0 to Since the integral is also bounded, part (ii) now follows from the dominated convergence theorem and change of variables.
This type of construction can be generalised to T ∈ N whose spectrum lies on a regular enough curve. However, it is much more straightforward in the general case to use the analytic properties of the resolvent. The next proposition does this and also holds for operators whose spectrum does not necessarily lie along a curve. Proposition 2.4 (Generalised Stone's formula). Let T ∈ N and γ be a rectifiable positively oriented Jordan curve with the following properties. The spectrum σ (T ) intersects γ at finitely many points z 1 , . . . , z m and in a neighbourhood of each of the z i , γ is formed of a line segment meeting σ (T ) only at z i , at which point σ (T ) has a local exterior cone condition with respect to γ (see Fig. 2). Let x ∈ l 2 (N). Then we can define the Cauchy principal value integral of the resolvent R(z, T )x along γ and have Proof. We will argue for the case m = 1, and the general case follows in exactly the same manner. Let > 0 be small so that in a neighbourhood of the −ball around z 1 , γ is given by a straight line. We then decompose γ into two disjoint parts where γ 2 denotes the line segment of γ at most away from z 1 (as shown in Fig. 2). We set We then consider the inner integral , consider the contour integral along γ in Fig. 2. We see that and hence f (z 1 ) = −iπ . We would like to apply the dominated convergence theorem. Clearly, away from z 1 , f is bounded as ↓ 0. Now let 0 < δ < then By rotating and translating, we can assume that w = 1 and z 1 = 0 without loss of generality. Let λ 1 = Re(λ) and λ 2 = Im(λ). Using the cone condition gives α |λ 1 | ≤ |λ 2 | for some α > 0. Assume λ 1 = 0 then where x = /λ 1 and y = λ 2 /λ 1 . Note that y 2 ≥ α 2 and without loss of generality we take y ≥ α. Define Note that h(x, y) → 0 as |x| 2 + |y| 2 → ∞. We must show that h(x, y) is bounded above −1 for y ≥ α. It is enough to consider points where ∂h/∂ x = 0 which occur when x ± = ± 1 + y 2 . We have and hence we have proved the required boundedness. We then define The relation (2.3) now follows from the dominated convergence theorem.

Computation of Measures
For the sake of brevity, the analysis in the rest of this paper will consider the selfadjoint case T ∈ SA , which is the case most encountered in applications. However, the algorithms we build are based on

Full spectral measure.
We start by considering the computation of E T U x, where U ⊂ R is a non-trivial open set. In other words, U is not the whole of R or the empty set. The collection of these subsets will be denoted by U. To be precise, we assume that we have access to a finite or countable collection With an abuse of notation, we add this information as evaluation functions to i (defined in (1.18)). Then In other words, we can construct a convergent sequence of arithmetic algorithms for the problem. Remark 4. Essentially, this theorem tells us that if we can compute the action of the resolvent operator with asymptotic error control near the real axis, then we can compute the spectral measures of open sets in one limit. In the unitary case, this can easily be extended to relatively open sets of T if we can evaluate the resolvent near the unit circle. For any U ∈ U, the approximation of E T U x has finite support, and hence we can take inner products to compute μ T x,y (U ). Remark 5. One may wonder whether it is possible to upgrade the convergence of the algorithm in Theorem 3.1 from 2 to 1 . In other words, whether it is possible to compute the measure with error control. However, this is difficult because the measure may be singular. Theorem 5.2 shows this is impossible even for singleton sets and discrete Schrödinger operators acting on l 2 (N).
Proof of Theorem 3.1. Let T ∈ SA and z 1 , z 2 ∈ C\R. By the resolvent identity and self-adjointness of T , Hence, for z = u+i with > 0, the vector-valued function K H (u+i ; T, x) (considered with argument u) is Lipschitz continuous with Lipschitz constant bounded by −2 x /π . Now consider the class f,α,β × U and let (T, x, U ) ∈ f,α,β × U. From Corollary 2.2, we can construct a sequence of arithmetic algorithms, n , such that for all (T, x) ∈ f,α,β . It follows from standard quadrature rules and taking subsequences if necessary (using that {α n } and can be approximated to an accuracy C(T, x)/n using finitely many arithmetic operations and comparisons and the relevant set of evaluation functions 1 (the constant C now becomes C due to not knowing the exact value of x ). Recall that we assumed the disjoint union where a m , b m ∈ R∪{±∞} and the union is at most countable. Without loss of generality, we assume that the union is over m ∈ N. We then let a m,n , b m,n ∈ Q be such that a m,n ↓ a m and b m,n ↑ b m as n → ∞ with a m,n < b m,n and hence (a m,n , b m,n ) ⊂ (a m , b m ). Let then the proof of Stone's formula in Proposition 2.3 (essentially an application of the dominated convergence theorem) can be easily adapted to show that Note that we do not have to worry about contributions from endpoints of the intervals (a m , b m ) since we approximate strictly from within the open set U . To finish the proof, we simply let n (T, x, U ) be an approximation of the integral with accuracy C(T, x)/n. By the above remarks, such an approximation can be computed using finitely many arithmetic operations and comparisons from the relevant set of evaluation functions 1 .
This theorem can clearly be extended to cover the more general case of Proposition 2.4 if γ is regular enough to allow approximation of given the ability to compute R(z, T )x with asymptotic error control. Note that when it comes to numerically computing the integrals in Propositions 2.3 and 2.4, it is advantageous to deform the contour so that most of the contour lies far from the spectrum so that the resolvent has a smaller Lipschitz constant. The proof can also be adapted to compute E I x, where I = [a, b] is a closed interval, by considering intervals shrinking to [a, b] (a, b finite). A special case of this is the computation of the spectral measure of singleton sets. However, for these it much easier to directly use the formulae for T ∈ SA and T ∈ U respectively.

Measure decompositions and projections.
Recall from Sect. 1.5 that P T I denotes the orthogonal projection onto the space H T I , where I denotes a generic type (ac, sc, pp, c or s). We have included the continuous and singular parts denoted by c or s which correspond to H ac ⊕ H sc and H sc ⊕ H pp respectively. These are often encountered in mathematical physics. As in Sect. 3.1, we assume the decomposition in (3.1) and add the {a m , b m } as evaluation functions to i (defined in (1.18)). In this section, we prove the following theorem.
for I = ac, sc, pp, c or s. Then for i = 1, 2 To prove this theorem, it is enough, by the polarisation identity, to consider x = y (note that all the projections commute). We will split the proof into two parts: the A 3 inclusion, for which it is enough to consider 1 , and the G 2 exclusion, for which it is enough to consider 2 .

Proof of inclusion in Theorem 3.2
Proof of inclusion in Theorem 3.2. Since P T pp = I − P T c , P T ac = I − P T s and P T sc = P T s − P T pp , it is enough, by Theorem 3.1 and Remark 4, to consider only I = c and I = s.
Step 1: We first deal with I = c, where we shall use a similar argument to the proof of Theorem 4.1 (which is more general than what we need). We recall the RAGE theorem [2,53,119] as follows. Let Q n denote the orthogonal projection onto vectors in l 2 (N) with support outside the subset {1, . . . , n} ⊂ N. Then for any x ∈ l 2 (N), The proof of Theorem 4.1 is easily adapted to show that there exists arithmetic algorithms n,m using 1 such that Note that this bound can be made independent of s (as we have written above) by sufficiently approximating the function λ → exp(−iλs)χ U (λ) (it has known total variation for a given s and uniform bound). We now define Using the fact that for a, b ∈ l 2 (N), it follows that where g n (s) = Q n e −i T s χ U (T )x 2 . Clearly the first term converges to 0 as m → ∞, so we only need to consider the second. Using (3.4), it follows that for any > 0 But e −i T − I converges strongly to 0 as ↓ 0 and hence the quantity Step 2: Next we deal with the case I = s. Note that for z ∈ C\R, R(z, T )x, x is simply the Stieltjes transform (also called the Borel transform) of the positive measure μ T The Hilbert transform of μ T x,x is given by the limit with the limit existing (Lebesgue) almost everywhere. This object was studied in [108,109], where we shall use the result (since the measure is positive) that for any bounded continuous function f , 6 where a m , b m ∈ R∪{±∞} and the disjoint union is at most countable as in (3.1). Without loss of generality, we assume that the union is over m ∈ N. Due to the possibility of point spectra at the endpoints a m , b m , we cannot simply replace f by χ U in the above limit (3.5). However, this can be overcome in the following manner. Let ∂U denote the boundary of U defined by U \U and let ν denote the measure μ T x,x | ∂U . Let { f l } l∈N denote a pointwise increasing sequence of continuous functions, converging everywhere up to χ U , such that the support of each f l is contained in Such a sequence exists (and can easily be explicitly constructed) precisely because U is open. We first claim that To see this note that for any k ∈ N, the following inequalities hold where the last equality is due to (3.5). Taking k → ∞ yields so we are left with proving a similar bound for the limit supremum. Note that any point in the support of f l is of distance at least 1/ √ l from ∂U . It follows that there exists a constant C independent of t such that for any t ∈ supp( f l ), Now we let f ↓ χ U , with pointwise convergence everywhere. This is possible since the complement of U is open. By the dominated convergence theorem, and since was arbitrary, this yields lim sup where the last equality follows from the definition of ν. The claim (3.6) now follows. Let χ n be a sequence of non-negative continuous piecewise affine functions on R, bounded by 1 and such that χ n (t) = 0 if t ≤ n − 1 and χ n (t) = 1 if t ≥ n + 1. Consider the integrals where F m (t) is an approximation of The dominated convergence theorem implies that Note that continuity of χ n is needed to gain convergence almost everywhere and prevent possible oscillations about the level set {H μ T x,x (t) = n}. We also have The same arguments used to prove (3.6), therefore show that Hence, completing the proof of inclusion in Theorem 3.2.

Proof of exclusion in Theorem 3.2
To prove the exclusion, we need two results which will also be used in Sect. 5. First, we consider a result connected to Anderson localisation (Theorem 3.3) and, second, we consider a result concerning sparse potentials of discrete Schrödinger operators (Theorem 3.4). The free Hamiltonian H 0 acts on l 2 (N) via the tridiagonal matrix representation We define a Schrödinger operator acting on l 2 (N) to be an operator of the form where v is a bounded (real-valued) multiplication operator with matrix diag(v (1), v(2), . . .). Since Anderson's introduction of his famous model 60 years ago [3], there has been a considerable amount of work by both physicists and mathematicians aiming to understand the suppression of electron transport due to disorder (Anderson localisation).
A full discussion of Anderson localisation is beyond the scope of this paper, and we refer the reader to [29,42,90] for broader surveys. When considering Anderson localisation, we will assume that v = v ω = {v (1), v(2), . . .} is a collection of independent identically distributed random variables. Following [67], we assume that the single-site probability distribution has a density ρ ∈ L 1 (R) ∩ L ∞ (R) with ρ 1 = 1 (with respect to the standard Lebesgue measure). For such a potential, a measure of disorder is given by the quantity ρ −1 ∞ . The first result we need is the following theorem which follows straightforwardly from the technique of [67], and hence we have not provided a proof which would be almost verbatim to [67]. [67]). There exists a constant C > 0 such that if ρ ∞ ≤ C and ρ has compact support, then the operator H v ω + A has only pure point spectrum with probability 1 for any fixed self-adjoint operator A of the form

Theorem 3.3 (Anderson Localisation for Perturbed Operator
In other words, the operator A's matrix with respect to the canonical basis has only finitely many non-zeros. The second result we need is the following from [92]. Proof of exclusion in Theorem 3.2. Since P T pp = I − P T c , P T ac = I − P T s and P T sc = P T s − P T pp , it is enough, by Theorem 3.1 and Remark 4, to consider I = pp, ac and sc. We restrict the proof to considering bounded Schrödinger operators H v acting on l 2 (N), which are clearly a subclass of f,0 for f (n) = n + 1. Note that since the evaluation functions in 2 can be recovered from those in 1 in this special case, we can assume that we are dealing with 1 . We also set x = e 1 , with the crucial properties that this vector is cyclic and hence μ H v e 1 ,e 1 has the same support as σ (H v ), and that x ∈ V 0 . Throughout, we also take U = (0, 4).
Step 1: We begin with P T pp . Suppose for a contradiction that there does exist a sequence of general algorithms n such that We take a general algorithm, denoted n , from Theorem 3.1 which has e 1 ((0, 4)).
Since e 1 is cyclic, this limit is non-zero if (0, 4) ∩ σ (H v ) = ∅. We therefore define n (H v ) otherwise. We will use Theorem 3.3 and the following well-known facts: The strategy will be to construct a potential v such that (0, 4) ⊂ σ (H v ), yet n (H v ) does not converge. This is a contradiction since by our assumptions, for such a v we must have e 1 ((0, 4)) .
To do this, choose ρ = χ [−c,c] /(2c) for some constant c such that the conditions of Theorem 3.3 hold and define the potential v inductively as follows.
Let v 1 be a potential of the form v ω (with the density ρ) such that σ (H v 1 ) is pure point. Such a v 1 exists by Theorem 3.3 and we have P e 1 ((0, 4)). Hence for large enough n it must hold that n (H v 1 ) > 3/4. Fix n = n 1 such that this holds. Then n 1 (H v 1 ) only depends on {v 1 ( j) : j ≤ N 1 } for some integer N 1 by (i) of Definition A.2. Define the potential v 2 by v 2 ( j) = v 1 ( j) for all j ≤ N 1 and v 2 ( j) = 0 otherwise. Then by fact (2) above, P e 1 ((0, 4)) = 0, and hence n (H v 2 ) < 1/4 for large n, say for n = n 2 > n 1 . But then n 2 (H v 2 ) only depends on {v 2 ( j) : j ≤ N 2 } for some integer N 2 .
We repeat this process inductively switching between potentials which induce n k (H v k ) < 1/4 for k even and potentials which induce n k (H v k ) > 3/4 for k odd. Explicitly, if k is even then define a potential v k+1 by v k+1 ( j) = v k ( j) for all j ≤ N k and v k+1 ( j) = v ω ( j) (with the density ρ) otherwise such that the spectrum of H v k is pure point. Such a ω exists from Theorem 3.3 applied with the perturbation A to match the potential for j ≤ N k . If k is odd then we define v k+1 by v k+1 ( j) = v k ( j) for all j ≤ N k and v k+1 ( j) = 0 otherwise. We can then choose n k+1 such that the above inequalities hold and N k+1 such that n k+1 (H v k+1 ) only depends on {v k+1 ( j) : j ≤ N k+1 }. We also ensure that N k+1 ≥ N k + k.
Finally set v( j) = v k ( j) for j ≤ N k . It is clear from (iii) of Definition A.2, that n k (H v ) = n k (H v k ) and this implies that n k (H v ) cannot converge. However, since N k+1 ≥ N k + k, for any k odd we have v(N k + 1) = v(N k + 2) = · · · = v(N k + k) = 0. e 1 ((0, 4)) = 0 and therefore n (H v ) converges. This provides the required contradiction.
Step 2: Next we deal with I = ac. To prove that one limit will not suffice, our strategy will be to reduce a certain decision problem to the computation of ac . Let (M , d ) be the discrete space {0, 1}, let denote the collection of all infinite sequence {a j } j∈N with entries a j ∈ {0, 1} and consider the problem function ({a j }) : 'Does {a j } have infinitely many non-zero entries?', which maps to (M , d ). In Appendix A, it is shown that SCI( , ) G = 2 (where the evaluation functions consist of component-wise evaluation of the array {a j }). Suppose for a contradiction that n is a height one tower of general algorithms such that We will gain a contradiction by using the supposed tower to solve { , }.

Then by Theorem 3.4, P H
Since e 1 is cyclic, this limit is non-zero for H v , where v is of the form (3.10). We therefore define The above comments show that each of these is a general algorithm and it is clear that it converges to ({a j }) as n → ∞, the required contradiction.
Step 3: Finally, we must deal with I = sc. The argument is the same as

Two Important Applications
Two important applications of our techniques are the computation of the functional calculus and of the Radon-Nikodym derivative of μ T x,y,ac with respect to Lebesgue measure, denoted by ρ T x,y . Both of these have applications throughout mathematics and the physical sciences, some of which are explored numerically in Sect. 6. For example, suppose that we wish to solve the Schrödinger equation where H is some self-adjoint Hamiltonian. We can express the solution at time t as For example, the quantity known as the autocorrelation function [141], is simply the Fourier transform of the spectral measure dμ H u 0 ,u 0 . In particular, if dμ H u 0 ,u 0 is absolutely continuous, then ρ H u 0 ,u 0 and L form a Fourier transform pair. The computation of evolution generated by an operator is in some sense dual to the computation of the spectral measure. This interpretation of a time evolution can be adapted to describe many signals generated by PDEs [48,85,126] and stochastic processes [65,88] [118,Ch. 7]. In this section, we show how to compute the functional calculus and ρ T x,y .

F(T ) is a densely defined closed normal operator with dense domain given by
For simplicity, we will only deal with the case that F is a bounded continuous function on R, that is, F ∈ C b (R). In this case D(F(T )) is the whole of l 2 (N) (the variations |μ T x,y | are finite) and we can use standard properties of the Poisson kernel. We assume that given F ∈ C b (R), we have access to piecewise constant functions F n supported in [−n, n] such that F − F n L ∞ ([−n,n]) ≤ n −1 . Clearly other suitable data also suffices and, as usual, we abuse notation slightly by adding this information to the evaluation sets i (recall that i are defined in (1.18)). Then The inner integral is bounded since F is bounded and the Poisson kernel integrates to 1 along the real line. It also converges to F(λ) everywhere. Hence by the dominated convergence theorem We now use the same arguments used to prove Theorem 3.1. Using Corollary 2.2, together with K H (u + i/n; T, x) L ∞ (R) ≤ nC 1 and the fact that K H (u + i/n; T, x) is Lipschitz continuous with Lipschitz constant n 2 C 2 for some (possibly unknown) constants C 1 and C 2 , we can approximate this integral with an error that vanishes in the limit n → ∞.
If σ (T ) is bounded, then, with slightly more information available to our algorithms, a simpler proof holds using the Stone-Weierstrass theorem. Suppose that given x, the vectors T n x can be computed to arbitrary precision. There exists a sequence of polynomials p m (z) converging uniformly to F(z) on σ (T ). Assuming such a sequence can be explicitly constructed (for example using Bernstein or Chebyshev polynomials), we can take p m (T )x as approximations of F(T )x. If we can bound p m (z) − F(z) σ (T ) ≤ m with m null, then the vector F(T )x can be computed with error control. However, computing T n x for large n (even if x = e 1 ) may be computationally expensive as was found in the example in Sect. 6.4. We will also see in Sect. 6

.4 that if σ (T ) is bounded and F is analytic in an open neighbourhood of σ (T ), then F(T )x can be computed
with error control by deforming the integration contour away from the spectrum. Such a deformation is useful since the resolvent does not blow up along such a contour, and we can bound its Lipschitz constant.

Computation of the Radon-Nikodym derivative.
Recall the definition of the Radon-Nikodym derivative in (1.16) and note that ρ T x,y ∈ L 1 (R) for T ∈ SA . We consider its computation in the L 1 sense in the following theorem, where, as before, we assume (3.1), adding the approximations of U to our evaluation set 1 (defined in (1.18)), along with component-wise evaluations of a given vector y. However, we must consider the computation away from the singular part of the spectrum.
We restrict this map to the quadruples (T, x, y, U ) such that U is strictly separated from supp(μ T x,y,sc )∪supp(μ T x,y,pp ) and denote this subclass by f,α,β . Then { RN , f,α,β , 1 } ∈ A 2 . Furthermore, each output n (T, x, y, U ) of the algorithms constructed in the proof consists of a piecewise affine function, supported in U with rational knots and taking (complex) rational values at these knots. Remark 6. Essentially, this theorem tells us that if we can compute the action of the resolvent operator with asymptotic error control, then we can compute the Radon-Nikodym derivative of the absolutely continuous part of the measures on open sets which are a positive distance away from the singular support of the measure. The assumption that U is separated from supp(μ T x,y,sc )∪supp(μ T x,y,pp ) may seem unnatural but is needed to gain L 1 convergence of the approximation. However, without it, the proof still gives almost everywhere pointwise convergence.
Proof. Let (T, x, y, U ) ∈ f,α,β . For u ∈ U we decompose as follows (4.1) The first term converges to ρ T x,y | U in L 1 (U ) as ↓ 0 since ρ T x,y | U ∈ L 1 (U ). Since we assumed that U is separated from supp(μ T x,y,sc ) ∪ supp(μ T x,y,pp ), it follows that the second term of (4.1) converges to 0 in L 1 (U ) as ↓ 0. Hence we are done if we can approximate K H (u + i/n; T, x), y in L 1 (U ) with an error converging to zero as n → ∞.
Recall that K H (u + i/n; T, x) is Lipschitz continuous with Lipschitz constant at most n 2 x /π . By assumption, and using Corollary 2.2, we can approximate K H (u + i/n; T, x) to asymptotic precision with vectors of finite support. Hence the inner product f n (u) := K H (u + i/n; T, x), y can be approximated to asymptotic precision (now with a possibly unknown constant also depending on y ) and f n is Lipschitz continuous with Lipshitz constant at most n 2 x y /π .
Recall that U can be written as the disjoint union where a m , b m ∈ R ∪ {±∞} and the union is at most countable. Without loss of generality, we assume that the union is over m ∈ N. Given an interval (a m , b m ), let a m < z m,1,n < z m,2,n < · · · < z m,r m ,n < b m such that z m, j,n ∈ Q and z m, j,n − z m, j+1,n ≤ Here C is some unknown constant which occurs from the asymptotic approximation of f n that arises from Corollary 2.2 and we can always compute such f m,n in finitely many arithmetic operations and comparisons. Let n (T, x, y, U ) be the function that agrees with f m,n on (a m , b m ) for m ≤ n and is zero elsewhere. Clearly the nodes of n (T, x, y, U ) can be computed using finitely many arithmetic operations and comparisons and the relevant set of evaluation functions 1

. A simple application of the triangle inequality implies that
where the last term arises due to the piecewise affine interpolant. The bound clearly converges to zero as required.

Computing Spectra as Sets
We now turn to computing the different types of spectra as sets in the complex plane. Specifically, define the problem functions C I (T ) = σ I (T ) for I = ac, sc or pp. Note also that σ pp (T ) = σ p (T ), the closure of the set of eigenvalues. Recalling the definition of a computational problem in Appendix A, we compute these quantities in a metric space M with metric d M . Since we wish to include unbounded operators, we use the Attouch-Wets metric defined by for C 1 , C 2 ∈ Cl(C), where Cl(C) denotes the set of closed non-empty subsets of C. When considering bounded T , whose spectrum is necessarily a compact subset of C, where d(x, y) = |x − y| is the usual Euclidean distance. Note that for compact sets (and hence for bounded operators), the topological notions of convergence according to d H and d AW coincide. To allow the possibility that the different spectral sets σ I (T ) are empty, we add the empty set to our metric space as an isolated point (the space remains metrisable). 7 The main theorem of this section is the following: If f (n) − n ≥ √ 2n + 1 2 , then { C sc , f,0 , i } ∈ G 3 also holds.
In order to prove Theorem 5.1, we only need to prove the lower bounds for 2 and the upper bounds for 1 (recall that i are defined in (1.18)). These results show that despite the results of Sects. 2-4, in general it is very hard to compute the decomposition of the spectrum in the sense of (1.17). We also answer the question posed in Sect. 2.2 and prove that the spectral measures, while computable in one limit, cannot, in general, be computed with error control (see Theorem 5.2), unless one has additional regularity assumptions such as in [39] (computation with error control is made precise in [33,Ch. 4]).

Proof for point spectra. Proof that
To prove this, it is enough to consider bounded Schrödinger operators acting on l 2 (N) (defined in Sect. 3.2.2), which are clearly a subclass of f,0 for f (n) = n + 1. Note that since the evaluation functions in 2 can be recovered from those in 1 in this special case, we can assume that we are dealing with 1 . Suppose for a contradiction that there does exist a sequence of general algorithms, n , with We will construct a potential v such that n (H v ) does not converge. To do this, choose ρ = χ [−c,c] /(2c) for some constant c such that the conditions of Theorem 3.3 hold. We will use Theorem 3.3 and the following well-known facts: We will define the potential v inductively as follows. Let v 1 be a potential of the form v ω (with density ρ) such that [−c, 4 + c] ⊂ σ (H v 1 ) and σ (H v 1 ) is pure point. Such a v 1 exists by Theorem 3.3 and fact (2) above. Then for large enough n there exists z n ∈ n (H v 1 ) such that |z n − 2| ≤ 1. Fix n = n 1 such that this holds. Then n 1 (H v 1 ) only depends on {v 1 ( j) : j ≤ N 1 } for some integer N 1 by (i) of Definition A.2. Define the potential v 2 by v 2 ( j) = v 1 ( j) for all j ≤ N 1 and v 2 ( j) = 0 otherwise. Then by fact (1) above, n (H v 2 ) ∩ [1/2, 7/2] = ∅ for large n, say for n 2 . But then n 2 (H v 2 ) only depends on {v 2 ( j) : j ≤ N 2 } for some integer N 2 .
We repeat this process inductively switching between potentials which induce n k (H v k )∩[1/2, 7/2] = ∅ for k even and potentials which induce (H v k+1 ) and σ (H v k+1 ) is pure point. Such a ω exists from Theorem 3.3 and fact (2) above applied with the perturbation A to match the potential for j ≤ N k . If k is odd then we define v k+1 by v k+1 ( j) = v k ( j) for all j ≤ N k and v k+1 ( j) = 0 otherwise. We can then choose n k+1 such that the above intersections hold and N k+1 such that n k+1 ( (H v k ). But then this implies that n k (H v ) cannot converge, the required contradiction. A similar argument gives the following theorem, where V is used to denote the set of bounded real-valued potentials on N and 3 denotes the pointwise evaluations of such potentials.

Theorem 5.2 (Impossibility of computing spectral measures with error control). Consider the problem function
In other words, can be computed in one limit, but it cannot be computed with error control.
Proof. The result { , V×N, 3 } ∈ A 2 follows directly from the remarks after Theorem 3.1 and Theorem 2.1. Suppose for a contradiction that { , V × N, 3 } ∈ G 1 and that n is a sequence of general algorithms solving the problem with error control. It follows that for each j ∈ N, there exists a sequence of general algorithms It is clear that this is a general algorithm using 3 . Furthermore, with convergence from below. Now we may choose a v such that 1 ∈ σ p (H v ) (this can be achieved for example by taking a potential which induces pure point spectrum and shifting the operator accordingly). It follows that for large n, we have n (v) = 1. But the computation of n (v) is only dependent on v( j) for j < N for some N ∈ N. Define v 0 ∈ V by v 0 ( j) = v( j) if j < N and v 0 ( j) = 0 otherwise. It follows that n (v 0 ) = 1. But since the potential has compact support, 1 / ∈ σ p (H v 0 ) and hence n (v 0 ) = 0, the required contradiction.
We now prove that C pp can be computed using a height two arithmetical tower. The first step is the following technical lemma, whose proof will also be used later when considering C ac . Lemma 5.3. Let a < b with a, b ∈ R and consider the decision problem Then there exists a height two arithmetical tower n 2 ,n 1 (with evaluation functions 1 ) for a,b,pp . Furthermore, the final limit is from below in the sense that n 2 (T ) := lim n 1 →∞ n 2 ,n 1 (T ) ≤ a,b,pp (T ).

Proof.
Step 1 of the proof of Theorem 3.2 yields a height two arithmetical tower j n 2 ,n 1 (T ) for the computation of μ T e j ,e j ,c ((a, b)). Note that the final limit is from above and using the fact that μ T e j ,e j ,c ({a, b}) = 0, we obtain a height two tower for μ T e j ,e j ,c ([a, b]). We can then use the height one tower for μ T e j ,e j ([a, b]) constructed in Sect. 2.2, denoted by This provides a height two arithmetical tower for μ T e j ,e j ,pp ([a, b]) with the final limit from below. Without loss of generality (by taking successive maxima), we can assume that these towers are non-decreasing in n 2 . Now set Then it is clear that the limit lim n 1 →∞ ϒ n 2 ,n 1 (T ) = ϒ n 2 (T ) exists. Furthermore, the monotonicity of a j,n 2 ,n 1 (T ) in n 2 implies that To convert this to a height two tower for the decision problem a,b,pp , that maps to the discrete space {0, 1}, we use the following trick. Consider the intervals J n 2 1 = [0, 1/n 2 ], and J n 2 2 = [2/n 2 , ∞). Let k(n 2 , n 1 ) ≤ n 1 be maximal such that ϒ n 2 ,n 1 (T ) ∈ J n 2 1 ∪ J n 2 2 . If no such k exists or ϒ n 2 ,k (T ) ∈ J n 2 1 then set n 2 ,n 1 (T ) = 0. Otherwise set n 2 ,n 1 (T ) = 1. These can be computed using finitely many arithmetic operations and comparisons using 1 . The point of the intervals J n 2 1 and J n 2 2 is that we can show lim n 1 →∞ n 2 ,n 1 (T ) = n 2 (T ) exists. This is because lim n 1 →∞ ϒ n 2 ,n 1 (T ) = ϒ n 2 (T ) exists and hence we cannot oscillate infinitely often between the separated intervals J n 2 1 and J n 2 2 . Now suppose that a,b,pp (T ) = 0, then lim n 1 →∞ n 2 ,n 1 (T ) = 0 and hence lim n 1 →∞ n 2 ,n 1 (T ) = 0 for all n 2 . Now suppose that a,b,pp (T ) = 1, then for large enough n 2 we must have that ϒ n 2 (T ) > 2/n 2 and hence n 2 (T ) = 1. Together, these prove the convergence and that n 2 (T ) ≤ a,b,pp (T ).
To construct a height two arithmetical tower for C pp , we will use Lemma 5.3 repeatedly. Let n 2 ,n 1 (·, I ) denote the height two tower constructed in the proof of Lemma 5.3 for the closed interval I (I = [a, b]), where without loss of generality by taking successive maxima in n 2 , we can assume that this tower is non-decreasing in n 2 (this is where we use convergence from below in the final limit in the statement of the lemma). For a given n 1 and n 2 , we construct n 2 ,n 1 (T ) as follows (we will use some basic terminology from graph theory).
Define the intervals I 0 n 2 ,n 1 , j = [ j, j + 1] for j = −n 2 , . . . , n 2 − 1 so that these cover the interval [−n 2 , n 2 ]. Now suppose that I k n 2 ,n 1 , j are defined for j = 1, . . . , r k (n 2 , n 1 , T ). Compute each n 2 ,n 1 (T, I k n 2 ,n 1 , j ) and if this is 1, bisect I k n 2 ,n 1 , j via its midpoint into two equal halves consisting of closed intervals. We then take all these bisected intervals and label them as I k+1 n 2 ,n 1 , j for j = 1, . . . , r k+1 (n 2 , n 1 , k, T ). This is repeated until we have no further bisections or the intervals I n 2 n 2 ,n 1 , j have been computed. By adding the interval [−n 2 , n 2 ] as a root with children I 0 n 2 ,n 1 , j , this creates a finite tree structure where a non-root interval I is a parent of two intervals precisely if those two intervals are formed from its bisection and n 2 ,n 1 (T, I ) = 1. We then prune this tree by discarding all leaves I which have n 2 ,n 1 (T, I ) = 0 to form the tree T n 2 ,n 1 (T ). Finally, we let n 2 ,n 1 (T ) be the union of all the leaves of T n 2 ,n 1 (T ). Clearly this can be computed using finitely many arithmetic operations and comparisons using 1 . The construction is shown visually in Fig. 3.
In the above construction, the number of intervals considered (including those not in the tree T n 2 ,n 1 (T )) for a fixed n 2 is n 2 2 n 2 +1 + 1 and hence independent of n 1 . It follows that T n 2 ,n 1 (T ) and n 2 ,n 1 (T ) are constant for large n 1 (due to the convergence of the n 2 ,n 1 (T, I ) in {0, 1}). We denote these limiting values by T n 2 (T ) and n 2 (T ) respectively and also denote the corresponding intervals in the construction at the m−th level of this limit by I m n 2 , j . Note also that if C pp (T ) = ∅ then n 2 (T ) = ∅. Now suppose that z ∈ C pp (T ), then there exists a sequence of nested intervals I m = I m n 2 ,a m,n 2 containing z for m = 0, . . . , n 2 , where these intervals are independent of n 2 . Fix m, then for large n 2 we must have that n 2 (T, I j ) = 1 for j = 1, . . . , m. It follows that I m has a descendent interval I n 2 ,m contained in n 2 (T ) and hence we must have dist(z, n 2 (T )) ≤ 2 −m .
Since m was arbitrary it follows that dist(z, n 2 (T )) converges to 0 as n 2 → ∞.
Conversely, suppose that z m j ∈ m j (T ) with m j → ∞, then we must show that all limit points of {z m j } lie in C pp (T ). Suppose this were false, then by taking a subsequence if necessary, we can assume that z m j → z and dist(z m j , C pp (T )) ≥ δ for some δ > 0. We claim that it is sufficient to prove that the maximum length of the leaves of T n 2 (T ) intersecting a fixed compact subset of R, converges to zero as n 2 → ∞. Suppose this has been shown, then z m j ∈ I m j for some leaf I m j of T m j (T ). It follows that I m j ∩ C pp (T ) = ∅ and I m j → 0. But this contradicts z m j being positively separated from C pp (T ). We are thus left with proving the claim regarding the lengths of leaves. Suppose this were false, then there exists a compact set K ⊂ R and leaves I j in T b j (T ) such that the lengths of I j do not converge to zero and I j intersect K . By taking subsequences if necessary, we can assume that the lengths of each I j are constant. Then by the compactness of K and taking subsequences if necessary again, we can assume that each of the I j are equal to a common interval I . It follows that b j (T, I ) = 1 but that b j (T, I 1 ) = b j (T, I 2 ) = 0 since I is a leaf, where I 1 and I 2 form the bisection of I . Taking b j → ∞, this implies that I ∩ C pp (T ) = ∅ but I 1 ∩ C pp (T ) = I 2 ∩ C pp (T ) = ∅ which is absurd. Hence we have shown the required contradiction, and we have finished the proof.

Proof for absolutely continuous spectra.
To prove the lower bound (that one limit will not suffice), our strategy will be to reduce a certain decision problem to the computation of  Proof We are done if we prove the result for f (n) = n +1 and α = 0. In this case 1 and 2 are equivalent so we can restrict the argument to 1 . Suppose for a contradiction that n is a height one tower of general algorithms solving { C ac , f,0 , 2 }. We will gain a contradiction by using the supposed tower to solve { , }.
Given {a j } ∈ , consider the operator H = H 0 + v, where the potential is of the following form: The above comments show that each of these is a general algorithm, and it is clear that it converges to ({a j }) as n → ∞, the required contradiction.
To construct the height two (arithmetical) tower for C ac , we require the following lemma. with monotonic convergence from below. Using Corollary 2.2 (and the now standard argument of Lipschitz continuity of the resolvent), we can compute approximations of a m,n 2 ,n 1 (T ) to accuracy 1/n 1 in finitely many arithmetic operations and comparisons. Call these approximations a m,n 2 ,n 1 (T ) and set ϒ n 2 ,n 1 (T ) = max 1≤ j≤n 2 a j,n 2 ,n 1 (T ).
The proof now follows that of Lemma 5.3 exactly.
Proof that { C ac , f,α , 1 } ∈ A 3 . This is exactly the same construction as in the above proof of the inclusion { C pp , f,α , 1 } ∈ A 3 . We simply replace the tower constructed in the proof of Lemma 5.3 by the tower constructed in the proof of Lemma 5.4.

Proof for singular continuous spectra.
We first prove the lower bound for the singular continuous spectrum via Theorem 3.4. Note that the impossibility result { C sc , f,α , 2 } / ∈ G 2 follows from the same argument that was used to show { C ac , f,α , 2 } / ∈ G 2 . To show that two limits will not suffice for f (n) − n ≥ √ 2n + 1/2, our strategy will be to again reduce a certain decision problem to the computation of  M , d ). In [12], a Baire category argument was used to prove that SCI( , ) G = 3 (where the evaluation functions consist of component-wise evaluations of the array {a i, j }).
Assume that the function f satisfies f (n) − n ≥ √ 2n + 1/2. The proof will use a direct sum construction. Given {a i, j } ∈ , consider the operators H j = H 0 +v ( j) , where the potential is of the following form: We choose the following bijection (where m lists the canonical basis in each Hilbert space): · · · · · · A straightforward computation shows that H ∈ f,0 . We also observe that if ({a i, j }) = 1 then [0, 4] ⊂ σ sc (H ), otherwise σ sc (H ) ∩ (0, 4) = ∅. Suppose for a contradiction that n 2 ,n 1 is a height two tower of general algorithms that solves { C sc , f,0 , 1 }. We will gain a contradiction by using the supposed height two tower to solve { , }. We now set where we use the convention dist(3, ∅) = 1. The comments above show that each of these is a general algorithm. Furthermore, the convergence of n 2 ,n 1 implies that lim n 2 →∞ lim n 1 →∞ We are not quite done since the convergence here takes place on the interval [0, 1] with the usual metric as opposed to {0, 1} with the discrete metric. To get round this, we use the following, now familiar, trick.
Finally, we will use the following lemma to prove that the singular continuous spectrum can be computed in three limits. Let a < b with a, b ∈ R and consider the decision problem   a,b,sc :
Once this is proven, we can use the same construction that was used to prove , but with an additional limit. Namely, we replace (n 2 , n 1 ) by (n 3 , n 2 ) in the proof and use the tower constructed in the proof of Lemma 5.4 instead of n 2 ,n 1 (T, I ) for an interval I . We still gain the required convergence, since the only change is an additional limit in the finite number of decision problems that decide the appropriate tree. We see that each successive limit converges, with the second from above and the final from below. By taking successive maxima, minima of our base algorithms, we can assume that the second and final limits are monotonic and that ϒ n 3 ,n 2 ,n 1 (T ) is monotonic in both n 2 and n 3 . Define ϒ n 3 ,n 2 (T ) = lim n 1 →∞ ϒ n 3 ,n 2 ,n 1 (T ), ϒ n 3 (T ) = lim n 2 →∞ ϒ n 3 ,n 2 (T ) and ϒ(T ) = lim n 3 →∞ ϒ n 3 (T ). Then ϒ(T ) is zero if a,b,sc (T ) = 0, otherwise it is a positive finite number.
With a slight change to the previous argument (the monotonicity in n 2 and n 3 is crucial for this to work), consider the intervals These can be computed using finitely many arithmetic operations and comparisons using 1 , and, as before, the first limit exists with Note also that the second and third sequential limits exist through the use of maxima and minima. Now suppose that a,b,sc (T ) = 0 and fix n 3 . Then for large n 2 , we must have that ϒ m,n 2 (T ) < 1/(2n 3 ) for all m ≤ n 3 due to the monotonic convergence of ϒ p as p → ∞. It follows in this case that The conclusion of the lemma now follows.

Numerical Examples
We now demonstrate the applicability of the new algorithms. In particular, these are the first algorithms that compute their respective spectral properties for the whole class f,α , and even for the restricted case of tridiagonal self-adjoint matrices. The algorithms are also implicitly parallelisable, allowing large scale computations. We focus on discrete operators and the numerical implementation of the algorithms for ODEs, PDEs and integral operators will be the focus of a future software package. 8 6.1. Jacobi matrices and orthogonal polynomials. For our first set of examples, we consider the natural link between the spectral measures of Jacobi operators and orthogonal polynomials on R. Let J be a Jacobi matrix with a j , b j ∈ R and a j > 0. In this case, under suitable conditions, the probability measure μ J := μ J e 1 ,e 1 is the probability measure associated with the orthonormal polynomials defined by 3 ) are also shown. As discussed in the text, the rates reflect the local smoothness properties of the density f α,β and the spectral measure that appears in the multiplicative version of the spectral theorem (see, for example, [45,137,143]). Classically, one usually first considers the measure and then constructs the orthogonal polynomials and the corresponding J . In some sense, the algorithms constructed in this paper, and in particular Sect. 4, compute the inverse problem (note that J ∈ n→n+1,0 ). In other words, we compute the measure μ J given the recurrence coefficients defining the orthogonal polynomials. This is a very difficult problem to study theoretically, and, until now, there has been no numerical method able to tackle the problem in this generality (see Sect. 1.7 for comments on methods that tackle compact perturbations of Toeplitz operators). To verify our method, we will consider problems where the measure μ J is known analytically.
We begin with the well-known class of Jacobi polynomials defined for α, β > −1 which have and measure on the interval [−1, 1] given by where N (α, β) is a normalising constant, ensuring the measure is a probability measure.
To assess the convergence of the algorithm in Sect. 4.2 that approximates the Radon-Nikodym derivative f α,β , in this section we will plot various errors as a function of , the distance from the points at which we compute the resolvents to the real axis. Figure 4 (left) shows a typical error plot for α = 0.7 and β = 0.3. We plot both the L 1 error (computed using a large number of discrete points), and the pointwise errors at −1, 0 and 1. The procedure of Theorem 2.1 is used to determine adaptively how large our (rectangular) matrix truncations should be for a given . We see that both the L 1 error and error at 0 appear to decrease as O( ), 9 whereas the errors at −1 and  O( min{2,1+α,1+β} ) after extrapolation. These rates for pointwise and L 1 errors reflect the local Hölder exponent and integrability of the density and its first derivative respectively, and can be proven using a Taylor series argument for general measures. 10 Moreover, we found that increased rates of convergence could be obtained (and again proven) locally near smoother parts of the measure by using further iterates of extrapolation. Note also that we took a uniform value over the whole interval. However, could just have easily been a function of the position x, allowing it to be smaller for points where the resolvent is estimated more accurately for a given matrix size.
To demonstrate the algorithm on unbounded operators, we next consider the class of generalised Laguerre polynomials for α > −1 which have and measure on the interval [0, ∞) given by Results are shown in Fig. 5 for α = 0.5, where we have plotted the (relative) L 1 error over the interval [0, 1], as well as pointwise errors at 0 and 1. Similar conclusions can be drawn as before. Pointwise errors are also shown for this example and the Jacobi operator, but now using the 10th iterate of Richardson extrapolation, in Fig. 6. The errors decay at the expected rates (also shown), with O( 10 ) convergence at smooth 99 and x = 0.01 for the Jacobi and Laguerre cases, respectively), the prefactor in front of the O( 10 ) term is larger, and smaller is needed before the expected rate kicks in. Finally, we demonstrate the computation of measures for a Jacobi operator with non-empty discrete spectrum. The Charlier polynomials are generated by for α > 0, and have measure where δ m denotes a Dirac measure located at the point m. Figure 7 shows plots of π K H (x + i ; T, e 1 ), e 1 for = 10 −7 , computed using an (n + 1) × n matrix with n = 1000. The peaks clearly coincide with the atoms of the measure. The difference between the peak values and the weight exp(−α)α m /m! was of the order 10 −13 for both examples. This demonstrates that an effective way to compute eigenvalues (particularly the challenging case of those in gaps of the essential spectrum, where spectral pollution occurs, or even those embedded in the essential spectrum) and projections onto eigenspaces may be to find local maxima or spikes of π K H (x + i ; T, e 1 ), e 1 . Such a routine would only require the solution of shifted linear systems (the resolvent), without diagonalisation, and could be executed rapidly in parallel.

A global collocation approach.
Typically, the size of the linear system needed to approximate the resolvent accurately (Theorem 2.1) grows as ↓ 0 and we approach the spectrum. It is therefore beneficial to increase the rate of convergence of approximating the measures as ↓ 0. One local (in terms of x ∈ R) approach, Richardson extrapolation, was used in Sect. 6.1. Here we briefly outline a different, more global approach. Suppose that we know the spectral measure of an operator T has support included in I ⊂ R and is absolutely continuous. Alternatively, we may analytically know, or be able to estimate, the singular part of the measure and subtract this from what follows. In this case, a natural way to approximate the Radon-Nikodym derivative is through a formal basis expansion where φ m are functions with support I whose Cauchy's transforms are easy to compute. To approximate the coefficients a m , we collocate in the complex plane as follows. Let C be a finite collection of complex points in the upper half-plane and truncate the approximation of ρ T x,y to M terms. To generate a linear system for {a m } M m=1 , we evaluate the Cauchy transform at points z ∈ C. The Cauchy transform satisfies which can be computed with error control using the results of Sect. 2.1. Call this approximation b x,y (z) and define Then for each z ∈ C, an approximate linear relation can be written as Evaluating this relation at ≥ M points in C gives rise to a linear system, which can be inverted in the least-squares sense for the approximation of the coefficients {a m } M m=1 . If x = y, and the basis functions are real, then the coefficients are real. Hence, in this case, taking real and imaginary parts of the linear system gives further linear relations without having to compute further resolvents.
If our basis functions satisfy recursion relations of the form

Hence, given the values of the integrals
we can compute φ m (z) for all m ∈ N. The integrals in (6.4) have simple forms for all the bases used in this paper. This method of computation of φ m is fast, meaning that the most expensive part of the collocation method is the computation of the right hand side of the linear system, that is, computing the resolvent. Figure 8 (left) revisits the Jacobi polynomials example in Sect. 6.1 and shows the collocation method in the case of using Chebyshev polynomials as the basis functions φ m (taking the first 1500). For collocation points, we took M Chebyshev nodes offset by an additional i to lie just above the interval [−1, 1] in the complex plane. Note that the rate of convergence is faster than O( min{2,1+α,1+β} ) achieved with the methods of Sect. 4 after extrapolation. There are at least two prices to pay, though. First, there is no current proof that collocation converges and, second, global regularity of the measure is needed. At the very least, we need to be able to subtract off the singular part of the measure. In practice, even if the measure is absolutely continuous, a large number of basis functions may be needed to capture the Radon-Nikodym derivative. Examples are given in Sect. 6.3, where collocation performs worse than the techniques of Sect. 4 due to either of these regularity conditions. Figure 8 (right) revisits the (generalised) Laguerre polynomials example in Sect. 6.1 and shows the collocation method in the case of using Laguerre functions (the polynomials multiplied by the square root of the weight function) as the basis functions with M = 1000. The collocation points were {1 2 /M 2 , 2 2 /M 2 , . . . , 1} + i. Again, this method converges with faster rates than that in Sect. 6.1.

CMV matrices and extensions to unitary operators.
We now demonstrate that the algorithms extend to the unitary case through use of the functions K D , namely, the convolution with the Poisson kernel of the unit disk. We will consider the class of CMV matrices (named after Cantero, Moral and Velázquez [28]) linked with orthogonal polynomials on the unit circle. A full discussion of this subject is beyond the scope of this paper, and we refer the reader to the monographs of Simon [130,131]. However, the background for this example is as follows. Given a probability measure μ on the unit circle T, whose support is not a finite set, we can define a system of orthogonal polynomials { n } ∞ n=0 by applying the Gram-Schmidt process to {1, z, z 2 , . . .}. Given a polynomial Q n (z) of degree n, we define the reversed polynomial Q * n (z) via Q * n (z) = z n Q n (1/z). Szegő's recurrence relation [139] is given by (6.5) where the α n are known as Verblunsky coefficients [148] and satisfy α j < 1. Verblunsky's theorem [147] sets up a one-to-one correspondence between μ and the coefficients {α j } ∞ j=0 . Define also This matrix is unitary and banded (and lies in U n→n+2,0 ). This last property does not hold for the so-called GGT representation [62,68,142], which has infinitely many non-zero entries in each row. The GGT representation uses the basis { n } ∞ n=0 , whereas the CMV representation obtains a basis via applying Gram-Schmidt to {1, z, z −1 , z 2 , z 2 , . . .}. The key result is that μ C := μ C e 1 ,e 1 is precisely the measure μ on the unit circle. Hence our new algorithms can be considered as a computational tool for the correspondence {α j } ∞ j=0 → μ, in much the same way as for orthogonal polynomials on the real line in Sect. 6.1.
The first example we consider are the Rogers-Szegő polynomials [117] given by where q ∈ (0, 1). In this case, which can be expressed in terms of the theta function. Figure 9 (left) shows the convergence of the new algorithm for various q and we see algebraic convergence as before, with rates O( ) and O( 2 ) before and after extrapolation respectively (here is such Radon-Nikodym derivatives for larger values of q (shown in the right of Fig. 9) have larger derivatives and hence larger pre-factor in front of these rates. We can use a similar collocation method as in Sect. 6.2, using the standard Fourier basis {e imθ } m∈Z . Note that the relevant Cauchy transforms can be computed explicitly using Cauchy's residue theorem. Collocation points inside and outside the unit disk are needed (collocating inside the unit disk causes the Cauchy transforms of the basis functions with negative m to vanish). This achieved machine precision using just 41 basis functions and collocating at distance = 0.1 from T. The next example we consider are the Geronimus polynomials [63], which have α j = a with |a| < 1. In this case, for |a + 1/2| ≤ 1/2, where θ ∈ [−π, π], b = 2arg(1 + a) and θ a = 2 sin −1 (|a|). Figure 10 (right) shows some typical examples, and note that these density functions are not smooth. When |a + 1/2| > 1/2, there is also a singular part of the measure (see [130] for the exact formula). In the cases shown in Fig. 10, the collocation method struggles to gain an accuracy beyond 10 −4 owing to discontinuity in the derivative of the Radon-Nikodym derivative and the algorithm based on convolutions is able to gain much more accurate results. Finally, a typical example is shown in Fig. 11 (left) for a = 0.8 (there is now a singular part located at θ = 0), where we have shown the output of the algorithm and the exact convolution with the Poisson kernel for r = 1/1.01, as well as the collocation method using 21 basis functions and collocation points {(1 ± 0.1)2π/21, (1 ± 0.1)2π · 2/21, . . . , (1 ± 0.1)2π }. Here, we see an exact agreement between the algorithm and convolution. Unsurprisingly in the presence of point spectra, the collocation method is unstable. Consistent with Theorem 4.2, the algorithm in Sect. 4.2 converges to the density over the portion of the spectrum, which is purely absolutely continuous. This is shown in the right of Fig. 11.

Fractional diffusion on a 2D quasicrystal.
In this example, we demonstrate computation of the functional calculus and consider operators acting on the graph of a Penrose  Note that the algorithm's output agrees almost perfectly with the exact convolution of the measure. We have also shown the output of the collocation method, which is unstable in the presence of the point mass. Right: Convergence of the algorithm from Theorem 4.2 over the portion of the spectrum which is purely absolutely continuous tile-the canonical model of a two-dimensional quasicrystal [47,140,146] (aperiodic crystals which typically have anomalous spectra/transport properties). Quasicrystals were discovered in 1982 by Shechtman [125] who was awarded the 2011 Nobel Prize in Chemistry for his discovery. Since then, quasicrystals have generated considerable interest due to their often exotic physical properties [135], with a vast literature on the physics and spectral properties of such aperiodic systems. 11 Unlike the one-dimensional case, little is known about the spectral properties of two-dimensional quasicrystals. A finite portion of the infinite tile is shown in Fig. 12, and we consider the natural graph whose vertices are the vertices of the tiling and edges correspond to the edges of the rhombi. Such a graph posses no periodic structure, and it is generically impossible to study its spectral properties analytically with current methods. The free Hamiltonian H 0 with summation over nearest neighbour sites (vertices). The first rigorous computation of the spectrum of H 0 (also with additional error control) was performed in [40]. We chose a natural ordering of the vertices as in Fig. 12, which leads to an operator H 0 acting on l 2 (N). The local bandwidth grows for this operator (our ordering is asymptotically optimal) and hence computation of powers H m 0 is infeasible for m 50, rendering polynomial approximations of the functional calculus intractable. In the above notation, H 0 ∈ f,0 with f (n) − n = O( √ n), and so this example provides a demonstration of the algorithm for a non-banded matrix. Throughout, we take u 0 = e 1 , though different initial conditions are handled in the same manner.
The ability to compute the functional calculus allows the solution of linear evolution equations. Given A ∈ N , a function F (continuous and bounded on σ (A)) and u 0 ∈ l 2 (N), consider the evolution equation We consider fractional diffusion governed by for α > 0. If α is an integer, then the solution can be represented via contour deformation as where γ is a closed contour looping once around the spectrum. Typically we took the rectangular contour shown in Fig. 12 and approximated the integral via Gaussian quadrature. This allows us to compute the solution with error control (we known the minimal distance between γ and σ (−H 0 ) so can bound the Lipschitz constant of the resolvent) and clearly, this holds for other functions F, holomorphic on a neighbourhood of σ (−H 0 ). Note that other methods, such as direct diagonalisation of finite square truncations or discrete time stepping (which is difficult if α / ∈ N), do not give error control and are slower. In fact, for this example, direct diagonalisation was impractical.
When α / ∈ N, we can still deform the contour, but not at 0 since 0 ∈ σ (−H 0 ). Hence we deform the contour to that shown in Fig. 12. For a discussion of contour methods applied to finite matrices (in the case that the spectrum is strictly positive), we refer the reader to [74]. Unfortunately, such methods cannot be applied here since 0 ∈ σ (−H 0 ). Moreover, 0 appears to not be an isolated point of the spectrum, meaning that dealing with this point separately is also impossible. Figure 13 shows the convergence of the algorithm for α = 1/2 and α = 1. For α = 1/2, error control is not given by our algorithm, so we computed an error by comparing to a "converged" solution using larger n. The l 2 error refers to the error in l 2 (N). The method converges algebraically for α = 1/2 (owing to the contour touching the spectrum at 0) but converges exponentially for α = 1 with similar convergence observed over a large range of times t. Figure 14 shows the magnitude (log scale) of the computed solution for various times. Note that the techniques presented here can be applied to any evolution equation of the form (6.8) on infinite-dimensional Hilbert spaces. They may also be useful for splitting methods/exponential integrators, which require fast computation of matrix/operator exponentials (see [83,100] and the references therein) and more generally in the field of infinite-dimensional ODE/PDE systems. For methods that compute general semigroups on infinite-dimensional separable Hilbert spaces with error control, see [35].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

A. The SCI Hierarchy: A Framework for Computation
Cornerstones in the SCI hierarchy are the definitions of a computational problem, a general algorithm and towers of algorithms. The basic objects in a computational problem are as follows: is some set, called the domain. The set is the set of objects that give rise to our computational problems. The problem function : → M is what we are interested in computing. Moreover, the set is the collection of functions that provide us with the information we are allowed to read.
Remark 7. Throughout the paper we have relaxed the condition that λ ∈ maps into C by considering evaluation functions consisting of intervals exhausting a set, piecewise constant functions of compact support etc. These seemingly more complicated objects can be effectively encoded by functions that map into C. For example, when considering the decomposition (3.1): of an open set U , we consider λ 1,m , λ 2,m with λ 1,m (U ) = a m and λ 2,m (U ) = b m . For the sake of clarity of presentation of the proofs, such encodings are used implicitly throughout the paper.
This leads to the following definition. The goal is to find algorithms that approximate the function . More generally, we need the concept of a tower of algorithms, which is needed to describe problems that need several limits in the computation. However, first one needs the definition of a general algorithm.
Note that the definition of a general algorithm is more general than the definition of a Turing machine or a Blum-Shub-Smale (BSS) machine. A general algorithm has no restrictions on the operations allowed. The only restriction is that it can only take a finite amount of information, though it is allowed to adaptively choose the finite amount of information it reads depending on the input. Condition (iii) assures that the algorithm consistently reads the information. Note that the purpose of such a general definition is to get strong lower bounds. In particular, the more general the definition is, the stronger a lower bound will be.
With a definition of a general algorithm, we can define the concept of towers of algorithms. In addition to a general tower of algorithms (defined above), we will focus on arithmetic towers. The definition of a general algorithm allows for strong lower bounds; however, to produce upper bounds we must add structure to the algorithm and towers of algorithms. An arithmetic tower allows for arithmetic operations and comparisons. is countable, we define the following: An arithmetic tower of algorithms of height k for { , , M, } is a tower of algorithms where the lowest functions = n k ,...,n 1 : → M satisfy the following: For each A ∈ the mapping (n k , . . . , n 1 ) → n k ,...,n 1 (A) = n k ,...,n 1 ({A f } f ∈ ) is recursive, the action of on A consists of only finitely many arithmetic operations and comparisons on {A f } f ∈ (A) , and n k ,...,n 1 (A) is a finite string of complex numbers that can be identified with an element in M.
Remark 8 (Recursiveness). By recursive we mean the following. If f (A) ∈ Q for all f ∈ , A ∈ , and is countable, then n k ,...,n 1 ({A f } f ∈ ) can be executed by a Turing machine [145], that takes (n k , . . . , n 1 ) as input, and that has an oracle tape consisting of can be executed by a Blum-Shub-Smale machine [17] that takes (n k , . . . , n 1 ), as input, and that has an oracle that can access any A f for f ∈ .
Given the above definitions, we can now define the Solvability Complexity Index: If there exists a tower { n } n∈N of type α and height one such that = n 1 for some n 1 < ∞, then we define SCI( , , M, ) α = 0. We may sometimes write SCI( , ) α to simplify notation when M and are obvious.
The definition of the SCI immediately induces the SCI hierarchy: In certain cases, the SCI hierarchy can also be refined with notions of error controlsee [12,33,36,37]. We also need the following result. operations. If C is not positive definite, then at some point a pivot is zero or negative, at this point the algorithm aborts. An alternative is the Cholesky decomposition. Although forming the lower triangular L ∈ C n×n (if it exists) such that C = L L * requires the use of radicals, the existence of L can be determined using finitely many arithmetic operations. This follows from the standard Cholesky algorithm, and we omit the details.
Finally, we also need the following result. This provides a height two arithmetic tower for and hence { , , } ∈ A 3 . To prove the lower bound, suppose for a contradiction that { n } is a sequence of general algorithms, using , such that lim n→∞ n ({a j }) = ({a j }).
We will construct a sequence {a j } such that n ({a j }) does not converge, providing the required contradiction.
Finally set a j = a k j for j ≤ N k . It is clear from (iii) of Definition A.2, that n k ({a j }) = n k ({a k j }) = (1 + (−1) k )/2 and this implies that n ({a j }) cannot converge, the required contradiction.

B. Partial Differential Operators: Proof of Theorem 1.1
Before stating and proving a more mathematically precise version of Theorem 1.1, we need to be precise about the computational problems involved. Recall that we consider L formally defined on L 2 (R d ) by with the class PDE consisting of self-adjoint L satisfying certain conditions given in Sect. 1.2.1. We take PDE to consist of We also implicitly assume that for any given L, the dimension d and integer N (order) are known. As well as this, we need to consider the pairs of vectors (members of L 2 (R d )) with which we compute the spectral properties. By the polarisation identity, we can consider equal functions of norm 1. We let V PDE consist of all f ∈ L 2 (R d ) of norm 1 such that (1) There exists a positive constant C and an integer D (both possibly unknown) such that | f (x)| ≤ C 1 + |x| 2D , almost everywhere on R d , that is, f is polynomially bounded. (2) The restrictions f | [−r,r ] d ∈ A r for all r > 0.
We then add to PDE the evaluation functions Similarly, from Theorem 4.1 we consider We assume that given F ∈ C b (R), we have access to piecewise constant functions F n supported in [−n, n] such that F − F n L ∞ ([−n,n]) ≤ n −1 . Finally, from Theorem 4.2 we consider We restrict this map to the quadruples (L , f, U ) such that U is strictly separated from supp(μ L f, f,sc ) ∪ supp(μ L f, f,pp ) and denote this subclass by PDE . We can now state the precise form of Theorem 1.1.
Theorem B.1 (Precise form of Theorem 1.1). Given the above set-up, In other words, we can construct a convergent sequence of arithmetic algorithms for each problem. Remark 9. The proof will make clear that we can assume different conditions on the operator L and function f . We simply choose an appropriate basis so that we can apply Theorem B.2. We can also extend this to scalar measures (through inner products) and also to the towers of algorithms (of height ≥ 2) used to compute the decompositions into pure point, absolutely continuous and singular continuous parts of measures and spectra.
To prove Theorem B.1, we need the following theorem, which is similar to the results of Sect. 2.1. Proof. The proof is similar to that of Theorem 2.1. Let (T, x, z) ∈ SA × S l 2 (N) × C\R.
We have that n = rank(P n ) = rank((T − z I )P n ) for large n since σ 1 (T − z I ) > 0 (recall that z / ∈ σ (T )). Hence we can define Suppose that n is large enough so that σ 1 (P n (T * − z I )(T − z I )P n ) > 1/n. Then n (T, x, z) is a least-squares solution of the optimisation problem argmin y (T −z I )P n y− x . The linear space span{e n : n ∈ N} forms a core of T and hence also of T − z I . It follows by invertibility of T − z I that given any > 0, there exists an m = m( ) and a y = y( ) with P m y = y such that It follows that for all n ≥ m, (T − z I ) n (T, x, z) − x ≤ (T − z I )y − x ≤ and hence that n (T, x, z) − R(z, T )x ≤ |Im(z)| .
Since > 0 was arbitrary, we see that n (T, x, z) converges to R(z, T )x. For n, m ∈ N, define the finite matrices B n = P n (T * − z I )(T − z I )P n , C m,n = P n (T * − z I )P m .
Given the evaluation functions inˆ , we have access to the entries of these matrices to asymptotic accuracy (i.e. for a given diverging subsequence a n , to precision O(2 −a n )). It follows that we have access to approximations of B n and C m,n denoted B n and C m,n respectively with B n − B n = O(n −1 ), C m,n − C m,n = O(n −1 ).
Recall that the O(·) notation also means independently of z and other parameters (though it may depend on T and x). Note that B −1 n can be computed using finitely many arithmetic operations and comparisons and the resolvent identity implies that From the proof of Proposition A.7 and a simple search routine, we can also compute σ 1 ( B n ) to accuracy n −2 using finitely many arithmetic operations and comparisons. Denote the approximation via τ n . We then define We know that x has norm 1 and hence we must have However, we can compute 1 − m j=1 |x (m) | 2 , which converges to zero as m → ∞. Proof of Theorem B.1. We choose an orthonormal basis of L 2 (R d ) so that we can carry over the results for l 2 (N) proven in this paper. In [37] it was shown that (the orthonormal basis of) tensor products of Hermite functions form a core for any L ∈ PDE . Namely for d = 1 we choose the Hermite functions