Past Numerical Analysis Group Internal Seminar

Incomplete Cholesky (IC) factorizations have long been an important tool in the armoury of methods for the numerical solution of large sparse symmetric linear systems Ax = b. In this talk, I will explain the use of intermediate memory (memory used in the construction of the incomplete factorization but is subsequently discarded)  and show how it can significantly improve the performance of the resulting IC preconditioner. I will then focus on extending the approach to sparse symmetric indefinite systems in saddle-point form. A limited-memory signed IC factorization of the form LDLT is proposed, where the diagonal matrix D has entries +/-1. The main advantage of this approach is its simplicity as it avoids the use of numerical pivoting.  Instead, a global shift strategy is used to prevent breakdown and to improve performance. Numerical results illustrate the effectiveness of the signed incomplete Cholesky factorization as a preconditioner.

Ever wondered how the log function in your code is computed? This talk, which was prepared for the 400th anniversary of Napier's development of logarithms, discusses the computation of reciprocals, exponentials and logs, and also my own work on some special functions which are important in Monte Carlo simulation.

Abstract: We propose a term structure power price model that, in contrast to widely accepted no-arbitrage based approaches, accounts for the non-storable nature of power. It belongs to a class of equilibrium game theoretic models with players divided into producers and consumers. Consumers' goal is to maximize a mean-variance utility function subject to satisfying inelastic demand of their own clients (e.g households, businesses etc.) to whom they sell the power on. Producers, who own a portfolio of power plants each defined by a running fuel (e.g. gas, coal, oil...) and physical characteristics (e.g. efficiency, capacity, ramp up/down times, startup costs...), would, similarly, like to maximize a mean-variance utility function consisting of power, fuel, and emission prices subject to production constraints. Our goal is to determine the term structure of the power price at which production matches consumption. In this talk we outline that such a price exists and develop conditions under which it is also unique. Under condition of existence, we propose a tractable quadratic programming formulation for finding the equilibrium term structure of the power price. Numerical results show performance of the algorithm when modeling the whole system of UK power plants.

We study the behaviour of orthogonal polynomials on triangles and their coefficients in the context of spectral approximations of partial differential equations.  For spectral approximation we consider series expansions $u=\sum_{k=0}^{\infty} \hat{u}_k \phi_k$ in terms of orthogonal polynomials $\phi_k$. We show that for any function $u \in C^{\infty}$ the series expansion converges faster than with any polynomial order.  With these result we are able to employ the polynomials $\phi_k$ in the spectral difference method in order to solve hyperbolic conservation laws.

It is a well known fact that discontinuities can arise leading to oscillatory numerical solutions. We compare standard filtering and the super spectral vanishing viscosity methods, which uses exponential filters build from the differential operator of the respective orthogonal polynomials.  We will extend the spectral difference method for unstructured grids by using 
 classical orthogonal polynomials and exponential filters. Finally we consider some numerical test cases.

A fast method for computing eigenpairs of positive definite matrices using GPUs is presented. The method uses Chebyshev polynomial spectral transformations to map the desired eigenvalues of the original matrix $A$ to exterior eigenvalues of the transformed matrix $p(A)$, making them easily computable using existing Krylov methods. The construction of the transforming polynomial $p(z)$ can be done efficiently and only requires knowledge of the spectral radius of $A$. Computing $p(A)v$ can be done using only the action of $Av$. This requires no extra memory and is typically easy to parallelize. The method is implemented using the highly parallel GPU architecture and for specific problems, has a factor of 10 speedup over current GPU methods and a factor of 100 speedup over traditional shift and invert strategies on a CPU.

Partial differential equations defined on surfaces appear in various applications, for example image processing and reconstruction of non-planar images. In this talk, I will present a penalty method for evolution equations, based on an implicit representation of the surface. I will derive a simple equation in the surrounding space, formulated with an extension operator, and then show some analysis and applications of the method.

Several problems lead to the question of how well can a fine grid function be approximated by a coarse grid function, such as preconditioning in finite element methods or data compression and image processing. Particular challenges in answering this question arise when the functions may be only piecewise-continuous, or when the coarse space is not nested in the fine space. In this talk, we solve the problem by using a stable approximation from a space of globally smooth functions as an intermediate step, thereby allowing the use of known approximation results to establish the approximability by a coarse space. We outline the proof, which relies on techniques from the theory of discontinuous Galerkin methods and on the theory of Helmholtz decompositions. Finally, we present an application of our to nonoverlapping domain decomposition preconditioners for hp-version DGFEM.

