{ padding-top: 0px; }

Here are some examples of my research
(see papers for further details)


The Foundations of Infinite-Dimensional Spectral Computations

Spectral computations in infinite dimensions are generally very hard. For example, standard discretisations/truncations of self-adjoint operators suffer from spectral pollution in gaps of the essential spectrum. The aim of this work has been to assess which infinite-dimensional spectral computational problems are possible to solve (there are many which are not) and to develop efficient algorithms for such problems. This has led to the development of the Solvability Complexity Index (SCI) hierarchy - a computational tool that classifies computational problems and allows the proof that algorithms are optimal. For example, consider the problem of computing the spectrum, \(\sigma(A)\), of self-adjoint (bounded) operators \(A\) acting on \(l^2(\mathbb{N})\). There is no algorithm \(\Gamma_n\) which, given access to the matrix elements of such an \(A\), converges to \(\sigma(A)\) as \(n\rightarrow\infty\). However, it is possible to compute the spectrum of such operators using two limits (with convergence in the Hausdorff metric): \[\lim_{m\rightarrow\infty}\lim_{n\rightarrow\infty} \Gamma_{m,n}(A)=\sigma(A).\] Given a computational spectral problem, there are typically three tasks; determining the SCI classification of the problem, finding assumptions that lower the index, and devising optimal algorithms to solve the problem. I have classified a large number of problems in this hierarchy, such as computing spectra with error control, discrete/essential spectra, fractal dimensions/Lebesgue measure of spectra, spectral gaps, spectra of PDEs, spectral radii, spectral capacity and spectral measures/type of normal operators (as given by spectral theorem). I have also developed efficient resolvent based techniques for spectral problems.



Example 1: Computing spectra with error control. This algorithm solves the long-standing computational spectral problem. Consider, for example, the graphical Laplacian of a Penrose tile in Figure 1 (left), modelling a quasicrystal. One of the algorithms I have developed computes the spectrum of such operators (those with known off-diagonal matrix decay and controlled resolvent growth near the spectrum) with error control. Denoting the output of the algorithm as \(\Gamma_n(A)\), \(\Gamma_n(A)\rightarrow\sigma(A)\) in the Attouch-Wets topology (local uniform convergence - for bounded operators this is just the Hausdorff metric on compact subsets of \(\mathbb{C}\)) and the algorithm computes a function \(E(n,z)\) such that \[\mathrm{dist}(z,\sigma(A))\leq E(n,z),\quad \sup_{z\in\Gamma_n(A)}E(n.z)=E(n)\rightarrow 0.\] Errors and time taken for the algorithm for a given \(n\) (number of basis sites) are shown in Figure 1 (right), along with those for the finite section method (standard truncation method). The finite section method is much slower and suffers from spectral pollution. This is highlighted in Figure 2 which shows the output and errors of each method. These are taken from the paper "How to compute spectra with error control" that appeared on the front cover of Physical Review Letters. The new algorithm is also parallelisable and local, is proven to be optimal, and can be extended to partial differential operators with coefficients of locally bounded total variation.



Figure 1: Left: Portion of Penrose tile. Right: Errors/time for algorithm (blue) vs finite section (green = periodic BCs, red = open BCs).


Figure 2: Output of algorithm (blue) vs finite section (green = periodic BCs, red = open BCs).



Example 2: Computing spectral measures. Figure 3 shows the computed spectral measure (given by the projection-valued measure version of the spectral theorem) of a model of graphene subjected to a magnetic field. This is free from boundary effects and computed via resolvent techniques applied directly to the infinite dimensional operator. Again, the new algorithm is also parallelisable and local. The algorithm can be used for discrete/lattice operators, differential operators and integral operators. The algorithm can also be adapted to achieve arbitrarily high order of convergence (in terms of a smoothing parameter) with explicit bounds in terms of the local regularity of the measure. High order methods are implemented in the software package \(\texttt{SpecSolve}\) written with Andrew Horning and available here, with additional examples.



Figure 3: Spectral measure (Radon-Nikodym derivative) vs magnetic field strength computed using a fourth order kernel for sharpness (see log-scale).



Example 3: State-of-the-art computation of bound states of Dirac operators. Figure 4 shows the computation of bound states of a Dirac operator with Coulomb-type potential. This is an important problem in computational chemistry, yet is also very difficult due to the fact that the interval \((-1,1)\) is a gap of the essential spectrum and eigenvalues cluster geometrically near \(1\). Previous methods typically achieve a few digits of precision for the first few eigenvalues (and there is a whole literature devoted to ingenious methods to beat spectral pollution). In contrast, the computation of (smoothed approximations of) spectral measures allows thousands of eigenvalues to be computed rapidly to near machine precision.



Figure 4: Left: Rescaled smoothed spectral measure. The spikes correspond to discrete eigenvalues and the magnified section shows the extreme clustering (the dashed line shows the \(1000\)th eigenvalue \(E_{1000}\)). Right: The errors in the computed eigenvalues in terms of the smoothing parameter \(\epsilon\).



