## Journal Articles

### Spectral Computations

 Computing spectral measures of self-adjoint operators Matthew Colbrook, Andrew Horning, Alex Townsend SIAM Review, 2021, 63(3), 489-524. (chosen as cover article) pdf journal link Using the resolvent operator, we develop an algorithm for computing smoothed approximations of spectral measures associated with self-adjoint operators. The algorithm can achieve arbitrarily high-orders of convergence in terms of a smoothing parameter for computing spectral measures of general differential, integral, and lattice operators. Explicit pointwise and $$L^p$$-error bounds are derived in terms of the local regularity of the measure. We provide numerical examples, including a partial differential operator, a magnetic tight-binding model of graphene, and compute one thousand eigenvalues of a Dirac operator to near machine precision without spectral pollution. The algorithm is publicly available in $$\texttt{SpecSolve}$$, which is a software package written in MATLAB. Computing spectral measures and spectral types Matthew Colbrook* Communications in Mathematical Physics, 2021, 384, 433-501 pdf journal link 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 $$L^2(\mathbb{R}^d)$$. 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. How to compute spectra with error control (IMA Lighthill-Thwaites Prize Winner 2021) (Smith-Knight/Rayleigh-Knight Prize Class I 2018) Matthew Colbrook*, Bogdan Roman, Anders Hansen Physical Review Letters, 2019, 122(25), 250201. (chosen as cover article) pdf journal link Computing the spectra of operators is a fundamental problem in the sciences, with wide-ranging applications in condensed-matter physics, quantum mechanics and chemistry, statistical mechanics, etc. While there are algorithms that in certain cases converge to the spectrum, no general procedure is known that (a) always converges, (b) provides bounds on the errors of approximation, and (c) provides approximate eigenvectors. This may lead to incorrect simulations. It has been an open problem since the 1950s to decide whether such reliable methods exist at all. We affirmatively resolve this question, and the algorithms provided are optimal, realizing the boundary of what digital computers can achieve. Moreover, they are easy to implement and parallelize, offer fundamental speed-ups, and allow problems that before, regardless of computing power, were out of reach. Results are demonstrated on difficult problems such as the spectra of quasicrystals and non-Hermitian phase transitions in optics. Follow up article in IMA Mathematics Today: Unscrambling the Infinite: Can we Compute Spectra? On the infinite-dimensional QR algorithm Matthew Colbrook*, Anders Hansen Numerische Mathematik, 2019, 143(1), 17-83. pdf journal link Spectral computations of infinite-dimensional operators are notoriously difficult, yet ubiquitous in the sciences. Indeed, despite more than half a century of research, it is still unknown which classes of operators allow for the computation of spectra and eigenvectors with convergence rates and error control. Recent progress in classifying the difficulty of spectral problems into complexity hierarchies has revealed that the most difficult spectral problems are so hard that one needs three limits in the computation, and no convergence rates nor error control is possible. This begs the question: which classes of operators allow for computations with convergence rates and error control? In this paper, we address this basic question, and the algorithm used is an infinite-dimensional version of the QR algorithm. Indeed, we generalise the QR algorithm to infinite-dimensional operators. We prove that not only is the algorithm executable on a finite machine, but one can also recover the extremal parts of the spectrum and corresponding eigenvectors, with convergence rates and error control. This allows for new classification results in the hierarchy of computational problems that existing algorithms have not been able to capture. The algorithm and convergence theorems are demonstrated on a wealth of examples with comparisons to standard approaches (that are notorious for providing false solutions). We also find that in some cases the IQR algorithm performs better than predicted by theory and make conjectures for future study. Pseudoergodic operators and periodic boundary conditions Matthew Colbrook* Mathematics of Computation, 2020, 89(322), 737-766. pdf journal link There is an increasing literature on random non-self-adjoint infinite matrices with motivations ranging from condensed matter physics to neural networks. Many of these operators fall into the class of "pseudoergodic" operators, which allows the elimination of probabilistic arguments when studying spectral properties. Parallel to this is the increased awareness that spectral properties of non-self-adjoint operators, in particular stability, may be better captured via the notion of pseudospectra as opposed to spectra. Although it is well known that the finite section method applied to these matrices does not converge to the spectrum, it is often found in practice that the pseudospectrum behaves better with appropriate boundary conditions. We make this precise by giving a simple proof that the finite section method with periodic boundary conditions converges to the pseudospectrum of the full operator. Our results hold in any dimension (not just for banded bi-infinite matrices) and can be considered as a generalisation of the well-known classical result for banded Laurent operators and their circulant approximations. Furthermore, we numerically demonstrate a convergent algorithm for the pseudospectrum, including cases where periodic boundary conditions converge faster than the method of uneven sections. Finally, we show that the result carries over to pseudoergodic operators acting on $$l^p$$ spaces for $$p\in [1,\infty ]$$.