Given a model dynamical system, a model of any measuring apparatus relating states to observations, and a prior assessment of uncertainty, the probability density of subsequent system states, conditioned upon the history of the observations, is of some practical interest.

When observations are made at discrete times, it is known that the evolving probability density is a solution of the Bayesian filtering equations. This talk will describe the difficulties in approximating the evolving probability density using a Gaussian mixture (i.e. a sum of Gaussian densities). In general this leads to a sequence of optimisation problems and related high-dimensional integrals. There are other problems too, related to the necessity of using a small number of densities in the mixture, the requirement to maintain sparsity of any matrices and the need to compute first and, somewhat disturbingly, second derivatives of the misfit between predictions and observations. Adjoint methods, Taylor expansions, Gaussian random fields and Newton’s method can be combined to, possibly, provide a solution. The approach is essentially a combination of filtering methods and '4-D Var’ methods and some recent progress will be described.

Everybody has heard of the Faraday cage effect, in which a wire mesh does a good job of blocking electric fields and electromagnetic waves. For example, the screen on the front of your microwave oven keeps the microwaves from getting out, while light with its smaller wavelength escapes so you can see your burrito.  Surely the mathematics of such a famous and useful phenomenon has been long ago worked out and written up in the physics books, right?

Well, maybe.   Dave Hewett and I have communicated with dozens of mathematicians, physicists, and engineers on this subject so far, and we've turned up amazingly little.   Everybody has a view of why the Faraday cage mathematics is obvious, and most of their views are different.  Feynman discusses the matter in his Lectures on Physics, but so far as we can tell, he gets it wrong. 

For the static case at least (the Laplace equation), Hewett and I have made good progress with numerical explorations based on  Mikhlin's method backed up by a theorem.   The effect seems to much weaker than we had imagined -- are we missing something?  For time-harmonic waves (the Helmholtz equation), our simulations lead to further puzzles.  We need advice!  Where in the world is the literature on this problem? 

In 1890, G. H. Bryan demonstrated that when a ringing wine glass rotates, the shape of the vibration pattern precesses, and this effect is the basis for a family of high-precision gyroscopes. Mathematically, the precession can be described in terms of a symmetry-breaking perturbation due to gyroscopic effects of a geometrically degenerate pair of vibration modes.  Unfortunately, current attempts to miniaturize these gyroscope designs are subject to fabrication imperfections that also break the device symmetry. In this talk, we describe how these devices work and our approach to accurate and efficient simulations of both ideal device designs and designs subject to fabrication imperfections.

State-space models are a very popular class of time series models which have found thousands of applications in engineering, robotics, tracking, vision,  econometrics etc. Except for linear and Gaussian models where the Kalman filter can be used, inference in non-linear non-Gaussian models is analytically intractable.  Particle methods are a class of flexible and easily parallelizable simulation-based algorithms which provide consistent approximations to these inference problems. The aim of this talk is to introduce particle methods and to present the most recent developments in this area.

This talk introduces a random linear model to investigate the memory bandwidth barrier effect on current shared memory computers. Based on the fact that floating-point operations can be hidden by implicit compiling techniques, the runtime for memory intensive applications can be modelled by memory reference time plus a random term. The random term due to cache conflicts, data reuse and other environmental factors is proportional to memory reference volume. Statistical techniques are used to quantify the random term and the runtime performance parameters. Numerical results based on thousands representative matrices from various applications are presented, compared, analysed and validated to confirm the proposed model. The model shows that a realistic and fair metric for performance of iterative methods and other memory intensive applications should consider the memory bandwidth capability and memory efficiency.

The Euler-Maclaurin formula is a quadrature rule based on corrections to the trapezoid rule using odd derivatives at the end-points of the function being integrated. It appears that no one has ever thought about a related function approximation that will give us the Euler-Maclaurin quadrature rule, i.e., just like we can derive Newton-Cotes quadrature by integrating polynomial approximations of the function, we investigate, what function approximation will integrate exactly to give the corresponding Euler-Maclaurin quadrature. It turns out, that the right function approximation is a combination of a trigonometric interpolant and a polynomial.

To make the method more practical, we also look at the closely related Newton-Gregory quadrature, which is very similar to the Euler-Maclaurin formula but instead of derivatives, uses finite differences. Following almost the same procedure, we find another mixed function approximation, derivative free, whose exact integration yields the Newton-Gregory quadrature rule.