Stability and Accuracy: Computational Theory of Neural Networks

There is an increasing awareness that the high performance of neural networks may come at a high price - most high performing neural networks are severely unstable to adversarial attacks (well chosen tiny perturbations). This is a problem in applications such as medical diagnosis and fraud detection. A simple example for the AUTOMAP network, reported in Nature as a state-of-the-art network, is shown in Figure 5 where the instability test from this paper is shown. This shows severe instability.




Figure 5: Stability test for AUTOMAP. The top row shows the perturbed input (true image on left). The bottom row shows the output of the neural network.


In contrast, I have been working on developing provably stable neural networks for image problems. A stability test applied to these algorithms (which searches for the worst possible perturbations - i.e. those that cause the biggest difference in output for their size) is shown in Figure 6. This demonstrates stability.

Our recent paper explores these issues, providing fundamental barriers and describing conditions under which (stable) NNs with a given accuracy can be computed by an algorithm. We also introduce Fast Iterative REstarted NETworks (FIRENETs) which we prove are stable and also converge exponentially in the number of layers.




Figure 6: Stability test for new neural networks. The size of the perturbations are the same as those in Figure 8. For a fair comparison, the perturbations have been optimised so as to maximise the difference (instability) in the output of the new neural network. The new network is stable (this is proven).



Boundary-Based Spectral Methods

Spectral methods are a powerful way of solving many PDEs. However, typically they take a global approach and expand the solution in terms of basis functions. This can be inefficient if the solution has different intrinsic length scales and/or the domain is infinite or has complicated geometry. This project seeks to develop boundary-based spectral methods to solve elliptic PDEs in domains with polygonal boundaries (and even curved boundaries).

Example 1: Local Mathieu Function Expansions. A canonical problem in areas such as acoustic scattering is solving the Helmholtz equation in unbounded domains. The solution for the exterior of a single plate can be found using separation of variables \[{\phi}(\nu,\tau)=\sum_{m=1}^\infty a_m \mathrm{se}_{m}(\tau)\mathrm{Hse}_{m}(\nu),\] where \((\nu,\tau)\) are local elliptic coordinates, \(\mathrm{se}_m\) denote sine-elliptic functions and \(\mathrm{Hse}_{m}\) denote Mathieu-Hankel functions. The goal of this project has been to develop such formulae numerically (using collocation methods) for multiple plates and complicated boundary conditions (such as variable elasticity and even nonlinear corrections), and then to apply these techniques to physical problems of interest. As a fun example, Figure 7 shows a scattering problem for multiple elastic plates of variable elasticity (plate vibrations shown).


Figure 7: Total field (real part of incident + scattered) for multiple elastic plates of variable elasticity responding to a quadrupole source. This figure was made for the new Waves in Complex Continua (Wavinar) seminar series.



Example 2: Unified Transform. Based on the unified transform, this subproject has developed a spectral space collocation method. Importantly, the method is boundary based, yet does not need a Green's function or the evaluation of singular integrals often encountered in boundary integral/element methods. Conceptually, the method is based on global relations which can be thought of as similar (though not the same) to a Fourier transform of boundary integral equations. For example, for the Helmholtz equation with frequency \(2\beta\) satisfied by \(q\) in some domain \(\mathcal{D}\), the relation reads \[\int_{\partial \mathcal{D}}e^{-i\beta(\lambda z+\frac{\bar{z}}{\lambda})}\left[q_{n}+\beta\left(\lambda \frac{dz}{ds}-\frac{1}{\lambda}\frac{d\bar{z}}{ds}\right)q\right]ds=0, \qquad \lambda\in\mathcal{C}(\mathcal{D}),\] where \(\mathcal{C}(\mathcal{D})\) denotes a suitable collocation set in the complex plane and \(q_n\) the normal derivative along the boundary. Through choosing suitable expansions of \(q\) and \(q_n\) along the boundary \(\partial\mathcal{D}\) and evaluating at suitable \(\lambda\), we can build up a linear system for the unknown coefficients. The basis can be local to a given side and chosen to capture singular/expected behaviour. This gives a fast and easy to use boundary spectral method and can be extended to separable PDEs. Other extensions I have worked on include curved geometries, separable PDEs, fast evaluation in domain interior and dealing with singular solutions.

As a simple example, consider scattering of a single (rigid) plate shown in Figure 8. Here the solution can be computed analytically and the type of singularity at the endpoints predicted (a square root type). Figure 8 also shows the jump in \(q\) computed across the plate (accurate to machine precision) whereas Figure 9 highlights the exponential convergence of the method across a broad range of frequencies. A spectral boundary integral method is also shown for comparison, which struggles due to the need to evaluate singular integrals.


Figure 8: Left: Total field (incident + scattered) for a single rigid plate (Neumann BCs). Right: Typical solutions.


Figure 9: Maximum relative error over plate. UT denotes unified transform, BIM denotes a semi-analytical boundary integral method designed for this type of problem.