### Computational Fluid Mechanics

 A hybrid analytical-numerical method for solving advection-dispersion problems on a half-line Felipe de Barros, Matthew Colbrook, Athanassios Fokas International Journal of Heat and Mass Transfer, 2019, 139, 482-491. pdf journal link This paper employs the unified transform, also known as the Fokas method, to solve the advection-dispersion equation on the half-line. This method combines complex analysis with numerics. Compared to classical approaches used to solve linear partial differential equations (PDEs), the unified transform avoids the solution of ordinary differential equations and, more importantly, constructs an integral representation of the solution in the complex plane which is uniformly convergent at the boundaries. As a consequence, such solutions are well suited for numerical computations. Indeed, the numerical evaluation of the solution requires only the computation of a single contour integral involving an integrand which decays exponentially fast for large values of the integration variable. A novel contribution of this paper, with respect to the solution of linear evolution PDEs in general, and the implementation of the unified transform in particular, is the following: using the advection-dispersion equation as a generic example, it is shown that if the transforms of the given data can be computed analytically, then the unified transform yields a fast and accurate method that converges exponentially with the number of evaluations N yet only has complexity . Furthermore, if the transforms are computed numerically using $$\mathcal{O}(M)$$ evaluations, the unified transform gives rise to a method with complexity $$\mathcal{O}(MN)$$. Results are successfully compared to other existing solutions. Scaling laws of passive-scalar diffusion in the interstellar medium Matthew Colbrook*, Xiangcheng Ma, Philip Hopkins, Jonathan Squire Monthly Notices of the Royal Astronomical Society, 2017, 467(2), 2421-2429. pdf journal link Passive-scalar mixing (metals, molecules, etc.) in the turbulent interstellar medium (ISM) is critical for abundance patterns of stars and clusters, galaxy and star formation, and cooling from the circumgalactic medium. However, the fundamental scaling laws remain poorly understood in the highly supersonic, magnetized, shearing regime relevant for the ISM. We therefore study the full scaling laws governing passive-scalar transport in idealized simulations of supersonic turbulence. Using simple phenomenological arguments for the variation of diffusivity with scale based on Richardson diffusion, we propose a simple fractional diffusion equation to describe the turbulent advection of an initial passive scalar distribution. These predictions agree well with the measurements from simulations, and vary with turbulent Mach number in the expected manner, remaining valid even in the presence of a large-scale shear flow (e.g. rotation in a galactic disc). The evolution of the scalar distribution is not the same as obtained using simple, constant ‘effective diffusivity’ as in Smagorinsky models, because the scale dependence of turbulent transport means an initially Gaussian distribution quickly develops highly non-Gaussian tails. We also emphasize that these are mean scalings that apply only to ensemble behaviours (assuming many different, random scalar injection sites): individual Lagrangian ‘patches’ remain coherent (poorly mixed) and simply advect for a large number of turbulent flow-crossing times.

### Computation and Analysis of PDEs

 Kernel density estimation with linked boundary conditions Matthew Colbrook*, Zdravko Botev, Karsten Kuritz, Shev MacNamara Studies in Applied Mathematics, 2020, 145(3), 357-396. pdf journal link Kernel density estimation on a finite interval poses an outstanding challenge because of the well‐recognized bias at the boundaries of the interval. Motivated by an application in cancer research, we consider a boundary constraint linking the values of the unknown target density function at the boundaries. We provide a kernel density estimator (KDE) that successfully incorporates this linked boundary condition, leading to a non‐self‐adjoint diffusion process and expansions in nonseparable generalized eigenfunctions. The solution is rigorously analyzed through an integral representation given by the unified transform (or Fokas method). The new KDE possesses many desirable properties, such as consistency, asymptotically negligible bias at the boundaries, and an increased rate of approximation, as measured by the AMISE. We apply our method to the motivating example in biology and provide numerical experiments with synthetic data, including comparisons with state‐of‐the‐art KDEs (which currently cannot handle linked boundary constraints). Results suggest that the new method is fast and accurate. Furthermore, we demonstrate how to build statistical estimators of the boundary conditions satisfied by the target function without a priori knowledge. Our analysis can also be extended to more general boundary conditions that may be encountered in applications.