We discuss a new collocation-type method for numerically solving partial differential equations (PDEs) on the sphere.  The method uses radial basis function (RBF) approximations in a partition of unity framework for approximating spatial derivatives on the sphere.  High-orders of accuracy are achieved for smooth solutions, while the overall computational cost of the method scales linearly with the number of unknowns.  The discussion will be primarily limited to the transport equation and results will be presented for a few well-known test cases.  We conclude with a preliminary application to the non-linear shallow water wave equations on a rotating sphere.

Compressed sensing and matrix completion are techniques by which simplicity in data can be exploited for more efficient data acquisition. For instance, if a matrix is known to be (approximately) low rank then it can be recovered from few of its entries. The design and analysis of computationally efficient algorithms for these problems has been extensively studies over the last 8 years. In this talk we present a new algorithm that balances low per iteration complexity with fast asymptotic convergence. This algorithm has been shown to have faster recovery time than any other known algorithm in the area, both for small scale problems and massively parallel GPU implementations. The new algorithm adapts the classical nonlinear conjugate gradient algorithm and shows the efficacy of a linear algebra perspective to compressed sensing and matrix completion.

When applied to an inequality constrained optimization problem, interior point methods generate iterates that belong to the interior of the set determined by the constraints, thus avoiding/ignoring the combinatorial aspect of the solution. This comes at the cost of difficulty in predicting the optimal active constraints that would enable termination.  We propose the use of controlled perturbations to address this challenge. Namely, in the context of linear programming with only nonnegativity constraints as the inequality constraints, we consider perturbing the nonnegativity constraints so as to enlarge the feasible set. Theoretically, we show that if the perturbations are chosen appropriately, the solution of the original problem lies on or close to the central path of the perturbed problem and that a primal-dual path-following algorithm applied to the perturbed problem is able to predict the optimal active-set of the original problem when the duality gap (for the perturbed problem) is not too small. Encouraging preliminary numerical experience is obtained when comparing the perturbed and unperturbed interior point algorithms' active-set predictions for the purpose of cross-over to simplex.

Networks arise pervasively in biology, physics, technology, social science, and myriad other areas. An ordinary network consists of a collection of entities (called nodes) that interact via edges. "Multilayer networks" are a more general representation that can be used when nodes are connected to each other via multiple types of edges or a network changes in time.  In this talk, I will discuss how to find dense sets of nodes called "communities" in multilayer networks and some applications of community structure to problems in neuroscience and finance.

Sobolev spaces are the standard framework in which to analyse weak (variational) formulations of PDEs or integral equations and their numerical solution (e.g. using the Finite Element Method or the Boundary Element Method). There are many different ways to define Sobolev spaces on a given domain, for example via integrability of weak derivatives, completions of spaces of smooth functions with respect to certain norms, or restriction from spaces defined on a larger domain. For smooth (e.g. Lipschitz) domains things many of these definitions coincide. But on rough (e.g. fractal) domains the picture is much more complicated. In this talk I'll try to give a flavour of the sort of interesting behaviour that can arise, and what implications this behaviour has for a "practical" example, namely acoustic wave scattering by fractal screens. 

Quantifying the uncertainty in computational simulations is one of the central challenges confronting the field of computational science and engineering today. The uncertainty quantification of inverse problems is neatly addressed in the Bayesian framework, where instead of seeking one unique minimiser of a regularised misfit functional, the entire posterior probability distribution is to be characterised. In this talk I review the deep connection between deterministic PDE-constrained optimisation techniques and Bayesian inference for inverse problems, discuss some recent advances made in the Bayesian viewpoint by adapting deterministic techniques, and mention directions for future research.

Composite dilation wavelets are affine systems which extend the notion of wavelets by incorporating a second set of dilations.  The addition of a second set of dilations allows the composite system to capture directional information in addition to time and frequency information.  We classify admissible dilation groups at two extremes: frequency localization through minimally supported frequency composite dilation wavelets and time localization through crystallographic Haar-type composite dilation wavelets. 

Convolution is widely-used and fundamental mathematical operation
in signal processing, statistics, and PDE theory.

Unfortunately the CONV() method in Chebfun for convolving two chebfun 
objects has long been one of the most disappointingly slow features of 
the project. In this talk we will present a new algorithm, which shows 
performance gains on the order of a factor 100.

The key components of the new algorithm are:
* a convolution theorem for Legendre polynomials 
* recurrence relations satisfied by spherical Bessel functions
* recent developments in fast Chebyshev-Legendre transforms [1]

Time-permitting, we shall end with an application from statistics,
using the fact that the probability distribution of the sum of two 
independent random variables is the convolution of their individual 
distributions.

[1] N. Hale and A. Townsend, "A fast, simple, and stable Chebyshev-
Legendre transform using an asymptotic formula”, SISC (to appear).