John W. Pearson

Jemima M. Tabeart, Selime Gürol, John W. Pearson et al.

Covariance matrices are central to data assimilation and inverse methods derived from statistical estimation theory. Previous work has considered the application of an all-at-once diffusion-based representation of a covariance matrix operator in order to exploit inherent parallelism in the underlying problem. In this paper, we provide practical methods to apply block $α$-circulant preconditioners to the all-at-once system for the case where the main diffusion operation matrix cannot be readily diagonalized using a discrete Fourier transform. Our new framework applies the block $α$-circulant preconditioner approximately by solving an inner block diagonal problem via a choice of inner iterative approaches. Our first method applies Chebyshev semi-iteration to a symmetric positive definite matrix, shifted by a complex scaling of the identity. We extend theoretical results for Chebyshev semi-iteration in the symmetric positive definite setting, to obtain computable bounds on the asymptotic convergence factor for each of the complex sub-problems. The second approach transforms the complex sub-problem into a (generalized) saddle point system with real coefficients. Numerical experiments reveal that in the case of unlimited computational resources, both methods can match the iteration counts of the `best-case' block $α$-circulant preconditioner. We also provide a practical adaptation to the nested Chebyshev approach, which improves performance in the case of a limited computational budget. Using an appropriate choice of $α$ our new approaches are robust and efficient in terms of outer iterations and matrix--vector products.

Stefan Güttel, John W. Pearson

We explore the use of flexible GMRES as an outer wrapper for sketched GMRES. Building on a new bound for the residual of FGMRES in terms of the residual of the preconditioner, we derive a practical randomized solver that requires very little parameter tuning, while still being efficient and robust in the sense of generating non-increasing residual norms.

Patrick E. Farrell, John W. Pearson

We propose a new preconditioner for the Ohta--Kawasaki equation, a nonlocal Cahn--Hilliard equation that describes the evolution of diblock copolymer melts. We devise a computable approximation to the inverse of the Schur complement of the coupled second-order formulation via a matching strategy. The preconditioner achieves mesh independence: as the mesh is refined, the number of Krylov iterations required for its solution remains approximately constant. In addition, the preconditioner is robust with respect to the interfacial thickness parameter if a timestep criterion is satisfied. This enables the highly resolved finite element simulation of three-dimensional diblock copolymer melts with over one billion degrees of freedom.

Roland Herzog, John W. Pearson, Martin Stoll

Optimal transport problems pose many challenges when considering their numerical treatment. We investigate the solution of a PDE-constrained optimisation problem subject to a particular transport equation arising from the modelling of image metamorphosis. We present the nonlinear optimisation problem, and discuss the discretisation and treatment of the nonlinearity via a Gauss--Newton scheme. We then derive preconditioners that can be used to solve the linear systems at the heart of the (Gauss--)Newton method. With the optical flow in mind, we further propose the reduction of dimensionality by choosing a radial basis function discretisation that uses the centres of superpixels as the collocation points. Again, we derive suitable preconditioners that can be used for this formulation.

Sergey Dolgov, John W. Pearson

In this paper, we consider the fast numerical solution of an optimal control formulation of the Keller--Segel model for bacterial chemotaxis. Upon discretization, this problem requires the solution of huge-scale saddle point systems to guarantee accurate solutions. We consider the derivation of effective preconditioners for these matrix systems, which may be embedded within suitable iterative methods to accelerate their convergence. We also construct low-rank tensor-train techniques which enable us to present efficient and feasible algorithms for problems that are finely discretized in the space and time variables. Numerical results demonstrate that the number of preconditioned GMRES iterations depends mildly on the model parameters. Moreover, the low-rank solver makes the computing time and memory costs sublinear in the original problem size.

John W. Pearson, Margherita Porcelli, Martin Stoll

PDE-constrained optimization problems with control or state constraints are challenging from an analytical as well as numerical perspective. The combination of these constraints with a sparsity-promoting $\rm L^1$ term within the objective function requires sophisticated optimization methods. We propose the use of an Interior Point scheme applied to a smoothed reformulation of the discretized problem, and illustrate that such a scheme exhibits robust performance with respect to parameter changes. To increase the potency of this method we introduce fast and efficient preconditioners which enable us to solve problems from a number of PDE applications in low iteration numbers and CPU times, even when the parameters involved are altered dramatically.

Konstantin Rusch, John W. Pearson, Konstantinos C. Zygalakis

Recurrent neural networks (RNNs) have gained a great deal of attention in solving sequential learning problems. The learning of long-term dependencies, however, remains challenging due to the problem of a vanishing or exploding hidden states gradient. By exploring further the recently established connections between RNNs and dynamical systems we propose a novel RNN architecture, which we call a Hamiltonian recurrent neural network (Hamiltonian RNN), based on a symplectic discretization of an appropriately chosen Hamiltonian system. The key benefit of this approach is that the corresponding RNN inherits the favorable long time properties of the Hamiltonian system, which in turn allows us to control the hidden states gradient with a hyperparameter of the Hamiltonian RNN architecture. This enables us to handle sequential learning problems with arbitrary sequence lengths, since for a range of values of this hyperparameter the gradient neither vanishes nor explodes. Additionally, we provide a heuristic for the optimal choice of the hyperparameter, which we use in our numerical simulations to illustrate that the Hamiltonian RNN is able to outperform other state-of-the-art RNNs without the need of computationally intensive hyperparameter optimization.

John W. Pearson, Sheehan Olver, Mason A. Porter

The two most commonly used hypergeometric functions are the confluent hypergeometric function and the Gauss hypergeometric function. We review the available techniques for accurate, fast, and reliable computation of these two hypergeometric functions in different parameter and variable regimes. The methods that we investigate include Taylor and asymptotic series computations, Gauss-Jacobi quadrature, numerical solution of differential equations, recurrence relations, and others. We discuss the results of numerical experiments used to determine the best methods, in practice, for each parameter and variable regime considered. We provide 'roadmaps' with our recommendation for which methods should be used in each situation.