## Selected Submitted Articles/Preprints

 The foundations of spectral computations via the solvability complexity index hierarchy: Part I Matthew Colbrook*, Anders Hansen arXiv pdf On the computation of geometric features of spectra of linear operators on Hilbert spaces Matthew Colbrook* pdf Computing Spectra - On the solvability complexity index hierarchy and towers of algorithms Jonathan Ben-Artzi, Matthew Colbrook, Anders Hansen, Olavi Nevanlinna, Markus Seidel arXiv pdf Can stable and accurate neural networks be computed? - On the barriers of deep learning and Smale's 18th problem Matthew Colbrook*, Vegard Antun, Anders Hansen arXiv pdf Computing semigroups with error control Matthew Colbrook* pdf A contour method for time-fractional PDEs and an application to fractional viscoelastic beam equations Matthew Colbrook*, Lorna Ayton pdf Stratified sampling based compressed sensing for structured signals Theresa Loss*, Matthew Colbrook, Anders Hansen pdf Bulk Localised Transport States in Infinite and Finite Quasicrystals via Magnetic Aperiodicity Dean Johnstone*, Matthew Colbrook*, Anne Nielsen, Patrik Öhberg, Callum Duncan arXiv

## Thesis and Other Publications

 The Foundations of Infinite-Dimensional Spectral Computations Matthew Colbrook PhD thesis, University of Cambridge, 2020. pdf repository Spectral computations in infinite dimensions are ubiquitous in the sciences. However, their many applications and theoretical studies depend on computations which are infamously difficult. This thesis, therefore, addresses the broad question, "What is computationally possible within the field of spectral theory of separable Hilbert spaces?" The boundaries of what computers can achieve in computational spectral theory and mathematical physics are unknown, leaving many open questions that have been unsolved for decades. This thesis provides solutions to several such long-standing problems. To determine these boundaries, we use the Solvability Complexity Index (SCI) hierarchy, an idea which has its roots in Smale's comprehensive programme on the foundations of computational mathematics. The Smale programme led to a real-number counterpart of the Turing machine, yet left a substantial gap between theory and practice. The SCI hierarchy encompasses both these models and provides universal bounds on what is computationally possible. What makes spectral problems particularly delicate is that many of the problems can only be computed by using several limits, a phenomenon also shared in the foundations of polynomial root-finding as shown by McMullen. We develop and extend the SCI hierarchy to prove optimality of algorithms and construct a myriad of different methods for infinite-dimensional spectral problems, solving many computational spectral problems for the first time. For arguably almost any operator of applicable interest, we solve the long-standing computational spectral problem and construct algorithms that compute spectra with error control. This is done for partial differential operators with coefficients of locally bounded total variation and also for discrete infinite matrix operators. We also show how to compute spectral measures of normal operators (when the spectrum is a subset of a regular enough Jordan curve), including spectral measures of classes of self-adjoint operators with error control and the construction of high-order rational kernel methods. We classify the problems of computing measures, measure decompositions, types of spectra (pure point, absolutely continuous, singular continuous), functional calculus, and Radon–Nikodym derivatives in the SCI hierarchy. We construct algorithms for and classify; fractal dimensions of spectra, Lebesgue measures of spectra, spectral gaps, discrete spectra, eigenvalue multiplicities, capacity, different spectral radii and the problem of detecting algorithmic failure of previous methods (finite section method). The infinite-dimensional QR algorithm is also analysed, recovering extremal parts of spectra, corresponding eigenvectors, and invariant subspaces, with convergence rates and error control. Finally, we analyse pseudospectra of pseudoergodic operators (a generalisation of random operators) on vector-valued lp spaces. All of the algorithms developed in this thesis are sharp in the sense of the SCI hierarchy. In other words, we prove that they are optimal, realising the boundaries of what digital computers can achieve. They are also implementable and practical, and the majority are parallelisable. Extensive numerical examples are given throughout, demonstrating efficiency and tackling difficult problems taken from mathematics and also physical applications. In summary, this thesis allows scientists to rigorously and efficiently compute many spectral properties for the first time. The framework provided by this thesis also encompasses a vast number of areas in computational mathematics, including the classical problem of polynomial root-finding, as well as optimisation, neural networks, PDEs and computer-assisted proofs. This framework will be explored in the future work of the author within these settings. Unscrambling the Infinite: Can we Compute Spectra? Matthew Colbrook IMA Mathematics Today, 2021. pdf journal

## Publication Profiles and Code

Code for $$\texttt{SpecSolve}$$ can be found here.
Code for $$\texttt{FIRENETs}$$ can be found here.